专栏首页生信宝典如何火眼金睛鉴定那些单细胞转录组中的混杂因素

如何火眼金睛鉴定那些单细胞转录组中的混杂因素

识别基因表达检测影响因素

混杂因素简介

scRNA-seq数据会受到一些人为因素、操作偏差、批次等因素的影响。scRNA-seq分析的一个挑战是没有办法通过评估技术重复来区分生物和技术各自带来的变化有多大比例。前面的分析,我们考虑了批次效应,下面我们看下还有没有其它实验因素会影响单细胞基因表达检测并移除这些因素。scater包提供了一些评估实验因素和生物因素对表达数据影响的检测方法。我们用Blischak数据做例子展示其应用。

library(scater, quietly = TRUE)
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

umi_endog_genes <- !rowData(umi)$is_feature_control
umi_endog <- umi[umi_endog_genes,]
umi.qc <- umi[rowData(umi)$use, colData(umi)$use]
umi_qc_endog_genes <- !rowData(umi.qc)$is_feature_control
umi.qc_endog <- umi.qc[umi_qc_endog_genes,]

umi.qc数据集包含质控过滤后的细胞和基因。下一步是探索技术因素导致的表达变化以应用于下游的基因表达标准化分析中。

与主成分的相关性

质控后数据集的PCA展示

# plotPCA(
#     umi.qc[endog_genes, ],
#     exprs_values = "logcounts_raw",
#     colour_by = "batch",
#     size_by = "total_features"
# )
umi.qc_endog <- runPCA(umi.qc_endog, ncomponents=100, exprs_values = "logcounts_raw")
scater::plotPCA(
    umi.qc_endog,
    by_exprs_values = "logcounts_raw",
    colour_by = "batch",
    size_by = "total_features_by_counts",
    shape_by = "individual"
)

scater通过构建线性模型判断主成分与各个影响变量的相关性,从而判断哪些实验或质控变量导致细胞在主成分上的分布。

检测到的基因数与主成分的关系

# plotQC(
#     umi.qc[umi_qc_endog_genes, ],
#     type = "find-pcs",
#     exprs_values = "logcounts_raw",
#     variable = "total_features"
# )
umi.qc_endog <- runPCA(umi.qc_endog, ncomponents=500, exprs_values = "logcounts_raw")
explanatoryPCs <- getExplanatoryPCs(umi.qc_endog, variables = "total_features_by_counts")
#explanatoryPCs <- getExplanatoryPCs(umi.qc_endog)
plotExplanatoryPCs(explanatoryPCs, nvars_to_plot = 5, npcs_to_plot = 10)

确实,PC1几乎完全可以被检测到的基因数解释。从上面PCA图的结果也可以看出,延PC1从左至右,细胞检测到的基因数整体逐步降低的趋势。这也是scRNA-seq一个已经知道的现象,具体见http://biorxiv.org/content/early/2015/12/27/025528.

Explanatory variables (解释变量)

scater也可以把质控变量与所有基因分别进行线性模型拟合获取其边际 (marginal) ,绘制其概率密度分布图谱。

# umi.qc_endog <- normalize(umi.qc_endog)
ExplanatoryVariable <- getVarianceExplained(umi.qc_endog, exprs_values = "logcounts_raw",
                                variables=c("total_features_by_counts",
                                            "total_counts",
                                            "batch",
                                            "individual",
                                            "pct_counts_MT",
                                            "pct_counts_ERCC"))
plotExplanatoryVariables(ExplanatoryVariable)

# plotQC(
#     umi.qc[endog_genes, ],
#     type = "expl",
#     exprs_values = "logcounts_raw",
#     variables = c(
#         "total_features",
#         "total_counts",
#         "batch",
#         "individual",
#         "pct_counts_ERCC",
#         "pct_counts_MT"
#     )
# )

结果显示检测到的基因数(total_features_by_counts)和测序深度 (total_counts)对基因表达的贡献度很大。因此在基因表达标准化过程中需要考虑移除这些因素的影响或整合到下游的统计分析模型中。ERCC的表达也是重要的解释变量,另外一个显著的特征是batchindividual更多解释基因表达的差异。

其他影响因素

除了考虑校正批次影响 (依赖于实验记录的外部信息),还有其他技术因子需要考虑如何进行抵消。一个常用方法是scLVM (https://github.com/PMBio/scLVM) 是允许识别和移除细胞周期程序性死亡引入的影响。(Seurat+Scran也可以)

另外,不同的实验方案对转录本的覆盖偏好也不同,这一偏好依赖于A/T的平均含量或短的转录本的捕获能力。理想情况下,我们需要消除这些所有的差异和偏差。

本文分享自微信公众号 - 生信宝典(Bio_data),作者:宋小云

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

原始发表时间:2019-04-17

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 什么?你做的差异基因方法不合适?

    上一步,我们鉴定出了重要的干扰因素和解释变量可能对表达定量带来影响。scater允许在后续统计模型中引入这些变量来屏蔽技术操作带来的影响,或者可以给函数norm...

    生信宝典
  • Graphpad,经典绘图工具初学初探

    大多数科研文章都离不开图表,尤其是图,熟悉一些绘图软件,并将图在文章和PPT中展示出来,是科研训练的重要内容。漂亮的文章配图能给自己的工作加不少分,生信宝典推出...

    生信宝典
  • Weblogo |Seq logo 在线绘制工具

    seqlogo图可以直观清晰的反应序列偏好特征,每个位置出现的碱基或氨基酸类型反映了该位置序列的偏好性,每个字母的大小与该碱基在该位置上的出现频率成正相关。

    生信宝典
  • 【热点】大数据分析的八大趋势

    Intuit公司的数据工程副总裁Bill Loconzolo,双脚踏进了数据湖。.Smarter Remarketer的首席数据科学家Dean Abbott直接...

    小莹莹
  • qiime2+picrust1学习笔记

    一直迷惑于如何把qiime2和picrust结合起来用来分析16S的数据,直到这两天,看到了微生太公众号的视频教程,才有了眉目,原来如此。详细视频教程可以查找相...

    用户1075469
  • 个人博客网站接入来必力评论系统

    网易在7月6日正式发表了公告,通知用户即将停止服务(7月底),所以我们这些免费的使用者也不得不换窝了。 第三方评论 在之前,第三方评论系统主要有畅言、来比力、网...

    xiangzhihong
  • MSYS2:获取本机的ip地址

    MSYS2虽然是个linux shell环境,但如果要获取网卡的信息,还是需要windows平台提供的命令 参照这篇文章 《bat脚本 - 获取局域网内的本...

    用户1148648
  • BAT介入教育O2O的背后是人工智能战略

    大数据文摘
  • 10w定时任务,如何高效触发超时

    一、缘起 很多时候,业务有定时任务或者定时超时的需求,当任务量很大时,可能需要维护大量的timer,或者进行低效的扫描。 例如:58到家APP实时消息通道系统,...

    架构师之路
  • python yield函数深入浅出理解

    首先关于生成器的那些事: 1.通常的for…in…循环中,in后面是一个数组,这个数组就是一个可迭代对象,类似的还有链表,字符串,文件。它的缺陷是所有数据都...

    学到老

扫码关注云+社区

领取腾讯云代金券