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

单细胞系列教程

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

混杂因素简介

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的表达也是重要的解释变量，另外一个显著的特征是batch比individual更多解释基因表达的差异。

其他影响因素

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

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

易生信系列培训课程，扫码获取免费资料

更多阅读

心得体会 TCGA数据库 Linux Python 

高通量分析 免费在线画图 测序历史 超级增强子

生信学习视频 PPT EXCEL 文章写作 ggplot2

后台回复“生信宝典福利第一波”获取教程合集

听说分享到朋友圈的朋友会在公众号周年庆时中奖 (大家还记得去年的大放送吧，不记得查查历史)