OSCA单细胞数据分析笔记7—Feature selection

对应原版教程第8章http://bioconductor.org/books/release/OSCA/overview.htm 细胞的测序结果一般都有1万多的基因。从表达水平来看，其中大部分的信息可能都是较为“普通”的。如何从中抽取出能够最大程度表征细胞间异质性的基因是本小节的主要目的。

笔记要点

• 1、背景知识
  • 1.1 为什么要挑选特定的基因
  • 1.2 怎么挑HVG
• 2、基因表达的方差分解
  • 2.1 基于log-counts的mean-variance拟合
  • 2.2 基于normalization counts的CV2变异系数
• 3、如何根据计算指标筛选高变基因
  • 3.1 基于Top X排名
  • 3.2 其它方法
• 4、确定高变基因之后sce对象的取舍
  • 4.1 仅保留高变基因信息（不建议）
  • 4.2 标记高变基因，降维设置subset.row=参数（建议）
• 5、补充：关于“技术误差”的进一步分解

1、背景知识

1.1 为什么要挑选特定的基因

• 单细胞数据分析的主要在于考虑细胞（群）之间的异质性，最直接的方式就是通过相同基因在不同细胞间的差异表达来描述。
• 不同细胞间的基因表达波动差异可能会由技术误差和生物水平的异质性两方面导致，而我们的目的就是找出主要由后者原因导致的高变基因(HVG, highly variable genes)，减少其它不太相关基因的信息干扰，用于下一步的降维。

1.2 怎么挑HVG

• 根据上面描述，最简单的想法就是计算基因在不同细胞间表达值的方差，取方差最大的基因即可；
• 但基因表达的总方差里有多少是技术误差，有多少是生物水平差异是需要分解的；
• 而且需要考虑到mean-variance的影响，即基因表达均值越大，可能受到越大的技术误差。

2、基因表达的方差分解

• 通过统计方法，计算基因的方差在多大程度上是由于不同类型细胞差异表达造成的；
• 下面介绍的两种方法有个统一的假设就是：大部分的基因表达波动是由于技术误差造成的。
• 示例数据集如下，已经质控，标准化处理。具体步骤见原教程

sce.pbmc
#class: SingleCellExperiment 
#dim: 33694 3985 
#assays(2): counts logcounts

sce.416b
#class: SingleCellExperiment 
#dim: 46604 185 
#assays(2): counts logcounts

2.1 基于log-counts的mean-variance拟合

绘制基因平均表达量（log转换）和方差的点图，拟合一条曲线，近似表示该表达量的基因受技术误差导致的方差。

library(scran)
dec.pbmc <- modelGeneVar(sce.pbmc)
fit.pbmc <- metadata(dec.pbmc)
plot(fit.pbmc$mean, fit.pbmc$var, xlab="Mean of log-expression",
     ylab="Variance of log-expression")
curve(fit.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)

如下图所示，高于拟合曲线的点到拟合曲线的竖直距离表示有意义的生物水平造成的基因表达方差。

具体计算结果如下--bio = total - tech。（至于负数的结果本身是没有意义的，可忽略）

dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),] 

引申两个参数：

1. 拟合中的可能存在的过拟合问题：这里主要针对基因表达量高，且方差大的少数点（基因）。造成曲线的膨胀,很大概率上高估了技术误差，拟合曲线表现为高跷的尾巴。可通过modelGeneVar()函数的density.weights参数设置(默认为T)。如下图，一般红线的拟合情况是我们期望看到的(density.weights=FALSE)

1. 针对多批次的高变基因的挑选，最直接的方法就是对每个批次单独拟合，再挑选。具体实现是在拟合函数中，设置block参数指定批次的设计 modelGenemodelGeneVar(sce.416b, "ERCC", block=sce.416b$block)

2.2 基于normalization counts的CV2变异系数

• 变异系数(Coefficient of variation)也是均值水平归一化的方差比较方法；
• 与2.1不同的是，这里使用的表达量是未经log转换的normalized counts；
• 但同样是拟合不同表达量水平下，技术误差造成的变异系数大小，然后取高于拟合曲线的点，认为是能够反应生物水平差异的基因。

dec.cv2.pbmc <- modelGeneCV2(sce.pbmc)
fit.cv2.pbmc <- metadata(dec.cv2.pbmc)
plot(fit.cv2.pbmc$mean, fit.cv2.pbmc$cv2, log="xy")
curve(fit.cv2.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)

• 具体计算结果如下--ratio = total / trend。注意ratio指标是用除法，而不是2.1里描述的减法。因为CV2=(标准差/均值)2，需要排除分母均值的干扰，才能使不同基因的CV2指标具有可比性。

dec.cv2.pbmc[order(dec.cv2.pbmc$ratio, decreasing=TRUE),]

• 如下图结果，较2.1的方法，由于未进行log压缩，更有可能挑选出均值水平较低的高变基因。

关于选择2.1还是2.2计算基因的筛选指标，如教程所说，二者都能很好的捕获基因的生物水平造成的差异性。但后面计算细胞间的距离是基于log-counts的方法，因此使用方法2.1可能会更具有统一性。

3、如何根据计算指标筛选高变基因

3.1 基于Top X排名

• 最简单的方法就是选择上一步计算的高变基因指标最显著的一批基因；
• 这会引出第二个问题：选多少个高变基因合适？一般的选择范围可能从500到5000不等；
• 一方面，选出很多高变基因可能会包含部分不感兴趣的无关基因的干扰；另一方面，选择的高变基因过少，又可能漏掉一些重要基因的信息。
• 因此，最好先对测序的细胞群体的异质性有一个大致评估。异质性高的话，高变基因就定义多一点；相反，则定义少一点高变基因。在此基础上，选一个大致的数量，跑接下来的流程；如果发现不合适再回到这里修改。而不是浪费太多时间，纠结最佳取值这个难题。

#基于log-count的计算指标
dec.pbmc <- modelGeneVar(sce.pbmc)
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n=1000)
#getTopHVGs(dec.pbmc, prop=0.1)
str(hvg.pbmc.var)
#chr [1:1000] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" "IGKC" "CCL5" "HLA-DRB1" ...

#基于log-count的计算指标
dec.pbmc <- modelGeneVar(sce.pbmc, )
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n=1000)
str(hvg.pbmc.var)
#chr [1:1000] "PPBP" "PRTFDC1" "HIST1H2AC" "FAM81B" "PF4" "GNG11" "KIAA1211" "HIST2H2BE" ...

3.2 其它方法

根据教程，再简单列举几种其它方法

3.2.1 基于FDR显著性

• 根据显著性的P值筛选高变基因，具体是设置getTopHVGs()函数的fdr.threshold参数。

getTopHVGs(dec.pbmc, fdr.threshold=0.05)

• 对于log-count指标，零假设就是bio<=0
• 对于CV2指标，零假设就是ratio<=1
• 根据设定的阈值，可以控制假阳性高变基因的数目。

3.2.2 取全部合格的高变基因

getTopHVGs(dec.pbmc, var.threshold=0)

• 对于log-count指标，取所有bio > 0的基因；
• 对于CV2指标，取所有ratio > 1的基因。
• 不做推荐，因为可能含有非常多的噪音基因

之后教程还举出一种方法是预先定义一组基因集作为作为高变基因。虽然与高变基因最初定义不符，但在某些情况下还是挺有意思的，就不展开了(其实没怎么看明白)，可参看原教程。

4、确定高变基因之后sce对象的取舍

确定高变基因后，后续降维操作主要根据能够表示细胞间生物水平差异的高变基因展开

4.1 仅保留高变基因信息（不建议）

dec.pbmc <- modelGeneVar(sce.pbmc)
chosen <- getTopHVGs(dec.pbmc, prop=0.1)
str(chosen)
sce.pbmc.hvg <- sce.pbmc[chosen,]
dim(sce.pbmc.hvg)
#[1] 1274 3985

4.2 标记高变基因，降维设置subset.row=参数（建议）

# Performing PCA only on the chosen HVGs.
library(scater)
sce.pbmc <- runPCA(sce.pbmc, subset_row=chosen)
reducedDimNames(sce.pbmc)

• 小技巧：可以通过rowSubset()函数，将高变基因信息储存到sce对象里(可以添加多次)

rowSubset(sce.pbmc,"hvg")=chosen
colnames(rowData(sce.pbmc))
rowData(sce.pbmc)[,2][rowData(sce.pbmc)[,3]]

5、补充：关于“技术误差”的进一步分解

• 如笔记一开始所说可将基因表达值的总方差分解为技术误差和生物水平的误差。但是实际上，前者还包括由于转录扰动造成的总体细胞表达水平的变化；可称之为“uninteresting” biological variation；例如管家基因house keeping gene，一般可忽略，统称为技术误差。
• 但在某些情况下，例如某些特定类型细胞的特异基因受到调控作用显著上调，即出现了一批（均值大、方差也大）的高变基因(hvg)，会造成拟合曲线的膨胀。即提高了该均值水平下的拟合技术误差的标准，从而可能漏掉潜在的高变基因。
• 针对此种情况，可使用外参转录本拟合单纯的技术误差。前提假设就是外参转录本的含量不受细胞内的生物过程调控。
• 如果数据集中没有外参转录本信息，那么可使用泊松分布近似拟合技术误差曲线。
• 相关函数如下，具体使用可参考原教程。

modelGeneVarWithSpikes(sce.416b, "ERCC")
modelGeneVarByPoisson(sce.pbmc)