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.


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.

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.

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.

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

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

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

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:

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:

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:

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