Initial exploration of the pancreas data set

workflowr

Untracked files:
    Untracked:  data/Immune_ALL_human.h5ad
    Untracked:  data/human_pancreas_norm_complexBatch.h5ad
    Untracked:  data/pancreas.RData

The aim of this analysis is to take an initial look at the “pancreas” data set that was featured in in the Luecken et al 2022 benchmarking paper, and prepare the data in a convenient form for subsequent analyses in R.

In addition to being featured in the Luecken et al paper, it has been used in several papers on “data integration” methods for single-cell data (also known as “batch correction” or “harmonization” methods). See for example the MNN paper. The Supplementary Note in the Luecken et al paper has additional references.

See Supplementary Fig. 13, Supplementary Note 3 and Supplementary Data 7 of the Luecken et al paper for more information on this data set.

First, load the packages needed for this analysis. Note that MatrixExtra is also used in one of the steps below.

library(tools)
library(Matrix)
library(hdf5r)
library(rsvd)
library(Rtsne)
library(ggplot2)
library(cowplot)
library(Seurat)

Download the file “human_pancreas_norm_complexBatch.h5ad” from figshare and copy it to the “data” subdirectory of this git repository. Then load the count data, and encode them as a sparse matrix:

dat <- H5File$new("../data/human_pancreas_norm_complexBatch.h5ad",mode = "r")
counts <- dat[["layers"]][["counts"]][,]
counts <- t(counts)
counts <- as(counts,"CsparseMatrix")
sample_info <- data.frame(id          = dat[["obs"]][["_index"]][],
                          tech        = dat[["obs"]][["tech"]][],
                          celltype    = dat[["obs"]][["celltype"]][],
                          size_factor = dat[["obs"]][["size_factors"]][],
                          stringsAsFactors = FALSE)
sample_info <- transform(sample_info,
                         tech     = factor(tech),
                         celltype = factor(celltype))
levels(sample_info$tech)     <- dat[["obs"]][["__categories"]][["tech"]][]
levels(sample_info$celltype) <- dat[["obs"]][["__categories"]][["celltype"]][]
genes <- dat[["var"]][["_index"]][]
rownames(counts) <- sample_info$id
colnames(counts) <- genes

Note that some of the data are not actually counts, so perhaps calling this matrix “counts” is a bit misleading. Regardless, in some of our analyses we will model these data as counts.

Also note that in Luecken et al the counts were log-transformed, but here we taking the untransformed data.

The matrix has 16,382 rows (cells) and 19,093 columns (genes), and about 18% of the entries are nonzeros:

nrow(counts)
ncol(counts)
mean(counts > 0)
# [1] 16382
# [1] 19093
# [1] 0.1779012

The pancreas data are actually a combination of several scRNA-seq data sets that are from different sequencing technologies or were processed in different ways:

table(sample_info$tech)
# 
#     celseq    celseq2 fluidigmc1    inDrop1    inDrop2    inDrop3    inDrop4 
#       1004       2285        638       1937       1724       3605       1303 
#    smarter  smartseq2 
#       1492       2394

The cells were previously annotated by cell type:

table(sample_info$celltype)
# 
#             acinar activated_stellate              alpha               beta 
#               1669                464               5493               4169 
#              delta             ductal        endothelial            epsilon 
#               1055               2142                313                 32 
#              gamma         macrophage               mast quiescent_stellate 
#                699                 79                 42                193 
#            schwann             t_cell 
#                 25                  7

Some of the cell types occur in only a very small number of cells.

The “size factors” (here, total counts per cell) vary across a very wide range:

s <- rowSums(counts)
pdat <- data.frame(log_size_factor = log10(s))
ggplot(pdat,aes(log_size_factor)) +
  geom_histogram(bins = 64,col = "black",fill = "black") +
  labs(x = "log(size factor)") +
  theme_cowplot(font_size = 10)

Past versions of size-factors-1.png

Most genes are expressed in at least one cell:

p <- colSums(counts)/sum(s)
sum(p > 0)
# [1] 18771

The (relative) gene expression levels also vary across a very wide range:

p <- p[p > 0]
pdat <- data.frame(log_rel_expression_level = log10(p))
ggplot(pdat,aes(log_rel_expression_level)) +
  geom_histogram(bins = 64,col = "black",fill = "black") +
  labs(x = "log-expression level (relative)") +
  theme_cowplot(font_size = 10)

Past versions of gene-expression-levels-2-1.png

Let’s now generate a 2-d nonlinear embedding of the cells using t-SNE. First, transform the counts into “shifted log counts”:

a <- 1
s <- rowSums(counts)
s <- s/mean(s)
Y <- MatrixExtra::mapSparse(counts/(a*s),log1p)

Next, project the cells onto the top 50 PCs:

set.seed(1)
U <- rsvd(Y,k = 50)$u

Now run t-SNE on the 50 PCs:

tsne <- Rtsne(U,dims = 2,perplexity = 100,pca = FALSE,
              num_threads = 8,verbose = TRUE)
sample_info$tsne1 <- tsne$Y[,1]
sample_info$tsne2 <- tsne$Y[,2]

t-SNE with cells colored by cell-type:

tsne_colors <- rep(c("#E69F00","#56B4E9","#009E73","#F0E442",
                     "#0072B2","#D55E00","#CC79A7"),times = 2)
tsne_shapes <- rep(c(19,17),each = 7)
ggplot(sample_info,aes(x = tsne1,y = tsne2,color = celltype,
                       shape = celltype)) +
  geom_point(size = 1.5) +
  scale_color_manual(values = tsne_colors) +
  scale_shape_manual(values = tsne_shapes) +
  labs(x = "tSNE 1",y = "tSNE 2") + 
  theme_cowplot(font_size = 10)

Past versions of umap-vs-cell-type-1.png

t-SNE with cells colored by batch:

ggplot(sample_info,aes(x = tsne1,y = tsne2,color = tech,shape = tech)) +
  geom_point(size = 1.5) +
  scale_color_manual(values = tsne_colors) +
  scale_shape_manual(values = tsne_shapes) +
  labs(x = "tSNE 1",y = "tSNE 2") + 
  theme_cowplot(font_size = 10)

Past versions of umap-vs-batch-1.png

It is clear from these t-SNE plots that both batch and cell-type contribute to structure in the data.

For comparison, let’s run the default t-SNE pipeline in Seurat:

pancreas <- CreateSeuratObject(counts = t(counts),project = "pancreas",
                               meta.data = sample_info)
pancreas <- NormalizeData(pancreas)
pancreas <- ScaleData(pancreas)
pancreas <- FindVariableFeatures(pancreas)
pancreas <- RunPCA(pancreas,npcs = 50,features = VariableFeatures(pancreas))
pancreas <- RunTSNE(pancreas)

Seurat t-SNE with cells colored by cell-type:

DimPlot(pancreas,reduction = "tsne",group.by = "celltype",
        shape.by = "celltype",pt.size = 1.5) +
  scale_color_manual(values = tsne_colors) +
  scale_shape_manual(values = tsne_shapes) +
  theme_cowplot(font_size = 10) +
  labs(title = "")

Seurat t-SNE with cells colored by batch:

DimPlot(pancreas,reduction = "tsne",group.by = "tech",
        shape.by = "tech",pt.size = 1.5) +
  scale_color_manual(values = tsne_colors) +
  scale_shape_manual(values = tsne_shapes) +
  theme_cowplot(font_size = 10) +
  labs(title = "")

The default Seurat t-SNE shows does not show as much structure in the data. This plot does not seem to pick up much batch structure data; on the other hand, it also does not pick up some of the more subtyle cell types.

(Note that I had to spend some time customizing the plots; among other things, Seurat appears to be using the default ggplot color scheme which is terrible.)

Finally, save the data and t-SNE results to an .Rdata file for more convenient analysis in R:

save(list = c("sample_info","counts"),file = "pancreas.RData")
resaveRdaFiles("pancreas.RData")

Session information

sessionInfo()
# R version 4.3.3 (2024-02-29)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS Sonoma 14.6.1
# 
# Matrix products: default
# BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
# LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# time zone: America/Chicago
# tzcode source: internal
# 
# attached base packages:
# [1] tools     stats     graphics  grDevices utils     datasets  methods  
# [8] base     
# 
# other attached packages:
# [1] Seurat_5.0.3       SeuratObject_5.0.1 sp_2.1-4           cowplot_1.1.3     
# [5] ggplot2_3.5.0      Rtsne_0.17         rsvd_1.0.5         hdf5r_1.3.11      
# [9] Matrix_1.6-5      
# 
# loaded via a namespace (and not attached):
#   [1] RColorBrewer_1.1-3     jsonlite_1.8.8         magrittr_2.0.3        
#   [4] spatstat.utils_3.0-4   farver_2.1.1           rmarkdown_2.26        
#   [7] fs_1.6.3               vctrs_0.6.5            ROCR_1.0-11           
#  [10] spatstat.explore_3.2-7 htmltools_0.5.7        sass_0.4.8            
#  [13] sctransform_0.4.1      parallelly_1.37.1      KernSmooth_2.23-22    
#  [16] bslib_0.6.1            htmlwidgets_1.6.4      ica_1.0-3             
#  [19] plyr_1.8.9             plotly_4.10.4          zoo_1.8-12            
#  [22] cachem_1.0.8           whisker_0.4.1          igraph_2.0.3          
#  [25] mime_0.12              lifecycle_1.0.4        pkgconfig_2.0.3       
#  [28] R6_2.5.1               fastmap_1.1.1          fitdistrplus_1.1-11   
#  [31] future_1.33.2          shiny_1.8.0            digest_0.6.34         
#  [34] colorspace_2.1-0       patchwork_1.2.0        rprojroot_2.0.4       
#  [37] tensor_1.5             RSpectra_0.16-1        irlba_2.3.5.1         
#  [40] labeling_0.4.3         progressr_0.14.0       fansi_1.0.6           
#  [43] spatstat.sparse_3.0-3  httr_1.4.7             polyclip_1.10-6       
#  [46] abind_1.4-5            compiler_4.3.3         bit64_4.0.5           
#  [49] withr_3.0.0            fastDummies_1.7.3      highr_0.10            
#  [52] float_0.3-2            MASS_7.3-60.0.1        lmtest_0.9-40         
#  [55] httpuv_1.6.14          future.apply_1.11.2    goftest_1.2-3         
#  [58] glue_1.7.0             nlme_3.1-164           promises_1.2.1        
#  [61] grid_4.3.3             cluster_2.1.6          reshape2_1.4.4        
#  [64] generics_0.1.3         gtable_0.3.4           spatstat.data_3.0-4   
#  [67] tidyr_1.3.1            data.table_1.15.2      utf8_1.2.4            
#  [70] spatstat.geom_3.2-9    RcppAnnoy_0.0.22       ggrepel_0.9.5         
#  [73] RANN_2.6.1             pillar_1.9.0           stringr_1.5.1         
#  [76] spam_2.10-0            RcppHNSW_0.6.0         later_1.3.2           
#  [79] splines_4.3.3          dplyr_1.1.4            lattice_0.22-5        
#  [82] survival_3.5-8         bit_4.0.5              deldir_2.0-4          
#  [85] tidyselect_1.2.1       miniUI_0.1.1.1         pbapply_1.7-2         
#  [88] knitr_1.45             git2r_0.33.0           gridExtra_2.3         
#  [91] scattermore_1.2        RhpcBLASctl_0.23-42    xfun_0.42             
#  [94] matrixStats_1.2.0      stringi_1.8.3          workflowr_1.7.1       
#  [97] lazyeval_0.2.2         yaml_2.3.8             evaluate_0.23         
# [100] codetools_0.2-19       tibble_3.2.1           cli_3.6.2             
# [103] uwot_0.2.2.9000        xtable_1.8-4           reticulate_1.36.1     
# [106] munsell_0.5.0          jquerylib_0.1.4        Rcpp_1.0.12           
# [109] globals_0.16.3         spatstat.random_3.2-3  png_0.1-8             
# [112] parallel_4.3.3         ellipsis_0.3.2         dotCall64_1.1-1       
# [115] listenv_0.9.1          viridisLite_0.4.2      MatrixExtra_0.1.15    
# [118] scales_1.3.0           ggridges_0.5.6         leiden_0.4.3.1        
# [121] purrr_1.0.2            rlang_1.1.3