专栏首页单细胞天地表达矩阵处理—数据可视化

表达矩阵处理—数据可视化

7.清理表达矩阵

7.3数据可视化

7.3.1 · 简介

在本章中,我们将继续使用Tung前一章中生成的过滤数据集。我们将探索可视化数据的不同方法,以便您在质量控制步骤之后评估表达式矩阵发生的情况。scaterpackage提供了几个非常有用的功能来简化可视化。

单细胞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进行您的下游分析,而是作为最低合适的数据使用logcountsSingleCellExperiment对象,这不只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

本文分享自微信公众号 - 单细胞天地(sc-ngs),作者:温柔的坚果绷带

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-07-12

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • stLearn :空间轨迹推断

    空间信息在空间转录组中的运用 Giotto|| 空间表达数据分析工具箱 SPOTlight || 用NMF解卷积空间表达数据 Seurat新版教程:分析空间转录...

    生信技能树jimmy
  • 表达矩阵处理—表达QC(reads)

    通过比较图7.6和图7.13,很明显基于read的过滤比基于UMI的分析去除了更多的细胞。如果您返回并比较结果,您应该能够得出结论,ERCC和MT过滤器对于基于...

    生信技能树jimmy
  • 去除细胞效应和基因效应

    其实质量控制三部曲,还有一个很关键的点没有讲解,就是多个样本整合,并且区分批次效应和生物学差异。但是这个点很大程度是依赖于经验,就是说,要想搞清楚,需要写很多自...

    生信技能树jimmy
  • 最小的K个数:用快排的思想去解相关问题

    实现快速排序算法的关键在于先在数组中选择一个数字,接下来把数组中的数字分为两部分,比选择的数字小的数字移到数组的左边,比选择的数字大的数字移到数组的右边。 这个...

    猿人谷
  • 工具的使用-msfvenom

    Kali里面有一个集成好的工具是msfvenom,主要用于生成后门和软件捆绑后门,其免杀效果还算不错。

    字节脉搏实验室
  • <data>标签

    Html5知典
  • 组合类型与类型保护_TypeScript笔记9

    Object.assign能把source: U身上的可枚举属性浅拷贝到target: T上,因此返回值类型为T & U

    ayqy贾杰
  • centos7 install python3.7 with problem and how to fix it.

    <!-- p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 11.0px Menlo; color: #000000; ...

    超蛋lhy
  • 【翻译】Kotlin 1.1 新版本同样适合安卓开发者

    2017-04-29 by Liuqingwen | Tags: Kotlin 翻译 | Hits

    IT自学不成才
  • 加州无人车路测再添新玩家,全来自中国:Pony.ai和图森

    若朴 假装发自 湾区 量子位 报道 | 公众号 QbitAI ? 美国加州交通管理局(DMV)最新更新的文件显示,又有两家公司获准在加州展开无人车路测,而且这两...

    量子位

扫码关注云+社区

领取腾讯云代金券