扫码联系客服
单细胞RNA测序数据基因表达不同转化方法比较
研究背景
单细胞RNA测序从细胞层面研究肿瘤异质性，获得肿瘤组织内部单个细胞的基因结构和基因表达状态，同时捕获细胞之间通讯和发育轨迹，相比RNAseq，对基因表达定量更加精准。目前主流技术是基于barcode的单细胞测序，将单细胞悬液样品和带有barcode的水凝胶珠子，通过微流体芯片，包裹在一个油滴之中。在油滴中进行逆转录后，每个单细胞的cDNA文库，就带上了独一无二的barcode。由于每个细胞中RNA含量极低，对每个转录本加入UMI，可有效降低PCR扩增错误导致的计数误差。基因x细胞的表达计数矩阵是单细胞RNA测序数据分析基础，无论是聚类分析还是特定细胞亚群表达基因分析，或者细胞通讯分析，都需要输入这个计数表格数据，且基于统计推断的单细胞分析工具都有个假设前提是基因在不同细胞之间表达计数具有相似方差。
研究结论
来自德国海德堡大学的科研人员在模拟和真实单细胞基准数据基础上，比较了4类对基因x细胞的表达计数矩阵进行转化以获得不同细胞之间动态相似方差的统计方法：delta、残差、直接表达推测和因子分析。结果发现：一种简单的计数增量然后log转化加PCA分析，跟有着较好理论基础的残差方法、表达推测和因子分析等有着同样甚至更加优异表现。
研究结果
基础理论
转录组测序数据中高表达基因往往具有较大的计数方差，不同实验批次、测序数据质量、分析流程差异都可能带来较大干扰，导致分析结果偏离真实情况。目前RNA表达数据的标准统计方法对方差一致的数据表现最好，对于方差不齐的RNA数据，一般采用抽样分布建模来预估分布参数或者数据转化以获得方差齐性进行预处理。建模参数推断往往耗时长且需要大量计算资源，更有效的方法是对计数进行转化获得方差齐性。对于单细胞RNA测序数据，UMI标记转录本计数比较符合负二项分布，各种转换方法均在此基础上进行。
转化方法
delta
通过简单的非线性转化，如g(y) = log (y + y0)，其中y代表基因在单个细胞中的原始计数，y0为固定值，使得方差在动态范围内更加接近，其中y0的选择对结果准确性很重要。单细胞RNA测序数据中每个细胞UMI总数受采样效率和细胞大小影响，在进行方差稳定转化前，需要对每个细胞确定一个放缩因子s，常规方法是让所有细胞的放缩因子之和接近于1，所以s=单个细胞UMI总数/L，其中L表示所有细胞的UMI总数之和的平均值，然后再进行方差稳定转化，即g(y) = log(y/s + y0)。Seurat中L设置为固定值10,000，相当于y0=0.5，α=0.5，更接近于真实数据中观察到的过度离散状态(overdispersion)，CPM计数本质是L=106，相当于y0=0.005，α=50，比典型单细胞表达数据的离散度大了两个数量级。对于模拟数据，该类别分别基于acosh转化、y0=1、y0=1/4(α)、CPM计数、log转化后选择表达高变化基因（HVG）或者zscore进行额外处理等方法进行评估。
残差
前面提到，最接近UMI标记单细胞表达计数的统计模型是负二项分布，因此可以采用广义线性模型（GLM）对实际表达数据进行回归分析，以期达到Pearson残差最小化，而负二项分布参数均值通过每个基因拟合对应的斜率和截距以及单个基因放缩因子的log变换得到。对于模拟数据，该类别分别基于剪切和未剪切的Pearson残差、随机分位数残差、选择表达高变化基因（HVG）或者zscore进行log转化后处理等方法进行评估。
基因表达推测
该方法的基本思想是以观察到的基因计数为前提，直接推断该基因表达状态，贝叶斯模型是这个问题的典型应用。代表方法有基于全贝叶斯模型的Sanity，它另外包括两种细分方法，一种通过计算基因表达log转化的均值和方差来衡量细胞之间的距离，找出距离最近的k个细胞，另一种Santiy MAP是返回后验最大值。Dino拟合混合负二项分布，返回随机样本后验值，Normalisr返回每个计数的最小均方误差估计。对于模拟数据，在上述提到的四种方法中进行了评估。
因子分析
该方法不对计数进行转化，而是将多维单细胞数据降到低维空间，在低维空间中，同样符合负二项分布的UMI计数进行因子分析。包括GLM PCA和NewWave，后者是前者的优化版本，对以上两种方法分别进行模拟数据评估。
方法比较
对包裹了来自相同RNA的等分试样的液滴均匀溶液，分别比较原始计数、delta、剪切Pearson残差、Sanity MAP等转化方法后PCA前两个主成分分布、均值-方差关系分布，图1a中每个点代表一个液滴。从分析结果可以看出，对于delta方法，放缩因子仍代表强方差变异，不同级别放缩因子液滴明显分群，这也是delta方法的致命缺点，将具有较大计数和较大放缩因子的液滴和较小计数和较小放缩因子的液滴分到了一起，而这并不符合负二项分布的均值-关系分布。相比之下，GLM残差和Sanity MAP都能够对液滴进行很好的混合。
对取自造血细胞的10x数据2597个基因表达分析结果显示，delta方法还有一个缺陷是对低表达基因(表达均值<0.1)转化后方差几乎为0，不能体现低表达基因对结果的影响。相比之下，经过GLM残差转化的方法显示低表达基因对均值的弱变化，除非极低表达基因，而Sanity MAP转化后表达是多样的，显示这种方法不直接与稳定方差相关。
SFTPC基因在肺上皮细胞中呈现双峰表达模式，挑选该基因的单细胞表达分布来衡量3类转化方法差异。结果显示：delta和Sanity MAP方法均能很好的体现SFTPC的不同表达水平细胞亚群水平，而Pearson残差对于高表达细胞亚群的方差相比低表达细胞亚群的方差不能更大程度的缩小，所以高表达亚群细胞会被更加细分，而低表达细胞亚群会更加聚集。这会影响此类在不同细胞中差异表达基因的可视化，还有标记基因检测和细胞聚类和分类。
图1.不同转化方法对比
基准数据测评
方法评估取决于分析目的，如果为了比较细胞类型特异性标记基因检测，可通过图1c看细胞计数分布，而本文评估目标定位于基于低维数学结构（聚类、轨迹、分支和组合）理解细胞类型和状态的多样性，k-NN图体现了基本功能结构，通过计算k-NN图之间的重叠度以评估结果的一致性，值越大说明重叠度越高，准确性越好。
本文主要通过以下方法构建真实数据集：
从GEO数据库中下载10个10X转录组数据集，每份数据一分为二，理论上两份数据之间k-NN图重叠度较好，一致性较高。
基于已发表的4种数据模拟框架生成模拟数据，且本文对其中一种模拟方法进行了修正，对模拟数据不同转化方法生成k-NN图与参考k-NN图比较重叠度。
对来自mcSCRB和Smart-seq3的5个深度测序数据集，对其进行下采样到10x深度，分别比较深度测序k-NN图与下采样数据k-NN图的重叠度。
对以上评估数据集分别进行PCA分析、k-NN图绘制、k-NN与参考k-NN重叠度计算，涉及包含在四大类中的22种转化方法，包括过度离散值固定为0、0.05和基因特异性评估，4-8个PCA亚群参数设置，k=10、50、100的最佳邻接细胞重叠计算。对k=50的最佳邻接图形重叠度计算结果显示：在三种数据集中，delta方法都得到了较高的k-NN图重叠，特别在一分为二的数据集上体现比较明显，而在模拟框架数据集上，不同转化方法差异不明显，但针对不同的模拟方法还是有很大的差异。下采样数据集代表了真实数据情况，更接近于10x转录组数据，结果显示不同转化方法k-NN重叠度差异不大。
之前文章研究发现，在random walk模拟方法产生的数据集上，Sanity是识别细胞k-NN的最佳方法，本文研究人员在先将多维单细胞数据降到低维，然后利用delta和残差转化方法分析细胞k-NN时取得同样优异表现。通过PCA亚群数量和k-NN重叠度分布曲线可以看出，数据空间维度大小对性能评估影响很大，随着空间维度增大（>1000），k-NN重叠度下降，特别是在两种模拟方法数据集上体现比较明显，在维度较低(<100)的一致性评估数据集和下采样数据集中也有所体现。
图2.基准数据测试结果
对10,064个辅助T细胞的10x转录组测序数据进行方差稳定转化和k-NN计算，不同方法耗费的CPU时间差距比较大，基于表达推测的方法整体比较耗时，其中Sanity distance时间最长，由于要计算细胞之间的欧式距离，GLM PCA和NewWave耗时次之，delta方法最快，保留了数据的稀疏性，除了Pearson残差，其它残差方法速度跟delta方法一样快。随着细胞个数增加，所有方法的消耗时间同步增加，大部分方法的增长曲线斜率为1，而Sanity和GLM PCA的增长斜率为1.5。
图3.计算资源消耗统计
所有转化方法综合对比
在细胞亚群分类方面，log(y/s+1)转化始终优于其它方法，而且当PCA降维到合适的目标维度后，log(y/s+1)转化在所有三个基准中都比更复杂的基于表达推断的转换表现得更好。额外的处理步骤，如重新缩放移位对数的输出，选择hvg或使用zscore使所有基因的方差相等等并没有太大改善识别细胞k-NN结果。对于皮Pearson残差和acosh变换，设置α > 0是有用的，但与使用固定值如0.05相比从输入数据估计α值，性能没有明显的提高。随着每个细胞测序深度增加，所有方法通常都有更好的k-NN重叠度，显示出更好的分析性能，主要是由于随着测序深度增加，随机采样带来的噪音就会减小。
图4.所有转化方法评估结果
总结
本文从三类基准数据，包括10个单细胞测序数据一分为二的数据集、4种不同模拟方法产生的数据集以及从高深度单细胞RNA测序数据下采样得到深度10的数据集，评估来自于4大类（delta、残差、表达推测、因子分析）共22种细分计数转化方法对细胞亚群结构的分析准确度，研究发现采用简单的log(y/s+1)转化，加上数据降维到低维度空间就可以跟有着较好理论基础的复杂模型一样的优异表现。
