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

单细胞系列教程

收藏北大生信平台” 单细胞分析、染色质分析” 视频和PPT分享
Science: 小鼠肾脏单细胞转录组+突变分析揭示肾病潜在的细胞靶标
10X单细胞测序分析软件:Cell ranger，从拆库到定量
Hemberg-lab单细胞转录组数据分析（一）- 引言
Hemberg-lab单细胞转录组数据分析（二）- 实验平台
Hemberg-lab单细胞转录组数据分析（三）- 原始数据质控
Hemberg-lab单细胞转录组数据分析（四）- 文库拆分和细胞鉴定
Hemberg-lab单细胞转录组数据分析（五）- STAR, Kallisto定量
Hemberg-lab单细胞转录组数据分析（六）- 构建表达矩阵，UMI介绍
Hemberg-lab单细胞转录组数据分析（七）- 导入10X和SmartSeq2数据Tabula Muris
Hemberg-lab单细胞转录组数据分析（八）- Scater包输入导入和存储
Hemberg-lab单细胞转录组数据分析（九）- Scater包单细胞过滤
Hemberg-lab单细胞转录组数据分析（十）- Scater基因评估和过滤
Hemberg-lab单细胞转录组数据分析（十一）- Scater单细胞表达谱PCA可视化
Hemberg-lab单细胞转录组数据分析（十二）- Scater单细胞表达谱tSNE可视化
单细胞分群后，怎么找到Marker基因定义每一类群？
DESeq2差异基因分析和批次效应移除

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

混杂因素简介

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_controlumi_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的平均含量或短的转录本的捕获能力。理想情况下，我们需要消除这些所有的差异和偏差。

往期精品(点击图片直达文字对应教程)

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集