10X Pipeline

Read Data : PBMC 27K

Read data. Keep the gene names separately.

orig = readMM('data/unnecessary_in_building/pbmc3k/matrix.mtx')
orig_genenames = read.table('data/unnecessary_in_building/pbmc3k/genes.tsv',
                            stringsAsFactors = FALSE)

Quality Control and Cell Selections

Get summary of the cells. The function “cellFilter” removes abnormal cells based on the read counts. The arguments minGene and maxGene restrict the number of genes detected in each cell. In 10X and Drop-seq data, having lower limit of 500 and upper limit of 2000 are generally appropriate. The cells with greater than 2000 detected genes are likeliy to be doublets, and those with less than 500 have too many dropouts. The default values are -Inf and Inf, so users are recommended to inspect the histogram of gene counts and determine the bounds. The function also requires the gene names to discover the mito-genes. Cells with high mitochondrial read proportion can indicate apoptosis. The default is 0.1.

nGene = colSums(orig > 0)
hist(nGene)

summaryX = cellFilter(X = orig, 
                      genenames = orig_genenames$V2,
                      minGene = 500, 
                      maxGene = 2000, 
                      maxMitoProp = 0.1)
tmpX = summaryX$X
nUMI = summaryX$nUMI
nGene = summaryX$nGene
percent.mito = summaryX$percent.mito
det.rate = summaryX$det.rate
dim(tmpX)

[1] 32738  2471

Gene Filtering using normalized dispersion

Next find variable genes using normalized dispersion. First, remove the genes where counts are 0 in all the cells, so that we use genes with at least one UMI count detected in at least one cell are used. Then genes are placed into a number of bins (user’s choice in “bins” parameter in “dispersion” function, default is 20) based on their mean expression, and normalized dispersion is calculated as the absolute difference between dispersion and median dispersion of the expression mean, normalized by the median absolute deviation within each bin. (Grace Zheng et al., 2017)

X = tmpX[rowSums(tmpX) > 0, ]
genenames = orig_genenames[rowSums(tmpX) > 0, ]

#gene filter by dispersion
disp = dispersion(X, bins = 20)
plot(disp$z ~ disp$genemeans)

select = which(abs(disp$z) > 1)
X = X[select, ]
genenames = genenames[select,]
dim(X)

[1]  833 2471

UMI Normalization

Normalize by the library size so that each cell has the same read count. First divide UMI counts by the total UMI counts in each cell, and then multiply with the median of the total UMI counts across cells. Note that entries with 0 are unaffected by the normalization and we can further clean the data based on the number of detected genes.

X = UMI_normalize(as.matrix(X))

Correct Detection Rate

#take log
logX = as.matrix(log(X + .1))

#check dependency
out = correct_detection_rate(logX, det.rate)

#regress out
log.cpm = out$residual

Dimension reduction on the data for visualization.

pc.base = irlba(log.cpm, 20)

plot(pc.base$d, ylab = "singular values")

tsne.base = Rtsne(pc.base$v[,1:10], dims=2, perplexity = 100, pca=FALSE)

rm(pc.base)

Run SLSL on the log.cpm matrix.

out.base = SLSL(log.cpm, log=FALSE,
                filter = FALSE,
                correct_detection_rate = FALSE,
                klist = c(200,250,300),
                sigmalist = c(1,1.5,2),
                kernel_type = "pearson",
                verbose=FALSE)
tab = table(out.base$result)
plot(tsne.base$Y, col=rainbow(length(tab))[out.base$result],
     xlab = 'tsne1', ylab='tsne2', main="base SLSL", cex = 0.5)

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.5

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/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] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2             SCNoisyClustering_0.1.0   
 [3] plotly_4.8.0               gplots_3.0.1              
 [5] diceR_0.5.1                Rtsne_0.13                
 [7] igraph_1.2.2               scatterplot3d_0.3-41      
 [9] pracma_2.1.4               fossil_0.3.7              
[11] shapefiles_0.7             foreign_0.8-71            
[13] maps_3.3.0                 sp_1.3-1                  
[15] caret_6.0-80               lattice_0.20-35           
[17] reshape_0.8.7              dplyr_0.7.6               
[19] quadprog_1.5-5             inline_0.3.15             
[21] matrixStats_0.54.0         irlba_2.3.2               
[23] Matrix_1.2-14              ggplot2_3.0.0             
[25] MultiAssayExperiment_1.6.0

loaded via a namespace (and not attached):
  [1] colorspace_1.3-2            class_7.3-14               
  [3] mclust_5.4.1                rprojroot_1.3-2            
  [5] XVector_0.20.0              RcppArmadillo_0.8.600.0.0  
  [7] GenomicRanges_1.32.6        pls_2.6-0                  
  [9] DRR_0.0.3                   prodlim_2018.04.18         
 [11] lubridate_1.7.4             codetools_0.2-15           
 [13] splines_3.5.1               R.methodsS3_1.7.1          
 [15] robustbase_0.93-2           knitr_1.20                 
 [17] RcppRoll_0.3.0              jsonlite_1.5               
 [19] workflowr_1.1.1             broom_0.5.0                
 [21] ddalpha_1.3.4               kernlab_0.9-26             
 [23] R.oo_1.22.0                 sfsmisc_1.1-2              
 [25] httr_1.3.1                  compiler_3.5.1             
 [27] backports_1.1.2             assertthat_0.2.0           
 [29] lazyeval_0.2.1              htmltools_0.3.6            
 [31] tools_3.5.1                 gtable_0.2.0               
 [33] glue_1.3.0                  GenomeInfoDbData_1.1.0     
 [35] reshape2_1.4.3              Rcpp_0.12.18               
 [37] Biobase_2.40.0              gdata_2.18.0               
 [39] nlme_3.1-137                iterators_1.0.10           
 [41] timeDate_3043.102           gower_0.1.2                
 [43] stringr_1.3.1               gtools_3.8.1               
 [45] DEoptimR_1.0-8              zlibbioc_1.26.0            
 [47] MASS_7.3-50                 scales_0.5.0               
 [49] ipred_0.9-6                 parallel_3.5.1             
 [51] SummarizedExperiment_1.10.1 yaml_2.2.0                 
 [53] rpart_4.1-13                stringi_1.2.4              
 [55] S4Vectors_0.18.3            foreach_1.4.4              
 [57] caTools_1.17.1.1            BiocGenerics_0.26.0        
 [59] BiocParallel_1.14.2         lava_1.6.2                 
 [61] geometry_0.3-6              GenomeInfoDb_1.16.0        
 [63] rlang_0.2.1                 pkgconfig_2.0.1            
 [65] bitops_1.0-6                evaluate_0.11              
 [67] purrr_0.2.5                 bindr_0.1.1                
 [69] htmlwidgets_1.2             recipes_0.1.3              
 [71] CVST_0.2-2                  tidyselect_0.2.4           
 [73] plyr_1.8.4                  magrittr_1.5               
 [75] R6_2.2.2                    IRanges_2.14.10            
 [77] dimRed_0.1.0                DelayedArray_0.6.2         
 [79] pillar_1.3.0                whisker_0.3-2              
 [81] withr_2.1.2                 survival_2.42-6            
 [83] abind_1.4-5                 RCurl_1.95-4.11            
 [85] nnet_7.3-12                 tibble_1.4.2               
 [87] crayon_1.3.4                KernSmooth_2.23-15         
 [89] rmarkdown_1.10              grid_3.5.1                 
 [91] data.table_1.11.4           git2r_0.23.0               
 [93] ModelMetrics_1.1.0          digest_0.6.15              
 [95] tidyr_0.8.1                 R.utils_2.6.0              
 [97] stats4_3.5.1                munsell_0.5.0              
 [99] viridisLite_0.3.0           magic_1.5-8