Introduction

This code will load in the combined scRNA seq data from external datasets combined with the iPSC-chondrocytes from the current study and fit a topic model (k=7) to the data. It then interprets the topics.

Load data and packages

The merged data were created by running the code located in Aggregation of external scRNA-seq data for topic model.

library(fastTopics)
library(Seurat)
library(Matrix)

#load data (stored in a seurat object)
merged_data <- readRDS("data/external_scRNA/merged_external_scRNA.rds")

#Extract raw count matrix from seurat object and get it in correct format for fastTopics
#need to fit the model to the count data (unintegrated)
raw_counts <- merged_data@assays$RNA@counts
#remove genes without any counts in droplets
raw_counts <- raw_counts[rowSums(raw_counts > 0) > 0,] 
#get into correct orientation (barcodes x features)
raw_counts <- t(raw_counts)
dim(raw_counts)

[1] 33218 36470

Use fastTopics functions to fit a topic model k=7 to the data

if (file.exists("output/topicModel_k=7.rds")) {
  fit <- readRDS("output/topicModel_k=7.rds")
} else {
  fit <- fit_poisson_nmf(raw_counts,k = 7,numiter = 150)
  saveRDS(fit, "output/topicModel_k=7.rds")
}

#compute weights and topics (rescale each of l and f to add up to 1)
l <- fit$L
f <- fit$F
weights <- sweep(l, MARGIN = 2, colSums(f), `*`)
scale <- rowSums(weights)
weights <- weights / scale
topics <- f / colSums(f) # add up to 1

Heatmap

Here we select a few marker genes to represent iPSCs, MSCs, Hepatocytes, and Chondrocytes and visualize the relative loadings of these marker genes in each of the 7 fitted topics using a scaled heatmap.

library(stringr)
library(pheatmap)
library(dummies)

dummies-1.5.6 provided by Decision Patterns

library(tidyverse)

Registered S3 method overwritten by 'cli':
  method     from    
  print.boxx spatstat

── Attaching packages ────────────────────────────────── tidyverse 1.3.0 ──

✓ ggplot2 3.3.3     ✓ purrr   0.3.4
✓ tibble  3.0.4     ✓ dplyr   1.0.2
✓ tidyr   1.1.2     ✓ forcats 0.4.0
✓ readr   1.3.1     

── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
x tidyr::expand() masks Matrix::expand()
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x tidyr::pack()   masks Matrix::pack()
x tidyr::unpack() masks Matrix::unpack()

#selected markers for different cell types
MSC_markers <- c("THY1", "NT5E", "ENG")
Chondrocyte_markers <- c("COL2A1", "ACAN", "SOX9", "SOX5", "SOX6", "COL9A1")
Hepatocyte_markers <- c("ALB")
iPSC_markers <- c("POU5F1", "SOX2", "NANOG")
markers <- c(iPSC_markers, Hepatocyte_markers, Chondrocyte_markers, MSC_markers)
#annotate the markers with the cell type they represent
markers_description <- data.frame(marker_type = c(rep("iPSC", 3), rep("Hepatocyte", 1), rep("Chondrocyte", 6), rep("MSC", 3)))
rownames(markers_description) <- markers

#plot heatmap of relative expression of marker genes in each topic
topics_markers <- topics[markers,]
pheatmap(topics_markers, cluster_cols = FALSE, cluster_rows = FALSE, annotation_row = markers_description, scale = "row")

Past versions of “plot correlations external-1.png”

Structure Plot

The structure plot represents the topic membership proportions of individual cells in each of the cell types using stacked bar plots.

#get labels of cells (cell type) and reorder to group them
sample <- as.data.frame(merged_data@meta.data$Cell.Type)
sample_labels <- sample %>% 
     dplyr::transmute(Cell.Type = stringr::word(`merged_data@meta.data$Cell.Type`, start = 1))
rownames(sample_labels) <- rownames(merged_data@meta.data)
sample <- sample_labels %>% 
     dplyr::arrange(Cell.Type)


set.seed(1)
topic_colors <- c("gold", "grey", "firebrick", "forestgreen", "royalblue", "turquoise", "salmon")
topics_order <- c(2, 3, 1, 7, 5, 4, 6)
rows_keep <- sort(c(sample(which(sample_labels$Cell.Type == "Hepatocyte"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC-MSC"), 800),
               which(sample_labels$Cell.Type == "iPSC-Chondrocyte"),
               sample(which(sample_labels$Cell.Type == "iPSC-Chondro_GAH"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC-Osteo"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d7"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d14"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d28"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d42"), 800),
               sample(which(sample_labels$Cell.Type == "chou_chondrocyte"), 800),
               sample(which(sample_labels$Cell.Type == "ji_chondrocyte"), 800)))
structure_plot <- structure_plot(select(poisson2multinom(fit),loadings = rows_keep),
                      grouping = factor(sample_labels[rows_keep,"Cell.Type"], 
                                        c("Hepatocyte", "iPSC", "iPSC-MSC", 
                                          "iPSC-Osteo", "iPSC-Chondro_GAH",
                                          "iPSC-Chondrocyte", "Wu_chondrogenic_pellet_d7", "Wu_chondrogenic_pellet_d14", "Wu_chondrogenic_pellet_d28", "Wu_chondrogenic_pellet_d42", "chou_chondrocyte", "ji_chondrocyte")),
                      topics = topics_order,
                      colors = topic_colors[topics_order],
                      perplexity = c(50),
                      n = 6043,gap = 100,num_threads = 4,verbose = FALSE)
print(structure_plot)

Past versions of “structure plot external-1.png”

Differential Expression analysis

First, calculate differential occurrence of genes in individual topics vs all other topics.

diff_count_topics <- diff_count_analysis(fit, raw_counts)

Fitting 36470 x 7 = 255290 univariate Poisson models.
Computing log-fold change statistics.
Stabilizing log-fold change estimates using adaptive shrinkage.

Then plot the results using a volcanoplot for each topic

library(ggrepel)
library(cowplot)
volcano_plot_with_highlighted_genes <- function (diff_count_res, k, 
                                                 genes, ...) {
  dat <- data.frame(beta  = diff_count_res$beta[genes,k],
                    y     = abs(diff_count_res$Z[genes,k]),
                    label = genes) 
  rows <- match(genes,rownames(diff_count_res$beta))
  rownames(diff_count_res$beta)[rows] <- ""
  rownames(diff_count_res$Z)[rows] <- ""
  return(volcano_plot(diff_count_res,k = k,
                      ggplot_call = volcano_plot_ggplot_call,...) +
         geom_text_repel(data = dat,
                         mapping = aes(x = beta,y = y,label = label),
                         inherit.aes = FALSE,color = "black",size = 5,
                         fontface = "italic",segment.color = "black",
                         segment.size = 0.25,
                         na.rm = TRUE))
}


# This is used by volcano_plot_with_highlighted_genes to create the
# volcano plot.
ggplot_call_for_volcano_plot <- function (dat, y.label, topic.label) {
  ggplot(dat,aes_string(x = "beta",y = "y",fill = "mean",label = "label")) +
    geom_point(color = "white",stroke = 0.3,shape = 21,na.rm = TRUE) +
    scale_y_continuous(trans = "sqrt",
      breaks = c(0,1,2,5,10,20,50,100,200,500,1e3,2e3,5e3,1e4,2e4,5e4)) +
    scale_fill_gradient2(low = "deepskyblue",mid = "gold",high = "orangered",
                         midpoint = mean(range(dat$mean))) +
    geom_text_repel(color = "gray",size = 5,fontface = "italic",
                    segment.color = "gray",segment.size = 0.25,na.rm = TRUE) +
    labs(x = "log-fold change (\u03b2)",y = y.label,fill = "log10 mean") +
    theme_cowplot(font_size = 18) +
    theme(plot.title = element_text(size = 18,face = "plain"))}

#plots
p1 <- volcano_plot_with_highlighted_genes(diff_count_res = diff_count_topics,
                                          k = "k2",
                                          genes = NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,5000)) +
     guides(fill = "none") +
     ggtitle("topic 2 (Hepatocytes)")
p2 <- volcano_plot_with_highlighted_genes(diff_count_topics,
                                          "k3",
                                          NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,225)) +
     guides(fill = "none") +
     ggtitle("topic 3 (iPSCs)")
p3 <- volcano_plot_with_highlighted_genes(diff_count_topics,
                                          "k1",
                                          NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,225)) +
     guides(fill = "none") +
     ggtitle("topic 1 (iPSC-MSCs)")
p4 <- volcano_plot_with_highlighted_genes(diff_count_topics,
                                          "k7",
                                          NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,225)) +
     guides(fill = "none") +
     ggtitle("topic 7 (iPSC-Chondrocytes)")
p5 <- volcano_plot_with_highlighted_genes(diff_count_topics,
                                          "k5",
                                          NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,225)) +
     guides(fill = "none") +
     ggtitle("topic 5 (Chondrocyte 1)")
p6 <- volcano_plot_with_highlighted_genes(diff_count_topics,
                                          "k4",
                                          NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,225)) +
     guides(fill = "none") +
     ggtitle("topic 4 (Chondrocyte 2)")
p7 <- volcano_plot_with_highlighted_genes(diff_count_topics,
                                          "k6",
                                          NA,
                                          label_above_quantile = 0.9995) +
     #ylim(c(0,225)) +
     guides(fill = "none") +
     ggtitle("topic 6 (Chondrocyte 3)")
plot_grid(p1,p2,p3,p4,p5,p6,p7,nrow = 3,ncol = 3)

Past versions of “DE plots-1.png”

For Supplement: Fit Topic models with k = 6 and k = 8 and plot a Structure plot in each case

if (file.exists("output/topicModel_k=6.rds")) {
  fit <- readRDS("output/topicModel_k=6.rds")
} else {
  fit <- fit_poisson_nmf(raw_counts,k = 6,numiter = 100)
  saveRDS(fit, "output/topicModel_k=6.rds")
}

#compute weights and topics (rescale each of l and f to add up to 1)
l <- fit$L
f <- fit$F
weights <- sweep(l, MARGIN = 2, colSums(f), `*`)
scale <- rowSums(weights)
weights <- weights / scale
topics <- f / colSums(f) # add up to 1


#get labels of cells (cell type) and reorder to group them
sample <- as.data.frame(merged_data@meta.data$Cell.Type)
sample_labels <- sample %>% 
     dplyr::transmute(Cell.Type = stringr::word(`merged_data@meta.data$Cell.Type`, start = 1))
rownames(sample_labels) <- rownames(merged_data@meta.data)
sample <- sample_labels %>% 
     dplyr::arrange(Cell.Type)

set.seed(1)
topic_colors <- c("turquoise", "firebrick", "grey", "salmon", "royalblue", "forestgreen")
topics_order <- c(3, 2, 4, 6, 5, 1)
rows_keep <- sort(c(sample(which(sample_labels$Cell.Type == "Hepatocyte"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC-MSC"), 800),
               which(sample_labels$Cell.Type == "iPSC-Chondrocyte"),
               sample(which(sample_labels$Cell.Type == "iPSC-Chondro_GAH"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC-Osteo"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d7"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d14"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d28"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d42"), 800),
               sample(which(sample_labels$Cell.Type == "chou_chondrocyte"), 800),
               sample(which(sample_labels$Cell.Type == "ji_chondrocyte"), 800)))
structure_plot <- structure_plot(select(poisson2multinom(fit),loadings = rows_keep),
                      grouping = factor(sample_labels[rows_keep,"Cell.Type"], 
                                        c("Hepatocyte", "iPSC", "iPSC-MSC", "iPSC-Osteo", 
                                          "iPSC-Chondro_GAH", "iPSC-Chondrocyte", "Wu_chondrogenic_pellet_d7", "Wu_chondrogenic_pellet_d14", "Wu_chondrogenic_pellet_d28", "Wu_chondrogenic_pellet_d42", "chou_chondrocyte", "ji_chondrocyte")),
                      topics = topics_order,
                      colors = topic_colors[topics_order],
                      perplexity = c(50),
                      n = 6043,gap = 100,num_threads = 4,verbose = FALSE)
print(structure_plot)

Past versions of “fit topic model k6-1.png”

if (file.exists("output/topicModel_k=8.rds")) {
  fit <- readRDS("output/topicModel_k=8.rds")
} else {
  fit <- fit_poisson_nmf(raw_counts,k = 8,numiter = 150)
  saveRDS(fit, "output/topicModel_k=8.rds")
}

#compute weights and topics (rescale each of l and f to add up to 1)
l <- fit$L
f <- fit$F
weights <- sweep(l, MARGIN = 2, colSums(f), `*`)
scale <- rowSums(weights)
weights <- weights / scale
topics <- f / colSums(f) # add up to 1


#get labels of cells (cell type) and reorder to group them
sample <- as.data.frame(merged_data@meta.data$Cell.Type)
sample_labels <- sample %>% 
     dplyr::transmute(Cell.Type = stringr::word(`merged_data@meta.data$Cell.Type`, start = 1))
rownames(sample_labels) <- rownames(merged_data@meta.data)
sample <- sample_labels %>% 
     dplyr::arrange(Cell.Type)

set.seed(1)
topic_colors <- c("royalblue", "turquoise","firebrick","gold","grey", "black", "salmon", "forestgreen")
topics_order <- c(5, 3, 4, 7, 6, 1, 8, 2)
rows_keep <- sort(c(sample(which(sample_labels$Cell.Type == "Hepatocyte"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC-MSC"), 800),
               which(sample_labels$Cell.Type == "iPSC-Chondrocyte"),
               sample(which(sample_labels$Cell.Type == "iPSC-Chondro_GAH"), 800),
               sample(which(sample_labels$Cell.Type == "iPSC-Osteo"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d7"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d14"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d28"), 800),
               sample(which(sample_labels$Cell.Type == "Wu_chondrogenic_pellet_d42"), 800),
               sample(which(sample_labels$Cell.Type == "chou_chondrocyte"), 800),
               sample(which(sample_labels$Cell.Type == "ji_chondrocyte"), 800)))
structure_plot <- structure_plot(select(poisson2multinom(fit),loadings = rows_keep),
                      grouping = factor(sample_labels[rows_keep,"Cell.Type"], 
                                        c("Hepatocyte", "iPSC", "iPSC-MSC", "iPSC-Osteo", 
                                          "iPSC-Chondro_GAH", "iPSC-Chondrocyte", "Wu_chondrogenic_pellet_d7", "Wu_chondrogenic_pellet_d14", "Wu_chondrogenic_pellet_d28", "Wu_chondrogenic_pellet_d42", "chou_chondrocyte", "ji_chondrocyte")),
                      topics = topics_order,
                      colors = topic_colors[topics_order],
                      perplexity = c(50),
                      n = 6043,gap = 100,num_threads = 4,verbose = FALSE)
print(structure_plot)

Past versions of “fit topic model k8-1.png”

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cowplot_1.1.0     ggrepel_0.9.0     forcats_0.4.0    
 [4] dplyr_1.0.2       purrr_0.3.4       readr_1.3.1      
 [7] tidyr_1.1.2       tibble_3.0.4      ggplot2_3.3.3    
[10] tidyverse_1.3.0   dummies_1.5.6     pheatmap_1.0.12  
[13] stringr_1.4.0     Matrix_1.2-18     Seurat_3.2.3     
[16] fastTopics_0.4-35

loaded via a namespace (and not attached):
  [1] readxl_1.3.1          backports_1.1.10      workflowr_1.6.2      
  [4] plyr_1.8.6            igraph_1.2.4.1        lazyeval_0.2.2       
  [7] splines_3.6.1         listenv_0.8.0         scattermore_0.7      
 [10] digest_0.6.27         invgamma_1.1          htmltools_0.5.0      
 [13] fansi_0.4.1           SQUAREM_2020.4        gdata_2.18.0         
 [16] magrittr_2.0.1        tensor_1.5            cluster_2.1.0        
 [19] ROCR_1.0-7            globals_0.12.5        modelr_0.1.8         
 [22] RcppParallel_5.0.2    matrixStats_0.57.0    MCMCpack_1.4-9       
 [25] prettyunits_1.1.1     colorspace_2.0-0      rvest_0.3.6          
 [28] rappdirs_0.3.1        haven_2.3.1           xfun_0.8             
 [31] crayon_1.3.4          jsonlite_1.7.2        spatstat_1.64-1      
 [34] spatstat.data_1.7-0   survival_2.44-1.1     zoo_1.8-8            
 [37] glue_1.4.2            polyclip_1.10-0       gtable_0.3.0         
 [40] MatrixModels_0.4-1    leiden_0.3.1          future.apply_1.3.0   
 [43] abind_1.4-5           SparseM_1.78          scales_1.1.1         
 [46] DBI_1.1.0             miniUI_0.1.1.1        Rcpp_1.0.5           
 [49] viridisLite_0.3.0     xtable_1.8-4          progress_1.2.2       
 [52] reticulate_1.16       rsvd_1.0.1            truncnorm_1.0-8      
 [55] htmlwidgets_1.5.2     httr_1.4.2            gplots_3.0.1.1       
 [58] RColorBrewer_1.1-2    ellipsis_0.3.1        ica_1.0-2            
 [61] farver_2.0.3          pkgconfig_2.0.3       dbplyr_1.4.2         
 [64] uwot_0.1.10           deldir_0.1-23         labeling_0.4.2       
 [67] tidyselect_1.1.0      rlang_0.4.10          reshape2_1.4.3       
 [70] later_1.1.0.1         cellranger_1.1.0      munsell_0.5.0        
 [73] tools_3.6.1           cli_2.2.0             generics_0.0.2       
 [76] broom_0.7.0           ggridges_0.5.1        evaluate_0.14        
 [79] yaml_2.2.1            goftest_1.2-2         mcmc_0.9-7           
 [82] npsurv_0.4-0          knitr_1.23            fs_1.3.1             
 [85] fitdistrplus_1.0-14   caTools_1.17.1.2      RANN_2.6.1           
 [88] pbapply_1.4-0         future_1.18.0         nlme_3.1-140         
 [91] whisker_0.3-2         mime_0.9              quantreg_5.73        
 [94] xml2_1.3.2            rstudioapi_0.13       compiler_3.6.1       
 [97] plotly_4.9.2.1        png_0.1-7             lsei_1.2-0           
[100] spatstat.utils_1.17-0 reprex_0.3.0          stringi_1.4.6        
[103] lattice_0.20-41       vctrs_0.3.6           pillar_1.4.7         
[106] lifecycle_0.2.0       lmtest_0.9-37         RcppAnnoy_0.0.18     
[109] data.table_1.13.0     bitops_1.0-6          irlba_2.3.3          
[112] conquer_1.0.2         httpuv_1.5.1          patchwork_1.1.0      
[115] R6_2.5.0              promises_1.1.1        KernSmooth_2.23-15   
[118] gridExtra_2.3         codetools_0.2-16      assertthat_0.2.1     
[121] MASS_7.3-52           gtools_3.8.1          rprojroot_2.0.2      
[124] withr_2.3.0           sctransform_0.3.2     mgcv_1.8-28          
[127] parallel_3.6.1        hms_0.5.3             quadprog_1.5-8       
[130] grid_3.6.1            rpart_4.1-15          coda_0.19-4          
[133] rmarkdown_1.13        ashr_2.2-47           Rtsne_0.15           
[136] git2r_0.26.1          mixsqp_0.3-43         lubridate_1.7.9      
[139] shiny_1.3.2