【生信文献200篇】02 用米氏方程解决单细胞转录组dropout现象

0.文章简介

本次要介绍的这篇文章提出了一个算法，被包装到了R包，是：M3Drop ， 文章最开始 2017年发表在bioRxiv，后来（2019）发表在Bioinformatics。

2017

英文标题：Modelling dropouts for feature selection in scRNASeq experiments

中文标题：为了挑选scRNASeq后的重要基因而对缺失建模

期刊：《bioRxiv : the preprint server for biology》

发表时间：2017.03

研究领域：R包

DOI号：10.1093/bioinformatics/bty1044

2019

英文标题：M3Drop: dropout-based feature selection for scRNASeq

中文标题：M3Drop是一种scRNASeq中基于缺失的重要基因挑选新算法

期刊：《Bioinformatics》

影响因子：5.61 发表时间：2019.08

DOI号：10.1093/bioinformatics/bty1044

1.背景

dropout是scRNASeq数据的一大技术难点，即很多基因在某些细胞根本就不表达，但是在另外的细胞却高表达。单细胞转录组数据里面的dropouts可以达到50%，但是通常认为这个dropouts是因为在文库构建的过程中，有部分基因没有被成功的反转录，是一个酶促反应，继而可以与米氏方程联系起来。

单细胞转录组测序的确可以一次性对所有细胞都检测到上千个基因的表达，但是，大多数情况下，只有其中的少部分基因是有生物学意义的，比如可以区分不同的细胞类型，或者细胞应对外界刺激。而且大多数基因之所以在不同的细胞里面表达有差异，其实是技术限制，背景噪音。这些技术限制，包括批次效应，都会阻碍我们发现那些真正的有生物学意义的基因。

寻找重要基因（feature selection）分析就是用来去除那些技术噪音相关基因，可以显著的提高信噪比，降低后续分析的复杂度。挑选到的跟feature相关的基因集，有点类似于在某些组间差异表达的基因集，都需要后续功能注释。

作者提出了两种针对scRNASeq数据的无监督feature selection的方法，与以往的所有方法不同的是，作者提出的两种方法都是基于基因丢失率而不是方差：

1.用M3Drop(Michaelis-Menten equation)来拟合全长转录本数据（例如Smartseq2）;

2.用depth adjusted negative binomial (DANB) 来拟合基于UMI的表达量数据（例如10X Chromium）。

经过与现有的几种方法的对比，证明了这两种新方法的优势：比别的方法好，可以克服数据集之间的批次效应，并减少数据集内的技术噪音，以揭示潜在的生物学过程。

米氏方程（Michaelis-Menten equation）是表示一个酶促反应的起始速度与底物浓度关系的速度方程:

v=Vmax × [S] /（Km+[S])

在酶促反应中，在低浓度底物情况下，反应相对于底物是一级反应（first order reaction）；而当底物浓度处于中间范围时，反应（相对于底物）是混合级反应（mixed order reaction）。当底物浓度增加时，反应由一级反应向零级反应（zero order reaction）过渡。

在假定存在一个稳态反应条件下推导出来的，其中 Km 值称为米氏常数，Vmax是酶被底物饱和时的反应速度，[S]为底物浓度。

Km值的物理意义为反应速度（v）达到1/2Vmax时的底物浓度（即Km=[S]），单位一般为mol/L，只由酶的性质决定，而与酶的浓度无关。可用Km的值鉴别不同的酶。

知道了这些背景知识，应该就会想到，这个由米氏方程包装的M3Drop包是怎么与dropout rate（多0表达的情况叫做dropout rate），gene expression结合起来？且到底有什么优势？

2.研究方法

1. 实验方法

1. Negative Binomial Models 负二项模型

这个是被应用的最广泛的转录组表达数据分布模型。但是对单细胞转录组测序数据来说，因为有很高的dropout情况，导致模型失准，所以就提出来了zero-inflated negative binomial models。

2. Zero-Inflation Models

是在原始的负二项分布数据里面随机挑选一些低表达量基因，给它们人为赋值为0表达量值。

3. Feature Selection Methods，比较了9种 feature selection 方法，每种方法都被用来从生物重要性最高到最不重要地对基因进行排序，并与各自的Bulk-RNAseq数据中的基本事实进行了比较，算法如下：

• M3Drop dropouts-mean expression curve (M3Drop)mean expression relationship (HVG)the dispersion-mean expression relationship fit by DANB (NBDisp)the dropouts-mean expression relationship fit by DANB (NBDrop)by the magnitude of their loadings in principal component analysis (PCA)by their relative Gini index (Gini)by the strength of their most negative gene-gene correlation (Cor)the squared coefficient of variation (CV2)

4. 这些算法都不需要预先对样本进行分类，是无监督的算法。因为这些差异表达方法的性能在很大程度上取决于预定义分组的质量/可靠性。

• differentially variable (DV) geneshighly variable (HV) genesdifferentially expressed (DE) genes

单细胞转录组数据的batch effects比较严重，所以 feature selection 过程的一个主要目的就是降低技术误差的影响，集中在有生物学意义的差异上面。

2. 验证方法的可行性和优势

1. Accuracy using bulk RNASeq ground truth

作者比较了 2个公共数据集，都是小鼠的胚胎细胞，含有17~255个细胞的测序数据

• Tung et al. (2017) [12] considered iPSCs from three different individuals and performed three replicates of UMI-tagged scRNASeq and three replicates of bulk RNASeq for each. (GSE77288 ).For Kolodziejczyk et al. (2015)，we considered ESCs grown under two conditions: alternative 2i and serum for which there were three replicates of scRNASeq and two replicates of bulk RNASeq.( E-MTAB-2600 )

对bulk转录组数据用了3种方法找差异基因，分别是 DESeq2,edgeR,limma-voom

只有3种方法都是 5% FDR的差异基因才认为是阳性标准基因集，那些3种方法都在 20% FDR的非差异基因认为是阴性金标准。

• 1,915 positives, and 8,398 negatives for the iPSCs709 positives and 11,278 negatives for the ESCs

有了这些基因，就可以计算ROC。

2. Single-cell RNASeq datasets

• 13个公共的scRNASeq数据集（文中的表1）

3. Simulated datasets

We simulated UMI-tagged data using the depth-adjusted negative binomial model fit to one of the three UMI-tagged datasets

• Tung （GSE77288 ）Zeisel（GSE60361 ）Klein （GSE65525）

We simulated full-transcript data using a zero-inflated negative binomial model fit to each of three full-transcript datasets

• Pollen（SRP041736）Buettner （E-MTAB-2805 ）Kolodziejczyk（E-MTAB-2600）

3.检验是否去除了批次效应

1. Reproducibility

• 我们使用五个完整转录的单细胞RNA序列数据集，检测从受精到囊胚的小鼠胚胎发育用log后的表达式值计算Anxa2与S100a10之间的Pearson相关性。

2. Identification of TE and ICM

• 通过将log后的表达值重新转化到Z-scores，将未鉴定的胚泡样品分配给ICM或TE。计算3个ICM标记(Sox2，Pou5f1/Oct4，Nanog)和3个TE标记（Elf5, Eomes, Cdx2）的平均得分，并将细胞归入得分较高的类型。

4.Code/Data Availability

可以在作者的github上免费获取M3Drop和DANB两个包，且获得使用代码。

3.研究结果

1

Michaelis-Menten Modelling of Dropouts (M3Drop) fits full length transcript data

使用Michaelis-Menten函数对 dropout rate 和 mean expression之间的关系进行建模。米氏方程：

其中S是基因在所有细胞中的平均表达，K_M是米氏常数，P_{dropout}表示cDNA的数量在任何细胞中达到某些实验特定的检测阈值的概率。

• 我们比较了M3Drop模型与两个另外的模型SCDE（SCDE是第一个特意针对单细胞转录组测序数据的差异分析而设计的，用贝叶斯统计方法把表达矩阵拟合到 zero-inflated negative binomial 分布模型里面。）和ZIFA在六个完整的scRNASeq数据集上的拟合结果。发现M3Drop模型的错误率最低，且拟合效果比其他两种模型合适。

• 同样，在2019年的文章中，模拟数据集中ZINB模型中两种基于Dropout的feature selection方法NBDrop和M3Drop均比基于方差的方法要好得多，consensus features也是如此。

2

Depth-adjusted Negative Binomial fits UMI-tagged data

带UMI标记的数据服从负二项分布。我们引入了深度调整负二项模型(DANB)，该模型将特定细胞的检测率合并到负二项模型中。

• 使用 dropout rate 拟合的DANB模型与适用于归一化数据的基本负二项模型进行了比较。明显的可以看出DANB模型比基本模型更适合带有UMI标记的数据集，错误率更低。

• While M3Drop fits full-transcript scRNASeq datasets well, data quantified using unique molecular identifiers (UMI) to obtain raw transcript counts often does not fit the Michaelis-Menten function.

3

Feature Selection

在单细胞RNAseq数据集中，经常使用重要基因的挑选以降维和减少技术干扰。与差异表达不同，重要基因的挑选对生物群体是不可知的，这使得它成为一个更困难的问题。目前已有的寻找单细胞转录组测序数据中的feature selection的方法都不够好，比如 scLVM 主要是根据先验基因集，比如cell-cycle or apoptosis来区分细胞。与此相反，基于 highly variable genes (HVG) 的方法挑选到的变化量大的那些基因很可能是技术带来的误差。而且低表达量基因的变动往往大于高表达量基因，而且所谓的表达变化大也并没有很好的生物学解释。一个比较好理解的概念是差异基因方法例如 SCDE ，但是需要预先把细胞群体分组后进行比较才能得到，而很多时候细胞太相似了，没办法很好的分开。像PCA或者t-SNE这样的降维方法也可以用来挑选重要基因，但它们也受制于系统误差或者批次误差等等。

对于UMI标记的数据，PCA和Cor的表现最好，以ROC曲线下的总面积(AuC=0.70)衡量，将假阳性率和真阳性率汇总为单一质量分数，并以排名前2,000的基因中真阳性的数量，即精确度来衡量(图2Ai，Bi)。

• 相比之下，这些方法在处理全长转录组数据时表现最差(图2Aii，Bii)。而M3Drop和两种DANB方法在全长转录组数据上的表现最好。有趣的是，当考虑到前2,000个Feature相关的基因时，基于dropout rate的新方法在两种数据类型上都优于基于方差的同类型方法，NBDrop>NBDisp和M3Drop>HVG

• 为了测试除了差异表达(DE)基因之外，用于识别差异变量(DV)或一般高变量(HV)基因的Feature selection方法的性能，我们使用深度调整的负二项模型或基于M3Drop的零膨胀负二项模型来模拟数据。我们将每个模型分别适合于三个不同的用户界面标记的或完整的数据集(图2C)。M3Drop在全长转录组数据上的表现比在UMI标记的数据上要好得多（Cii），而基尼指数只在UMI标记的数据上表现良好（Ci）。与DE基因相比，在两种数据类型中都没有通过Cor检测到差异变量基因。

4

Feature selection makes it possible to overcome batch effects

scRNASeq对批次效应非常敏感，这些批次效应包括实验重复之间或数据集之间的各种技术混杂因素，如测序效率、试剂质量等。Feature selection的目标之一是通过只关注与生物最相关的基因来减少这种技术噪音对下游分析的影响。如果Feature捕捉到了数据中与生物相关的方面，它们在检查同一生物系统的数据集中应该是一致的，而主要受技术噪声影响的Feature预计会有所不同。我们通过考虑前2,000个特征在检查植入前小鼠胚胎的五个数据集上的重现性来测试Feature selection是否有助于克服批次效应。数据集的大小从17个细胞到255个细胞不等，涵盖了从受精卵到囊胚的发育过程。

所有的Feature selection方法在不同的发展数据集之间显示出显著的重叠，看图中柱形图的长度，基于方差的方法比基于dropout rate的方法重复性较差(图3A)。Cor、M3Drop和NBDrop是最具重复性的方法。

当使用主成分分析组合和可视化这五个数据集时，所有Feature selection方法都大大改善了按发育时间点而不是按起源数据集对细胞的聚类(图3B，C)。

通过检查这五个数据集最具重复性的基因，我们发现了Anxa2和S100a10。从4细胞期(r=0.32，p=0.0007)到8和16细胞期(r=0.43，p=2×10^{-7})，这些基因的表达显著相关，直到囊胚分化到滋养外胚层的单细胞水平，在滋养外胚层的单细胞水平上大致呈线性关系，在细胞内团块细胞内又恢复不相关，r=0.49，p<10^{-10}(图3 D-F)。

M3Drop使用Feature selection后的囊胚细胞聚类清楚地表明存在两种细胞类型，我们使用标记基因识别为TE和ICM(图4A)

我们对此数据集使用M3Drop。相反，在不执行 feature selection方法的情况下对单元进行群集，建议使用2到5种单元类型的混合(图4B)。在此数据集中，只有M3Drop、DANB和HVG能够区分ICM和TE。

4.总结

1.我们的两种方法M3Drop和DANB是基于dropout rate，比别的基于variance的方法好。

2.相比之下，我们的方法既不依赖技术突增，也不改变观察到的表达值。

3.可以克服数据集之间的批次效应，并减少数据集内的技术噪音，以揭示潜在的生物学过程。revealing biological signals otherwise hidden beneath technical noise.

5.补充

单细胞转录组数据文章一般分成下面两大类，有兴趣的可以自己去了解一下：

第一大类是：deep sequencing of full-transcripts for a relatively small number of cells代表性的文章如下:

• Accounting for technical noise in single­cell RNA­seq experiments. Nat. Methods 10, 1093–1095 (2013).
• Fast, scalable and accurate differential expression analysis for single cells. (2016). doi:10.1101/049734
• Single­cell RNA­seq reveals dynamic, random monoallelic gene expression in mammalian cells. Science 343, 193–196 (2014). 14. Brennecke, P. et al. Accounting for technical noise in single­cell RNA­seq experiments. Nat. Methods 10, 1093–1095 (2013).
• Dynamics of Global Gene Expression Changes during Mouse Preimplantation Development. Dev. Cell 6, 117–131 (2004).
• Roles of CDX2 and EOMES in human induced trophoblast progenitor cells. Biochem. Biophys. Res. Commun. 431, 197–202 (2013).

第二类是：high-cell number, low-depth sequencing of 3’ or 5’ ends of transcripts tagged with unique molecular identifiers 代表性的文章是：

• Quantification noise in single cell experiments. Nucleic Acids Res. 39, e124 (2011).
• Quantification of mRNA in single cells and modelling of RT­qPCR induced noise. BMC Mol. Biol. 9, 63 (2008).
• ZIFA: Dimensionality reduction for zero­inflated single­cell gene expression analysis. Genome Biol. 16, 241 (2015).
• DNA methylation dynamics during epigenetic reprogramming in the germline and preimplantation embryos. Genes Dev. 28, 812–828 (2014).
• Genetic programs in human and mouse early embryos revealed by single­cell RNA sequencing. Nature 500, 593–597 (2013).