表达矩阵处理—数据可视化 / 开普饭

单细胞RNA-seq的一个重要方面是控制批次效应。批量效应是在处理过程中添加到样品中的技术假象。例如，如果在不同实验室中或甚至在同一实验室中的不同日期制备两组样品，那么我们可以观察到一起处理的样品之间更大的相似性。在最坏的情况下，批量效应可能被误认为是真正的生物变异。理想情况下，我们期望看到来自同一个体的批次组合在一起，并且对应于每个个体的不同组。

library(SingleCellExperiment)
library(scater)
options(stringsAsFactors = FALSE)
umi <- readRDS("tung/umi.rds")
umi.qc <- umi[rowData(umi)$use, colData(umi)$use]
endog_genes <- !rowData(umi.qc)$is_feature_control

7.3.2 · PCA图

概述数据的最简单方法是使用主成分分析对其进行转换，然后可视化前两个主要成分。

主成分分析（PCA）（https://en.wikipedia.org/wiki/Principal_component_analysis）是一种统计程序，它使用转换将一组观察值转换为一组称为主成分（PC）的线性不相关变量值。主成分的数量小于或等于原始变量的数量。

在数学上，PC对应于协方差矩阵的特征向量。特征向量按特征值排序，因此第一主成分尽可能地考虑数据的可变性，并且每个后续成分在与前面的成分正交的约束下具有最高的方差（图中）以下是从这里（http://www.nlpca.org/pca_principal_component_analysis.html）获取的）。

7.3.2.1 · QC之前

没有log转换：

plotPCA(
umi[endog_genes, ],
exprs_values = "counts",
colour_by = "batch",
size_by = "total_features",
shape_by = "individual"
)

使用log转换：

plotPCA(
umi[endog_genes, ],
exprs_values = "logcounts_raw",
colour_by = "batch",
size_by = "total_features",
shape_by = "individual"
)

显然，对数转换对我们的数据是有益的-它减少了第一主成分的方差，并且已经分离了一些生物效应。而且，它使表达值的分布更正常。在下面的分析和章节中，我们将默认使用对数转换的原始计数。

但是，请注意，只有log转换不足以解释细胞之间的不同技术因素（例如测序深度）。因此，请不要使用logcounts_raw进行您的下游分析，而是作为最低合适的数据使用logcounts的SingleCellExperiment对象，这不只log数转换，也可由文库大小（例如CPM正常化）归一化。在课程中，我们logcounts_raw仅用于演示目的！

7.3.2.2 · 质量控制后

plotPCA(
umi.qc[endog_genes, ],
exprs_values = "logcounts_raw",
colour_by = "batch",
size_by = "total_features",
shape_by = "individual"
)

比较图7.17和图7.18，很明显，在质量控制后，NA19098.r2细胞不再形成一组异常值。

默认情况下，scater仅使用前500个最易变基因来计算PCA。这可以通过更改ntop参数来调整。

练习1如果使用所有14,214个基因，PCA图如何变化？或者只使用前50个基因？为什么第一个PC变化所引起的方差分数如此显着？

提示使用ntop函数的参数plotPCA。

我们的答案

如果您的答案不同，请将您的代码与我们的代码进行比较（您需要在打开的文件中搜索此练习）。

7.3.3 · tSNE图

用于可视化scRNASeq数据的PCA的替代方案是tSNE图。tSNE（https://lvdmaaten.github.io/tsne/）（t分布式随机邻域嵌入）将维数降低（例如PCA）与最近邻网络上的随机游走相结合，将高维数据（即我们的14,214维表达矩阵）映射到二维空间，同时保留细胞之间的局部距离。与PCA相比，tSNE是一种随机算法，这意味着在同一数据集上多次运行该方法将导致不同的图。由于算法的非线性和随机性，tSNE更难以直观地解释。为了确保可重复性，我们在下面的代码中修改随机数生成器的“种子”，以便我们始终获得相同的图。

7.3.3.1 · QC之前

plotTSNE(
umi[endog_genes, ],
exprs_values = "logcounts_raw",
perplexity = 130,
colour_by = "batch",
size_by = "total_features",
shape_by = "individual",
rand_seed = 123456
)

7.3.3.2 · 质量控制后

plotTSNE(
umi.qc[endog_genes, ],
exprs_values = "logcounts_raw",
perplexity = 130,
colour_by = "batch",
size_by = "total_features",
shape_by = "individual",
rand_seed = 123456
)

解释PCA和tSNE图通常具有挑战性，并且由于它们具有随机和非线性特性，因此它们不太直观。但是，在这种情况下，很明显它们提供了类似的数据图像。比较图7.21和7.22，再次清楚的是，在QC之后，来自NA19098.r2的样本不再是异常值。

此外，tSNE要求您提供的值perplexity反映用于构建最近邻网络的邻居数量; 高值会创建一个密集的网络，将细胞聚集在一起，而低值会使网络更稀疏，从而允许细胞群彼此分离。scater默认使用所有细胞的perplexity除以5（向下舍入）。

您可以在这里（http://distill.pub/2016/misread-tsne/）阅读更多关于使用tSNE的缺陷。

练习2当使用10或200的perplexity 时，tSNE图如何变化？perplexity 的选择如何影响结果的解释？

我们的答案

7.3.4 · 大练习

使用Blischak数据的reads执行相同的分析。使用tung/reads.rdsfile加载read SCE对象。完成后，请将您的结果与我们的结果进行比较（下一章）。

7.3.5 · sessionInfo( )

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
##
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-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=C
## [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] parallel stats4 methods stats graphics grDevices utils
## [8] datasets base
##
## other attached packages:
## [1] scater_1.6.3 ggplot2_2.2.1
## [3] SingleCellExperiment_1.0.0 SummarizedExperiment_1.8.1
## [5] DelayedArray_0.4.1 matrixStats_0.53.1
## [7] Biobase_2.38.0 GenomicRanges_1.30.3
## [9] GenomeInfoDb_1.14.0 IRanges_2.12.0
## [11] S4Vectors_0.16.0 BiocGenerics_0.24.0
## [13] knitr_1.20
##
## loaded via a namespace (and not attached):
## [1] viridis_0.5.0 httr_1.3.1 edgeR_3.20.9
## [4] bit64_0.9-7 viridisLite_0.3.0 shiny_1.0.5
## [7] assertthat_0.2.0 highr_0.6 blob_1.1.0
## [10] vipor_0.4.5 GenomeInfoDbData_1.0.0 yaml_2.1.17
## [13] progress_1.1.2 pillar_1.2.1 RSQLite_2.0
## [16] backports_1.1.2 lattice_0.20-34 glue_1.2.0
## [19] limma_3.34.9 digest_0.6.15 XVector_0.18.0
## [22] colorspace_1.3-2 cowplot_0.9.2 htmltools_0.3.6
## [25] httpuv_1.3.6.1 Matrix_1.2-7.1 plyr_1.8.4
## [28] XML_3.98-1.10 pkgconfig_2.0.1 biomaRt_2.34.2
## [31] bookdown_0.7 zlibbioc_1.24.0 xtable_1.8-2
## [34] scales_0.5.0 Rtsne_0.13 tibble_1.4.2
## [37] lazyeval_0.2.1 magrittr_1.5 mime_0.5
## [40] memoise_1.1.0 evaluate_0.10.1 beeswarm_0.2.3
## [43] shinydashboard_0.6.1 tools_3.4.3 data.table_1.10.4-3
## [46] prettyunits_1.0.2 stringr_1.3.0 munsell_0.4.3
## [49] locfit_1.5-9.1 AnnotationDbi_1.40.0 bindrcpp_0.2
## [52] compiler_3.4.3 rlang_0.2.0 rhdf5_2.22.0
## [55] grid_3.4.3 RCurl_1.95-4.10 tximport_1.6.0
## [58] rjson_0.2.15 labeling_0.3 bitops_1.0-6
## [61] rmarkdown_1.8 gtable_0.2.0 DBI_0.7
## [64] reshape2_1.4.3 R6_2.2.2 gridExtra_2.3
## [67] dplyr_0.7.4 bit_1.1-12 bindr_0.1
## [70] rprojroot_1.3-2 ggbeeswarm_0.6.0 stringi_1.1.6
## [73] Rcpp_0.12.15 xfun_0.1

7.4

数据可视化（reads）

library(scater)
options(stringsAsFactors = FALSE)
reads <- readRDS("tung/reads.rds")
reads.qc <- reads[rowData(reads)$use, colData(reads)$use]
endog_genes <- !rowData(reads.qc)$is_feature_control
plotPCA(
reads[endog_genes, ],
exprs_values = "counts",
colour_by = "batch",
size_by = "total_features",
shape_by = "individual"
)

plotPCA(
reads[endog_genes, ],
exprs_values = "logcounts_raw",
colour_by = "batch",
size_by = "total_features",
shape_by = "individual"
)
plotPCA(
reads.qc[endog_genes, ],
exprs_values = "logcounts_raw",
colour_by = "batch",
size_by = "total_features",
shape_by = "individual"
)

plotTSNE(
reads[endog_genes, ],
exprs_values = "logcounts_raw",
perplexity = 130,
colour_by = "batch",
size_by = "total_features",
shape_by = "individual",
rand_seed = 123456
)

plotTSNE(
reads.qc[endog_genes, ],
exprs_values = "logcounts_raw",
perplexity = 130,
colour_by = "batch",
size_by = "total_features",
shape_by = "individual",
rand_seed = 123456
)

sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
##
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-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=C
## [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] stats4 parallel methods stats graphics grDevices utils
## [8] datasets base
##
## other attached packages:
## [1] knitr_1.20 scater_1.6.3
## [3] SingleCellExperiment_1.0.0 SummarizedExperiment_1.8.1
## [5] DelayedArray_0.4.1 matrixStats_0.53.1
## [7] GenomicRanges_1.30.3 GenomeInfoDb_1.14.0
## [9] IRanges_2.12.0 S4Vectors_0.16.0
## [11] ggplot2_2.2.1 Biobase_2.38.0
## [13] BiocGenerics_0.24.0
##
## loaded via a namespace (and not attached):
## [1] viridis_0.5.0 httr_1.3.1 edgeR_3.20.9
## [4] bit64_0.9-7 viridisLite_0.3.0 shiny_1.0.5
## [7] assertthat_0.2.0 highr_0.6 blob_1.1.0
## [10] vipor_0.4.5 GenomeInfoDbData_1.0.0 yaml_2.1.17
## [13] progress_1.1.2 pillar_1.2.1 RSQLite_2.0
## [16] backports_1.1.2 lattice_0.20-34 glue_1.2.0
## [19] limma_3.34.9 digest_0.6.15 XVector_0.18.0
## [22] colorspace_1.3-2 cowplot_0.9.2 htmltools_0.3.6
## [25] httpuv_1.3.6.1 Matrix_1.2-7.1 plyr_1.8.4
## [28] XML_3.98-1.10 pkgconfig_2.0.1 biomaRt_2.34.2
## [31] bookdown_0.7 zlibbioc_1.24.0 xtable_1.8-2
## [34] scales_0.5.0 Rtsne_0.13 tibble_1.4.2
## [37] lazyeval_0.2.1 magrittr_1.5 mime_0.5
## [40] memoise_1.1.0 evaluate_0.10.1 beeswarm_0.2.3
## [43] shinydashboard_0.6.1 tools_3.4.3 data.table_1.10.4-3
## [46] prettyunits_1.0.2 stringr_1.3.0 munsell_0.4.3
## [49] locfit_1.5-9.1 AnnotationDbi_1.40.0 bindrcpp_0.2
## [52] compiler_3.4.3 rlang_0.2.0 rhdf5_2.22.0
## [55] grid_3.4.3 RCurl_1.95-4.10 tximport_1.6.0
## [58] rjson_0.2.15 labeling_0.3 bitops_1.0-6
## [61] rmarkdown_1.8 gtable_0.2.0 DBI_0.7
## [64] reshape2_1.4.3 R6_2.2.2 gridExtra_2.3
## [67] dplyr_0.7.4 bit_1.1-12 bindr_0.1
## [70] rprojroot_1.3-2 ggbeeswarm_0.6.0 stringi_1.1.6
## [73] Rcpp_0.12.15 xfun_0.1