Visualize and interpret topics in PBMC data

Last updated: 2020-09-11

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: 1d0f457

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 1d0f457. 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:    analysis/figure/
    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/rds/
    Ignored:    output/pbmc-purified/fits-pbmc-purified.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	1d0f457	Peter Carbonetto	2020-09-11	Added some volcano plots for 68k data to yet_more_pbmc.R.
Rmd	8d8e2e9	Peter Carbonetto	2020-09-11	Improved names of clusters in PBMC purified data for plots_pbmc analysis.
Rmd	88d80e1	Peter Carbonetto	2020-09-11	Added subclustering of cluster A2 in plots_pbmc.
Rmd	d749c17	Peter Carbonetto	2020-09-07	Added notes to plots_pbmc; committed pbmc-clustering.RData output.
html	a0aeead	Peter Carbonetto	2020-09-07	Added subclustering of A1 cluster, and added step to save clustering results.
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”.

TO DO:

Try 68k data with a larger number of topics (e.g., k=8).
Run differential expression analysis using the topic model (cache the results).
Run the differential expression analysis using the clusters (cache these results, too).
Give better names to the clusters in the 68k data.

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 “CD34+” and “CD14+”). 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("X",n)
pc1  <- pca[,1]
pc2  <- pca[,2]
x[pc1 < 0 & pc2 < 0.4] <- "CD34+"
x[pc1 > 0.5 & pc2 < 0.15] <- "CD14+"
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, CD34+ and CD14+, account for most of the samples in cluster C. We also define a third subset, X—a “background cluster”—containing all the samples that were not assigned to CD34+ or CD14+.

table(x)
# x
# CD14+ CD34+     X 
#  2909  7822   871

Now we turn to cluster A. Within this cluster, there is a large subcluster, which we label as “NK”; the subset of samples that are not assigned to this cluster are labeled “A2”. (The NK 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] <- "NK"
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

rows <- which(samples_purified$cluster == "A2")
fit  <- select(poisson2multinom(fit_purified),loadings = rows)
p7   <- pca_plot(fit,k = 4)
p8   <- pca_hex_plot(fit,c("PC1","PC2"),bins = c(0,1,10,20,100,Inf))
plot_grid(p7,p8)

pca <- prcomp(fit$L)$x
n   <- nrow(fit$L)
x   <- rep("T",n)
pc1 <- pca[,1]
pc2 <- pca[,2]
x[pc1 < 0.2 & pc2 < -0.3] <- "CD8+"
samples_purified[rows,"cluster"] <- x
p9 <- pca_plot_with_labels(fit,c("PC1","PC2"),x) +
      labs(fill = "cluster")
print(p9)

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

samples_purified$cluster <- factor(samples_purified$cluster)
table(samples_purified$cluster)
# 
#     B CD14+ CD34+  CD8+    NK     T     X 
# 10439  2909  7822  1685  8352 62577   871

The structure plot summarizes the topic proportions in each of these 7 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 == "NK"),250),
               sample(which(samples_purified$cluster == "CD8+"),250),
               sample(which(samples_purified$cluster == "T"),1200),
               sample(which(samples_purified$cluster == "B"),250),
               sample(which(samples_purified$cluster == "CD34+"),250),
               sample(which(samples_purified$cluster == "CD14+"),200),
               sample(which(samples_purified$cluster == "X"),200)))
p10 <- 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,50,100,70,70,50,50),
                      gap = 40,num_threads = 4,verbose = FALSE)
# Perplexity automatically changed to 82 because original setting of 100 was too large for the number of samples (250)
print(p10)

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 NK, B, CD34+, CD14+). And, indeed, they align closely with their inclusion in the individual FACS-purified data sets:

with(samples_purified,table(celltype,cluster))
#                               cluster
# celltype                           B CD14+ CD34+  CD8+    NK     T     X
#   CD19+ B                      10073     0     0     0     0     3     9
#   CD14+ Monocyte                   8  2369     1     3     0    27   204
#   CD34+                          352   539  7740    15     4    28   554
#   CD4+ T Helper2                   0     0    16     3     1 11180    13
#   CD56+ NK                         0     1    17    50  8285    28     4
#   CD8+ Cytotoxic T                 0     0     0  1538    60  8558    53
#   CD4+/CD45RO+ Memory              0     0    19    60     0 10141     4
#   CD8+/CD45RA+ Naive Cytotoxic     3     0     0    13     1 11932     4
#   CD4+/CD45RA+/CD25- Naive T       1     0    25     2     1 10438    12
#   CD4+/CD25 T Reg                  2     0     4     1     0 10242    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 NK, B, CD34+, CD14+. 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 (“CD34+”). 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 X 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] <- "CD14+"
x[pc1 < -0.3]  <- "D"
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 B. 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] <- "B"
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    A3     B 
# 59294   522  3589

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 well by clusters. Nonetheless, here we (somewhat arbitrarily) designate two subsets, NK and A1c, in which all samples in these subsets are mostly explained by topic 3 and 4, respectively:

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] <- "NK"
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
a0aeead	Peter Carbonetto	2020-09-07
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   A1c    A3     B    B3     C CD14+     D    NK 
# 28485 20842   522  3589   179   165  4011   819  9967

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 == "NK"),500),
               sample(which(samples_68k$cluster == "A1c"),800),
               sample(which(samples_68k$cluster == "B"),500),
               sample(which(samples_68k$cluster == "A3"),300),
               sample(which(samples_68k$cluster == "CD14+"),500),
               sample(which(samples_68k$cluster == "D"),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)
# Perplexity automatically changed to 98 because original setting of 100 was too large for the number of samples (300)
# Perplexity automatically changed to 58 because original setting of 70 was too large for the number of samples (179)
# Perplexity automatically changed to 53 because original setting of 80 was too large for the number of samples (165)
print(p17)

Version	Author	Date
a0aeead	Peter Carbonetto	2020-09-07
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   A1c    A3     B    B3     C CD14+     D
#   CD14+ Monocyte                   6     0     1     1    11     2  2840     1
#   CD19+ B                        541  1442   310  3577     1     0     0    37
#   CD34+                           10     0    47     5     1   163    30    21
#   CD4+ T Helper2                  46    21    16     0     0     0     6     8
#   CD4+/CD25 T Reg               5237   935    11     0     0     0     2     0
#   CD4+/CD45RA+/CD25- Naive T     491  1372     4     1     1     0     1     3
#   CD4+/CD45RO+ Memory           2842   215     0     0     0     0     2     0
#   CD56+ NK                       776     0    24     0     5     0    12     1
#   CD8+ Cytotoxic T             14338  4329    74     1     2     0    24     0
#   CD8+/CD45RA+ Naive Cytotoxic  4155 12494     9     0     3     0     0     5
#   Dendritic                       43    34    26     4   155     0  1094   743
#                               cluster
# celltype                          NK
#   CD14+ Monocyte                   0
#   CD19+ B                          0
#   CD34+                            0
#   CD4+ T Helper2                   0
#   CD4+/CD25 T Reg                  2
#   CD4+/CD45RA+/CD25- Naive T       0
#   CD4+/CD45RO+ Memory              2
#   CD56+ NK                      7958
#   CD8+ Cytotoxic T              2005
#   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 (CD14+) and dendritic cells (D).
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.

Differential expression analysis

Add code and text here.

B-cells

Add code and text here.

Natural killer cells

Add code and text here.

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-175 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

Differential expression analysis

B-cells

Natural killer cells