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

par(mfrow = c(1,4))
boxplot(nUMI, main='nUMI'); 
boxplot(nGene, main='nGene'); 
boxplot(percent.mito, main='percent'); 
boxplot(det.rate, main='det')

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,
     xlab = "mean expression",
     ylab = "normalized dispersion")

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

UMI Normalization

Use quantile-normalization to make the distribution of each cell the same.

X = quantile_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