Visualize and interpret topics in PBMC data

Last updated: 2020-09-07

Checks: 7 0

Knit directory: single-cell-topics/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2.9000). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(1)

The command set.seed(1) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: f06eda5

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version f06eda5. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    data/droplet.RData
    Ignored:    data/pbmc_68k.RData
    Ignored:    data/pbmc_purified.RData
    Ignored:    data/pulseseq.RData
    Ignored:    output/droplet/fits-droplet.RData
    Ignored:    output/droplet/rds/
    Ignored:    output/pbmc-68k/fits-pbmc-68k.RData
    Ignored:    output/pbmc-68k/postfit-pbmc-68k-scd-ex-k=6.RData
    Ignored:    output/pbmc-68k/rds/
    Ignored:    output/pbmc-purified/fits-pbmc-purified.RData
    Ignored:    output/pbmc-purified/postfit-pbmc-purified-scd-ex-k=6.RData
    Ignored:    output/pbmc-purified/rds/
    Ignored:    output/pulseseq/fits-pulseseq.RData
    Ignored:    output/pulseseq/rds/

Untracked files:
    Untracked:  analysis/pbmc-clustering.RData

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/plots_pbmc.Rmd) and HTML (docs/plots_pbmc.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	f06eda5	Peter Carbonetto	2020-09-07	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	87cfee1	Peter Carbonetto	2020-09-07	Added further clustering of cluster A1 to plots_pbmc.
html	daf30f2	Peter Carbonetto	2020-09-07	Added further clustering of cluster A1 to plots_pbmc.
Rmd	162977c	Peter Carbonetto	2020-09-07	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	3716962	Peter Carbonetto	2020-09-07	Build site.
Rmd	80fb4a2	Peter Carbonetto	2020-09-07	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	bbaea07	Peter Carbonetto	2020-09-07	Cleaned up the code in plots_pbmc.Rmd.
Rmd	d95bf0a	Peter Carbonetto	2020-09-06	Working on next steps in plots_pbmc.R analysis, in particular comparison of differential expression analysis in clusters and topics.
Rmd	cff02ac	Peter Carbonetto	2020-09-05	Working on volcano plots in plots_pbmc analysis.
Rmd	86fcfce	Peter Carbonetto	2020-08-31	Fixed bug in create_progress_plots.
Rmd	908efaa	Peter Carbonetto	2020-08-27	A couple small edits to my notes in plots_pbmc.Rmd.
html	11f4282	Peter Carbonetto	2020-08-27	Further refinement of notes in plots_pbmc.
Rmd	ea66c0b	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”, verbose = TRUE)
html	006c07f	Peter Carbonetto	2020-08-27	Improved the 68k structure plot in plots_pbmc a bit.
Rmd	8b37a44	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”, verbose = TRUE)
html	c516da1	Peter Carbonetto	2020-08-27	Added more notes/thoughts to plots_pbmc analysis.
Rmd	6959d21	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”, verbose = TRUE)
html	858a54b	Peter Carbonetto	2020-08-27	Fixed dimensions of some of the plots in plots_pbmc analysis.
Rmd	9625e81	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”, verbose = TRUE)
html	04a90b5	Peter Carbonetto	2020-08-27	Made several improvements to the analysis of the k=6 68k fit in
Rmd	070f612	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”, verbose = TRUE)
Rmd	5ee1a1d	Peter Carbonetto	2020-08-27	Made some refinements to the 68k pbmc analysis in plots_pbmc.Rmd.
html	25a24f7	Peter Carbonetto	2020-08-27	Build site.
Rmd	655dfae	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	208d263	Peter Carbonetto	2020-08-27	Added hex plots to plots_pbmc.
Rmd	8bfc908	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	d4f9ec7	Peter Carbonetto	2020-08-27	A couple small edits to the text in plots_pbmc.
html	30d2ef1	Peter Carbonetto	2020-08-27	Added hex plot to plots_pbmc analysis.
Rmd	6507d90	Peter Carbonetto	2020-08-27	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	3ab7da1	Peter Carbonetto	2020-08-25	A few minor edits.
Rmd	2defe6d	Peter Carbonetto	2020-08-25	Added crosstab plot to plots_tracheal_epithelium analysis.
html	e23f091	Peter Carbonetto	2020-08-25	Added a few more notes to plots_pbmc.
Rmd	0ca90f1	Peter Carbonetto	2020-08-25	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	94370c2	Peter Carbonetto	2020-08-25	Added notes adjoining results on 68k data in plots_pbmc.
Rmd	c6e754d	Peter Carbonetto	2020-08-25	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	399c597	Peter Carbonetto	2020-08-25	Added another few clusters to 68k fit, revised structure plot, and
Rmd	1a0c2d8	Peter Carbonetto	2020-08-25	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	eac2d23	Peter Carbonetto	2020-08-25	Revised the cross-tab plot, and added notes to plots_pbmc analysis.
html	2d156b8	Peter Carbonetto	2020-08-25	Revised the cross-tab plot, and added notes to plots_pbmc analysis.
Rmd	e8971d8	Peter Carbonetto	2020-08-25	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	9d8bb28	Peter Carbonetto	2020-08-25	Revised structure plot in plots_pbmc.
Rmd	27f6f51	Peter Carbonetto	2020-08-25	A couple small adjustments to the plotting settings in plots_pbmc.
html	abb846e	Peter Carbonetto	2020-08-25	Added dotplot to compare clusters vs. cell-types.
Rmd	45c3c32	Peter Carbonetto	2020-08-25	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	f53c86c	Peter Carbonetto	2020-08-24	Made a few small improvements to plots_pbmc.
Rmd	f7c157a	Peter Carbonetto	2020-08-24	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	e13c17e	Peter Carbonetto	2020-08-24	Removed reordering of 68k fit, and fixed downstream PCA analyses.
Rmd	01a0c8b	Peter Carbonetto	2020-08-24	Fixed colours in 68k structure plot.
html	a0cb7c6	Peter Carbonetto	2020-08-24	Added 68k fit structure plot to plots_pbmc analysis.
Rmd	5b2fb5c	Peter Carbonetto	2020-08-24	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	b6489db	Peter Carbonetto	2020-08-23	Re-built plots_pbmc with new PCA plots from 68k fit.
Rmd	89d98ed	Peter Carbonetto	2020-08-23	Removed basic_pca_plot in plots.R; working on PCA plots for 68k fit.
html	13ee038	Peter Carbonetto	2020-08-23	Revised notes in plots_pbmc.Rmd.
Rmd	e766da5	Peter Carbonetto	2020-08-23	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	dbd5882	Peter Carbonetto	2020-08-23	Added some explanatory text to plots_pbmc.Rmd.
html	97d7e86	Peter Carbonetto	2020-08-23	Fixed subclustering of cluster A in purified PBMC data.
Rmd	005b217	Peter Carbonetto	2020-08-23	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	59777e7	Peter Carbonetto	2020-08-22	Added table comparing Zheng et al (2017) labeling with clusters in
Rmd	9cd4a08	Peter Carbonetto	2020-08-22	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	c87ddf8	Peter Carbonetto	2020-08-22	Added PBMC purified structure plot with the new clustering (A1, A2, B,
Rmd	c6dd5f2	Peter Carbonetto	2020-08-22	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	ce421ae	Peter Carbonetto	2020-08-22	A few small revisions to plots_pbmc analysis.
Rmd	ddc6de4	Peter Carbonetto	2020-08-22	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	a406a2f	Peter Carbonetto	2020-08-22	Added back first PCA plot for 68k pbmc data.
Rmd	8daa131	Peter Carbonetto	2020-08-22	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	7900d17	Peter Carbonetto	2020-08-22	Revamping the process of defining clusters in plots_pbmc using PCA.
Rmd	14dd774	Peter Carbonetto	2020-08-22	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	944ce0c	Peter Carbonetto	2020-08-22	Working on new PCA plots for purified PBMC c3 subclusters.
html	b05232d	Peter Carbonetto	2020-08-22	Build site.
Rmd	d0e9f6e	Peter Carbonetto	2020-08-22	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	fbb0697	Peter Carbonetto	2020-08-21	Build site.
Rmd	cdca13e	Peter Carbonetto	2020-08-21	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	216027a	Peter Carbonetto	2020-08-21	Re-built plots_pbmc webpage with structure plots.
Rmd	98888b2	Peter Carbonetto	2020-08-21	Added structure plots to plots_pbmc.R.
html	9310993	Peter Carbonetto	2020-08-21	Revised the PCA plots in pbmc_plots.
Rmd	3a746ff	Peter Carbonetto	2020-08-21	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	1091dd0	Peter Carbonetto	2020-08-21	Added code for Ternary plots to plots_pbmc.
html	6d3d7e4	Peter Carbonetto	2020-08-20	Added 68k PCA plots to plots_pbmc.
Rmd	6ec82ca	Peter Carbonetto	2020-08-20	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	38f07a2	Peter Carbonetto	2020-08-20	A few small revisions to the plots_pbmc analysis.
Rmd	fd0316f	Peter Carbonetto	2020-08-20	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	606cd97	Peter Carbonetto	2020-08-20	Added basic PCA plots to plots_pbmc.
Rmd	1b41a60	Peter Carbonetto	2020-08-20	workflowr::wflow_publish(“plots_pbmc.Rmd”)
html	d6e5d39	Peter Carbonetto	2020-08-20	Added PCA plot with purified PBMC clustering to plots_pbmc analysis.
Rmd	99301a7	Peter Carbonetto	2020-08-20	workflowr::wflow_publish(“plots_pbmc.Rmd”)
Rmd	bf23ca0	Peter Carbonetto	2020-08-20	Added manual labeling of purified PBMC data to plots_pbmc analysis.
html	0c1b570	Peter Carbonetto	2020-08-20	First build on plots_pbmc page.
Rmd	eb7f776	Peter Carbonetto	2020-08-20	workflowr::wflow_publish(“plots_pbmc.Rmd”)

Here we examine and compare the topic modeling results for the two closely related data sets from Zheng et al (2017), the mixture of FACS-purified PBMC data sets, and the “unsorted” 68k PBMC data. The goal of this analysis is to illustrate how the topic models fitted to these data sets can be used to learn about structure in the data. In particular, we would like to identify clusters, and interpret clusters and topics as “cell types” or “gene expression programs”.

Load the packages used in the analysis below, as well as additional functions that will be used to generate some of the plots.

library(tools)
library(dplyr)
library(fastTopics)
library(ggplot2)
library(cowplot)
source("../code/plots.R")

Mixture of FACS-purified PBMC data

We begin with the mixture of FACS-purified PBMC data.

load("../data/pbmc_purified.RData")
samples_purified <- samples
rm(counts,samples,genes)

Load the \(k = 6\) Poisson NMF model fit.

fit_purified <-
  readRDS("../output/pbmc-purified/rds/fit-pbmc-purified-scd-ex-k=6.rds")$fit

Here, we explore the structure of the single-cell data as inferred by the topic model. Specifically, we use PCA to uncover structure in the estimated topic proportions of the multinomial topic model. Although PCA is simple, we will see that it works well, both for visualization and identifying clusters, and avoids the complications of the popular t-SNE and UMAP nonlinear dimensionality reduction methods. (Note that, since the topic proportions sum to 1, there are only 5 PCs to examine, not 6.)

fit <- poisson2multinom(fit_purified)
pca <- prcomp(fit$L)$x

Three large clusters are evident from first two PCs. We label the three large clusters as “A”, “B” and “C”.

Since there are so many samples, the scatterplot suffers from “overplotting”. So it also helpful to view this PC projection as a density plot.

n   <- nrow(pca)
x   <- rep("C",n)
pc1 <- pca[,"PC1"]
pc2 <- pca[,"PC2"]
x[pc1 + 0.2 > pc2] <- "A"
x[pc2 > 0.25] <- "B"
x[(pc1 + 0.4)^2 + (pc2 + 0.1)^2 < 0.07] <- "C"
samples_purified$cluster <- x
p1 <- pca_plot_with_labels(fit_purified,c("PC1","PC2"),
                           samples_purified$cluster) +
      labs(fill = "cluster")
p2 <- pca_hex_plot(fit_purified,c("PC1","PC2"))
plot_grid(p1,p2,rel_widths = c(9,10))

Version	Author	Date
30d2ef1	Peter Carbonetto	2020-08-27
7900d17	Peter Carbonetto	2020-08-22
38f07a2	Peter Carbonetto	2020-08-20

Most of the samples are in cluster A:

table(x)
# x
#     A     B     C 
# 72614 10439 11602

A small number of outlying data points do not seem to belong to any the three clusters, or they fall in between the clusters. For these data points, we assign them (rather arbitrarily) to one of the three clusters.

From these plots, there also also appears to be finer scale structure. For example, judging by the density plot, cluster A appears to split into two subclusters. We will examine this finer scale structure below.

Also note that other PCs beyond the first two may also sometimes reveal additional clustering, and we will see examples of this in the 68k PBMC data.

Within cluster C, there are two mostly well-defined subclusters (labeled “C1” and “C2”). There appear to be at least a couple other smaller, less well-defined subclusters, but in this analysis we focus on the largest, most obvious clusters.

rows <- which(samples_purified$cluster == "C")
fit  <- select(poisson2multinom(fit_purified),loadings = rows)
pca  <- prcomp(fit$L)$x
n    <- nrow(pca)
x    <- rep("C3",n)
pc1  <- pca[,1]
pc2  <- pca[,2]
x[pc1 < 0 & pc2 < 0.4] <- "C1"
x[pc1 > 0.5 & pc2 < 0.15] <- "C2"
samples_purified[rows,"cluster"] <- x
p3 <- pca_plot_with_labels(fit,c("PC1","PC2"),x) +
      labs(fill = "cluster")
p4 <- pca_hex_plot(fit,c("PC1","PC2"),bins = c(0,1,5,10,100,Inf))
plot_grid(p3,p4,rel_widths = c(9,10))

Version	Author	Date
208d263	Peter Carbonetto	2020-08-27
7900d17	Peter Carbonetto	2020-08-22

The two subclusters, C1 and C2, account for most of the samples in cluster C. We also define a third subset, C3—a “background cluster”—containing all the samples that were not assigned to C1 or C2.

table(x)
# x
#   C1   C2   C3 
# 7822 2909  871

Now we turn to cluster A. Within this cluster, there is a large subcluster, which we label as “A1”; the subset of samples that are not assigned to this cluster are labeled “A2”. (The A1 subcluster is much less distinct than the other clusters we have seen so far, and may not show up clearly in this scatterplot—it is more apparent from the density plot.) Otherwise, there is no obvious additional clustering of the samples within cluster A.

rows <- which(samples_purified$cluster == "A")
fit  <- select(poisson2multinom(fit_purified),loadings = rows)
pca  <- prcomp(fit$L)$x
n    <- nrow(fit$L)
x    <- rep("A2",n)
pc1  <- pca[,1]
pc2  <- pca[,2]
x[pc1 > 0.55 - pc2 | pc1 > 0.65] <- "A1"
samples_purified[rows,"cluster"] <- x
p5 <- pca_plot_with_labels(fit,c("PC1","PC2"),x) +
      labs(fill = "cluster")
p6 <- pca_hex_plot(fit,c("PC1","PC2"))
plot_grid(p5,p6,rel_widths = c(9,10))

Version	Author	Date
208d263	Peter Carbonetto	2020-08-27
97d7e86	Peter Carbonetto	2020-08-23
7900d17	Peter Carbonetto	2020-08-22

In summary, we have subdivided the data into 6 subsets:

samples_purified$cluster <- factor(samples_purified$cluster)
table(samples_purified$cluster)
# 
#    A1    A2     B    C1    C2    C3 
#  8352 64262 10439  7822  2909   871

We also inspected principal components individually in each of these 6 clusters and we did not find any of clear-cut examples of subclustering withing these clusters.

The structure plot summarizes the topic proportions in each of these 6 subsets:

set.seed(1)
pbmc_purified_topic_colors <- c("gold","forestgreen","dodgerblue",
                                "gray","greenyellow","magenta")
pbmc_purified_topics <- c(2,5,3,1,4,6)
rows <- sort(c(sample(which(samples_purified$cluster == "A1"),250),
               sample(which(samples_purified$cluster == "A2"),1200),
               sample(which(samples_purified$cluster == "B"),250),
               sample(which(samples_purified$cluster == "C1"),250),
               sample(which(samples_purified$cluster == "C2"),200),
               sample(which(samples_purified$cluster == "C3"),200)))
p7 <- structure_plot(select(poisson2multinom(fit_purified),loadings = rows),
                     grouping = samples_purified[rows,"cluster"],
                     topics = pbmc_purified_topics,
                     colors = pbmc_purified_topic_colors[pbmc_purified_topics],
                     n = Inf,perplexity = c(70,100,70,70,50,50),
                     gap = 40,num_threads = 4,verbose = FALSE)
print(p7)

Version	Author	Date
208d263	Peter Carbonetto	2020-08-27
eac2d23	Peter Carbonetto	2020-08-25
2d156b8	Peter Carbonetto	2020-08-25
abb846e	Peter Carbonetto	2020-08-25
f53c86c	Peter Carbonetto	2020-08-24
13ee038	Peter Carbonetto	2020-08-23
97d7e86	Peter Carbonetto	2020-08-23
59777e7	Peter Carbonetto	2020-08-22
c87ddf8	Peter Carbonetto	2020-08-22
7900d17	Peter Carbonetto	2020-08-22
fbb0697	Peter Carbonetto	2020-08-21
216027a	Peter Carbonetto	2020-08-21

Out of the 6 topics, 4 of them (\(k = 2, 3, 4, 5\)) align closely with clusters (labeled A1, B, C1, C2). And, indeed, they align closely with their inclusion in the individual FACS-purified data sets:

with(samples_purified,table(celltype,cluster))
#                               cluster
# celltype                          A1    A2     B    C1    C2    C3
#   CD19+ B                          0     3 10073     0     0     9
#   CD14+ Monocyte                   0    30     8     1  2369   204
#   CD34+                            4    43   352  7740   539   554
#   CD4+ T Helper2                   1 11183     0    16     0    13
#   CD56+ NK                      8285    78     0    17     1     4
#   CD8+ Cytotoxic T                60 10096     0     0     0    53
#   CD4+/CD45RO+ Memory              0 10201     0    19     0     4
#   CD8+/CD45RA+ Naive Cytotoxic     1 11945     3     0     0     4
#   CD4+/CD45RA+/CD25- Naive T       1 10440     1    25     0    12
#   CD4+/CD25 T Reg                  0 10243     2     4     0    14

Based on the above results, we make a few observations:

Because of their close correspondence, subsequent analysis of topics 2, 3, 4 and 5 should yield similar results to analyzing the clusters A1, B, C1, C2. For example, cluster B corresponds almost exactly to the B-cell data set. The largest cluster, cluster A2, is mostly comprised of T-cells.
Cluster A2—see also the PCA plot above—is an example where analyzing the most prevalent topics (\(k = 1, 6\)) will yield different results than a cluster-based analysis.
Many samples labeled as “CD34+” are not assigned to the CD34+ cluster (C1). This reflects the fact that this population was much less pure (45%) than the others, so we would expect some cells labeled as “CD34+” to not necessarily be CD34+ cells.
Cluster C3 is a heterogeneous cluster with a relatively small number of samples (790) that could potentially contain additional clusters of biological relevance, but will likely be more challenging to analyze and interpret than the other clusters, so we do not investigate this cluster further.

In summary, a cluster-based analysis and topic-based analysis should yield mostly similar results, except for the analysis of cluster A2, which should benefit from a topic-based analysis.

Unsorted 68k PBMC data

Next, we turn to the 68k data set. One feature of this data set is that it is not biased by the FACS purification, so we expect to observe a greater variety—or more continuous range—of cells states.

load("../data/pbmc_68k.RData")
samples_68k <- samples
rm(counts,samples,genes)

Load the \(k = 6\) Poisson NMF model fit.

fit_68k <- readRDS("../output/pbmc-68k/rds/fit-pbmc-68k-scd-ex-k=6.rds")$fit

Compute PCs from the topic proportions.

fit <- poisson2multinom(fit_68k)
pca <- prcomp(fit$L)$x

From the \(k = 6\) fit, we find least three distinct clusters in the projection onto PCs 3 and 4. We label these clusters “A”, “B” and “C”, as above, while cautioning that this labeling does not imply a connection with the purified PBMC clusters above.

n   <- nrow(pca)
x   <- rep("A",n)
pc3 <- pca[,"PC3"]
pc4 <- pca[,"PC4"]
x[pc4 < -0.12 | pc3/1.9 - 0.17 > pc4] <- "B"
x[pc4 < -0.75] <- "C"
samples_68k$cluster <- x
p8 <- pca_plot_with_labels(fit_68k,c("PC3","PC4"),x) +
      labs(fill = "cluster")
p9 <- pca_hex_plot(fit_68k,c("PC3","PC4"))
plot_grid(p8,p9,rel_widths = c(9,10))

Version	Author	Date
04a90b5	Peter Carbonetto	2020-08-27

The vast majority of the cells are in cluster A.

table(samples_68k$cluster)
# 
#     A     B     C 
# 63405  5009   165

The wide range in the sizes of these clusters is striking; the smallest cluster (C) is less than 1% the size of the largest (A). By contrast, community detection methods such as the Louvain algorithm are biased toward more uniformly sized clusters (this is a known limitation of community detection methods).

Examine the top two PCs in cluster B, we identify two large clusters, with the remaining assigned to a “background cluster”, B3.

rows <- which(samples_68k$cluster == "B")
fit  <- select(poisson2multinom(fit_68k),loadings = rows)
pca  <- prcomp(fit$L)$x
n    <- nrow(pca)
x    <- rep("B3",n)
pc1  <- pca[,"PC1"]
x[pc1 > -0.12] <- "B1"
x[pc1 < -0.3]  <- "B2"
samples_68k[rows,"cluster"] <- x
p10 <- pca_plot_with_labels(fit,c("PC1","PC2"),x) +
       labs(fill = "cluster")
p11 <- pca_hex_plot(fit,c("PC1","PC2"),bins = c(0,1,5,10,20,Inf))
plot_grid(p10,p11)

Version	Author	Date
858a54b	Peter Carbonetto	2020-08-27
04a90b5	Peter Carbonetto	2020-08-27

Cluster A subdivides into two large clusters, labeled as A1 and A2. The remaining samples are assigned to a subset labeled “A3”.

rows <- which(samples_68k$cluster == "A")
fit  <- select(poisson2multinom(fit_68k),loadings = rows)
pca  <- prcomp(fit$L)$x
n    <- nrow(pca)
x    <- rep("A3",n)
pc2  <- pca[,"PC2"]
pc3  <- pca[,"PC3"]
x[2.5*pc3 < 0.4 - pc2] <- "A1"
x[pc3 > 0.75 - pc2] <- "A2"
samples_68k[rows,"cluster"] <- x
p12 <- pca_plot_with_labels(fit,c("PC2","PC3"),x) +
       labs(fill = "cluster")
p13 <- pca_hex_plot(fit,c("PC2","PC3"),bins = c(0,1,5,10,100,Inf))
plot_grid(p12,p13)

Version	Author	Date
858a54b	Peter Carbonetto	2020-08-27
04a90b5	Peter Carbonetto	2020-08-27
f53c86c	Peter Carbonetto	2020-08-24
b6489db	Peter Carbonetto	2020-08-23

Within cluster A, the vast majority of the samples are assigned to the A1 subcluster:

table(x)
# x
#    A1    A2    A3 
# 59294  3589   522

Although we do not find additional clusters within the large A1 subset, the continuous structure is nonetheless quite interesting, and worth investigating.

rows <- which(samples_68k$cluster == "A1")
fit  <- select(poisson2multinom(fit_68k),loadings = rows)
p14 <- pca_plot(fit,k = 3:4)
print(p14)

Version	Author	Date
858a54b	Peter Carbonetto	2020-08-27
04a90b5	Peter Carbonetto	2020-08-27
399c597	Peter Carbonetto	2020-08-25

From these two plots, we observe that topics 3 and 4 exist on a continuous spectrum, but that mixtures of topics 3 and 4 are relatively rare. This is particularly evident from a density plot:

p15 <- pca_hex_plot(fit,c("PC1","PC2"),bins = c(0,1,10,20,100,Inf))
print(p15)

Version	Author	Date
858a54b	Peter Carbonetto	2020-08-27
04a90b5	Peter Carbonetto	2020-08-27

Topic 3 appears to characterize natural killer cells, and topic 4 has yet to be characterized, but may represent some subset of naive T-cells. This is an example with interesting substructure that can’t be captured by clusters. This is also an example of substructure that is much better captured by linear representations such as PCA as opposed to nonlinear dimensionality reduction methods such as t-SNE and UMAP.

Add text here.

fit  <- select(poisson2multinom(fit_68k),loadings = rows)
pca  <- prcomp(fit$L)$x
n    <- nrow(pca)
x    <- rep("A1a",n)
q3   <- fit$L[,"k3"]
q4   <- fit$L[,"k4"]
x[q3 > 0.5] <- "A1b"
x[q4 > 0.6] <- "A1c"
samples_68k[rows,"cluster"] <- x
p16 <- pca_plot_with_labels(fit,c("PC1","PC2"),x) +
  labs(fill = "cluster")
print(p16)

Version	Author	Date
87cfee1	Peter Carbonetto	2020-09-07
daf30f2	Peter Carbonetto	2020-09-07

In summary, we have subdivided these data into 9 subsets:

samples_68k$cluster <- factor(samples_68k$cluster)
table(samples_68k$cluster)
# 
#   A1a   A1b   A1c    A2    A3    B1    B2    B3     C 
# 28485  9967 20842  3589   522  4011   819   179   165

Again, the wide range in cluster sizes is striking.

The structure plot summarizes the topic proportions in each of these 9 subsets:

set.seed(1)
pbmc_68k_topic_colors <- c("yellow","lightskyblue","salmon",
                           "firebrick","royalblue","olivedrab")
pbmc_68k_topics <- c(2,5,1,3,4,6)
rows <- sort(c(sample(which(samples_68k$cluster == "A1a"),1000),
               sample(which(samples_68k$cluster == "A1b"),500),
               sample(which(samples_68k$cluster == "A1c"),800),
               sample(which(samples_68k$cluster == "A2"),500),
               sample(which(samples_68k$cluster == "A3"),300),
               sample(which(samples_68k$cluster == "B1"),500),
               sample(which(samples_68k$cluster == "B2"),300),
               which(samples_68k$cluster == "B3"),
               which(samples_68k$cluster == "C")))
p17 <- structure_plot(select(poisson2multinom(fit_68k),loadings = rows),
                      grouping = samples_68k[rows,"cluster"],
                      topics = pbmc_68k_topics,
                      colors = pbmc_68k_topic_colors[pbmc_68k_topics],
                      perplexity = c(100,100,100,100,70,80,50,50,50),
                      n = Inf,gap = 40,num_threads = 4,verbose = FALSE)
print(p17)

Version	Author	Date
87cfee1	Peter Carbonetto	2020-09-07
daf30f2	Peter Carbonetto	2020-09-07
006c07f	Peter Carbonetto	2020-08-27
04a90b5	Peter Carbonetto	2020-08-27
399c597	Peter Carbonetto	2020-08-25
abb846e	Peter Carbonetto	2020-08-25
f53c86c	Peter Carbonetto	2020-08-24
a0cb7c6	Peter Carbonetto	2020-08-24
7900d17	Peter Carbonetto	2020-08-22
fbb0697	Peter Carbonetto	2020-08-21
216027a	Peter Carbonetto	2020-08-21

These subsets do not align as closely with the cell-type labeling inferred by Zheng et al (2017). This is not surprising considering that the Zheng et al labeling is based on the FACS-purified data set.

with(samples_68k,table(celltype,cluster))
#                               cluster
# celltype                         A1a   A1b   A1c    A2    A3    B1    B2    B3
#   CD14+ Monocyte                   6     0     0     1     1  2840     1    11
#   CD19+ B                        541     0  1442  3577   310     0    37     1
#   CD34+                           10     0     0     5    47    30    21     1
#   CD4+ T Helper2                  46     0    21     0    16     6     8     0
#   CD4+/CD25 T Reg               5237     2   935     0    11     2     0     0
#   CD4+/CD45RA+/CD25- Naive T     491     0  1372     1     4     1     3     1
#   CD4+/CD45RO+ Memory           2842     2   215     0     0     2     0     0
#   CD56+ NK                       776  7958     0     0    24    12     1     5
#   CD8+ Cytotoxic T             14338  2005  4329     1    74    24     0     2
#   CD8+/CD45RA+ Naive Cytotoxic  4155     0 12494     0     9     0     5     3
#   Dendritic                       43     0    34     4    26  1094   743   155
#                               cluster
# celltype                           C
#   CD14+ Monocyte                   2
#   CD19+ B                          0
#   CD34+                          163
#   CD4+ T Helper2                   0
#   CD4+/CD25 T Reg                  0
#   CD4+/CD45RA+/CD25- Naive T       0
#   CD4+/CD45RO+ Memory              0
#   CD56+ NK                         0
#   CD8+ Cytotoxic T                 0
#   CD8+/CD45RA+ Naive Cytotoxic     0
#   Dendritic                        0

A few more notes about these results:

As in the purified fit, here we identify a B-cells cluster (B) and topic (5) that closely match. Other close matches include CD34+ cells (C), CD14+ monocytes (B1) and dendritic cells (B2).
As we found above, unlike the FACS-purified data, we don’t identify a clear-cut cluster for NK cells from the 68k data; the NK cells are mixed in with the T-cells (subset A1). NK cells will therefore only emerge only after subsequent analysis of topic 3.
The small number (~1%) in the A3 and B3 subsets aren’t distinct subpopulations per se but rather appear to be in some intermediate developmental state, and the mixture of topics in these subsets captures this.

In summary, the topics and clusters seem to offer very much complementary biological insights, although subsequent analysis is needed to determine what these insights are.

save(list = c("samples_purified","samples_68k"),
     file = "pbmc-clustering.RData")
resaveRdaFiles("pbmc-clustering.RData")

sessionInfo()
# R version 3.6.2 (2019-12-12)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS Catalina 10.15.6
# 
# Matrix products: default
# BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] tools     stats     graphics  grDevices utils     datasets  methods  
# [8] base     
# 
# other attached packages:
# [1] cowplot_1.0.0      ggplot2_3.3.0      fastTopics_0.3-165 dplyr_0.8.3       
# 
# loaded via a namespace (and not attached):
#  [1] ggrepel_0.9.0        Rcpp_1.0.5           lattice_0.20-38     
#  [4] tidyr_1.0.0          prettyunits_1.1.1    assertthat_0.2.1    
#  [7] zeallot_0.1.0        rprojroot_1.3-2      digest_0.6.23       
# [10] R6_2.4.1             backports_1.1.5      MatrixModels_0.4-1  
# [13] evaluate_0.14        coda_0.19-3          httr_1.4.1          
# [16] pillar_1.4.3         rlang_0.4.5          progress_1.2.2      
# [19] lazyeval_0.2.2       data.table_1.12.8    irlba_2.3.3         
# [22] SparseM_1.78         hexbin_1.28.0        whisker_0.4         
# [25] Matrix_1.2-18        rmarkdown_2.3        labeling_0.3        
# [28] Rtsne_0.15           stringr_1.4.0        htmlwidgets_1.5.1   
# [31] munsell_0.5.0        compiler_3.6.2       httpuv_1.5.2        
# [34] xfun_0.11            pkgconfig_2.0.3      mcmc_0.9-6          
# [37] htmltools_0.4.0      tidyselect_0.2.5     tibble_2.1.3        
# [40] workflowr_1.6.2.9000 quadprog_1.5-8       viridisLite_0.3.0   
# [43] crayon_1.3.4         withr_2.1.2          later_1.0.0         
# [46] MASS_7.3-51.4        grid_3.6.2           jsonlite_1.6        
# [49] gtable_0.3.0         lifecycle_0.1.0      git2r_0.26.1        
# [52] magrittr_1.5         scales_1.1.0         RcppParallel_4.4.2  
# [55] stringi_1.4.3        farver_2.0.1         fs_1.3.1            
# [58] promises_1.1.0       vctrs_0.2.1          glue_1.3.1          
# [61] purrr_0.3.3          hms_0.5.2            yaml_2.2.0          
# [64] colorspace_1.4-1     plotly_4.9.2         knitr_1.26          
# [67] quantreg_5.54        MCMCpack_1.4-5

Visualize and interpret topics in PBMC data

Peter Carbonetto

Mixture of FACS-purified PBMC data

Unsorted 68k PBMC data