不同表达矩阵对筛选差异长链非编码RNA的影响

魏 豪，邱家俊，颜景斌

上海市儿童医院，上海交通大学医学院附属儿童医院医学遗传研究所，上海市胚胎与生殖工程重点实验室，上海 200040

长 链 非 编 码RNA （long non-coding RNA，lncRNA）是指一类核苷酸长度大于200nt 的非编码RNA，与其他的非编码RNA 曾一度被认为是“转录噪声”［1］。随着第二代测序技术的发展，人们发现在许多模式生物以及人类的基因组中，lncRNA 均有广泛的表达。同时越来越多的功能研究认为lncRNA 可以参与调节多种重要的细胞活动，如基因表达、招募染色质修饰物、调节X染色体失活、基因组印记、蛋白质折叠和蛋白质活性等［2］。

以第二代测序技术为基础的全转录组测序逐渐成为研究lncRNA 表达和功能的重要手段。通过对转录组测序数据的处理和分析，人们可以更加具体地在转录本异构体、核苷酸变异、转录后碱基修饰等方面进行研究。其中的一个主要环节就是对不同组别之间的测序数据进行差异表达分析，表达水平的准确估计对这个过程至关重要［3］。

在转录组测序中，计算一个转录本的表达水平往往要考虑到测序深度、基因长度等因素。起初人们选择使用每百万reads 中每1 000 碱基长度的reads 数（reads per kilobase per million mapped reads，RPKM）等标准化数值来表示相对表达水平，以减少非研究因素的影响，但仍有一些不足，例如在差异分析时部分基因表达水平过高就会产生大量的假阳性结果［4-6］。此外人们还开发出许多R 语言软件包，例如DESeq2［7］、edgeR［8］、PoissonSeq［9］、CuffDiff［10］等，这些软件往往采用自研的一套算法直接对原始表达数据进行标准化处理，同时分析差异基因表达水平。

在分析全转录组测序中lncRNA 的差异表达水平时，人们往往直接对全转录组表达矩阵进行差异分析，再从差异分析结果中将lncRNA 分选出来。但许多研究已经表明lncRNA 的表达水平普遍低于编码蛋白质的mRNA［1］，这就意味着lncRNA 和mRNA 的表达水平分布或许也不相同，以此为背景对全部RNA的表达数据进行标准化和差异分析时，就可能导致一部分没有表达差异的lncRNA 被认为具有差异，产生假阳性结果。为了减少差异分析时各种背景因素的影响，人们选择使用DESeq2 或edgeR 中的标准化模型进行分析。已有不少研究表示，DESeq2 和edgeR 差异表达分析的假阳性率（false positive rate，FPR）更低［11］。因此，选择合适的分析方法对于筛选差异的lncRNA至关重要。

本研究对全转录组测序中仅具有lncRNA 的表达矩阵进行差异表达分析，并将之与总体RNA 表达矩阵的分析结果作比较，从而探究出在筛选差异lncRNA方面误差更小的方法。

1 数据与方法

1.1 数据来源

从NCBI_GEO 数据库（https：//www.ncbi.nlm.nih.gov/geo）中下载2个不同来源的人类全转录组测序数据（GSE49712），共10个样本。A组为5个人类通用参考RNA（universal human reference RNA，UHRR）样本，每个样本掺入了2%体积的外源RNA对照物联盟（external RNA control consortium，ERCC）混合物1；B组为5个人脑参考RNA（human brain reference RNA，HBRR）样本，每个样本掺入了2%体积的ERCC混合物2。ERCC spike-in 对照是由92 种含polyA 序列的人工合成的外源寡核苷酸组成，它们在A组和B组之间的掺入比例有以下4种：1∶2、2∶3、1∶1和4∶1。

从赛默飞网站（https：//www.thermofisher.cn）中下载ERCC spike-in 对照的fastq 文件和GTF 文件。从GENCODE 网站（https：//www.gencodegenes.org）中下载人类GRCh37 完全注释基因组以及相应的lncRNA注释基因组。

1.2 获取表达矩阵

在获得原始SRA 格式的数据后，对其进行格式转换、质控（去除低质量reads和接头）、比对和计数等处理，具体步骤如下。

（1）使用sratoolkit（version 2.10.8）程序将原始数据转换为fastq格式。

（2）使用trim_galore（version 0.6.6）程序过滤掉低质量reads；使用FastQC程序查看质控结果。

（3）使用hisat2（version 2.2.1）程序中的hisat2-build 命令构建含有spike-in 序列的参考基因组，然后使用hisat2命令对reads进行比对。

（4）使用samtools （version 1.7）程序对获得的二进制sam文件进行排序、建立索引等处理。

（5）使用awk 命令从人类GRCh37 完全注释基因组中提取出仅含有编码蛋白基因的注释基因组，再将spike-in的注释信息分别加入各个注释基因组中。

（6）使用featureCounts（version 2.0.1）程序对比对后的片段进行计数和注释，分别获取总体RNA 表达矩阵、mRNA表达矩阵和lncRNA表达矩阵。

1.3 表达水平差异分析

分别使用DESeq2 和edgeR 软件包对A 组和B 组之间的差异表达水平进行分析，获取以3 种表达矩阵为背景的差异spik-in 信息。使用DESeq2 分别对A 组和B组组内的差异表达水平进行分析，获取总体RNA表达矩阵和lncRNA表达矩阵的差异lncRNA信息。

本研究要比较的2 种分析方法分别为使用总体RNA 表达矩阵筛选差异lncRNA 的方法（简称总体RNA 法）和使用lncRNA 表达矩阵筛选差异lncRNA的方法（简称lncRNA法）。

1.4 样本层级聚类分析

分别过滤掉每个表达矩阵中没有表达的基因，使用R 语言中的hclust 功能对各个表达矩阵中所有样本进行层级聚类分析。

1.5 受试者操作特征曲线（ROC）分析

使用pROC 软件包对3个表达矩阵的差异spike-in信息进行分析，获取在不同P值下的FPR和真阳性率（true positive ratio，TPR），并计算其曲线下面积（area under the curve，AUC）值。其中在A 组和B 组之间掺入比例为1 的可认定为没有差异，其余3 种比例则认定为有差异。

1.6 统计学方法

对不同方法获取差异spike-in 信息效果的评价采用ROC 曲线，使用AUC 值预测特异性和准确性。部分所得数据使用Origin 2022 软件进行统计学分析并作图。P＜0.05表示差异有统计学意义。

2 结果

2.1 表达矩阵的获取

本研究对2个组别共10个样本的测序数据进行处理之后，使用不同的注释基因组对其计数，获得了3个表达矩阵（总体RNA、mRNA、lncRNA）。经统计，A 组样本和B 组样本中分别有47 583 个和47 177 个已知基因得到了表达，其中A 组样本包括19 061 个mRNA 基因和15 711 个lncRNA 基因，B 组样本包括18 884个mRNA基因和15 652个lncRNA基因。A组样本中mRNA 和lncRNA 总体表达水平的平均占比分别为95.58%和3.12%，B 组样本中mRNA 和lncRNA 总体表达水平的平均占比分别为93.87%和3.71%。我们对各个表达矩阵中的样本进行了层级聚类分析，结果表明所有表达矩阵均能很好地分辨样本差异（图1）。

图1 A组和B组样本表达矩阵的层级聚类图Fig 1 Hierarchical clustering of the expression profiles from sample A and B

2.2 Spike-in表达水平差异分析

使用DESeq2软件包对3个表达矩阵进行差异表达分析，从中提取出spike-in在A组和B组之间的差异表达信息。为了确认在广泛采用的差异基因筛选标准（P＜0.05） 下，总体RNA 法和lncRNA 法筛选差异lncRNA的效果，我们统计了spike-in在总体RNA表达矩阵和lncRNA 表达矩阵中假阳性和假阴性的数目（表1），并计算了FPR 和假阴性率（false negative rate，FNR）（图2A）。结果发现，使用总体RNA法分析的spike-in FPR 为0.52（12/23），而使用lncRNA 法分析的spike-in FPR 为0.30（7/23），显然后者筛选差异spike-in的特异性更好。我们对总体RNA 表达矩阵中多出的5个假阳性spike-in进行了进一步分析，除了其中一个spike-in 表达水平极低以外，其余4 个spikein的表达数据在经过标准化处理之后，总体RNA表达矩阵比lncRNA表达矩阵的组间差异更大（图2B）。

图2 Spike-in RNAs差异表达分析的结果Fig 2 Differential expressions of spike-in RNAs

表1 2个表达矩阵预测的差异表达spike-in RNAs数目Tab 1 Number of differential expression spike-in RNAs predicted by two different expression profiles

2.3 ROC分析

为了进一步观察其他P值条件下的差异基因筛选效果，本研究采用spike-in 在不同表达矩阵背景下的P值作ROC 曲线。同时为了排除算法因素的影响，我们分别使用DESeq2 和edgeR 进行ROC 分析（图3），其中横坐标和纵坐标分别用特异性和准确性表示，它们在数值上分别与1-FPR 和TPR 相等。通过每条ROC 曲线下的面积大小，即AUC 值来量化不同分析方法筛选差异spike-in 的效果。在DESeq2 的分析结果中，Spike-in 在3 个表达矩阵中的AUC 值大小关系为AUC（lncRNA）=0.852＞AUC（allRNA）=0.768＞AUC（mRNA）=0.750，而在edgeR 的分析结果中为AUC（lncRNA）=0.878＞AUC（mRNA）=0.798＞AUC（allRNA）=0.787。显然，spike-in 以lncRNA 表达矩阵为背景的AUC值显著高于其他2个表达矩阵。

图3 使用DESeq2和edgeR对spike-in RNAs的ROC分析结果Fig 3 ROC curve of ERCC spike-in RNAs analyzed by DESeq2 and edgeR

2.4 差异表达lncRNA的筛选

由于组内各个样本的测序数据来源完全相同，理论上可以认为组内样本之间基因表达水平没有差异，因此可以用来评估差异分析方法的FPR大小。本研究使用DESeq2 对总体RNA 表达矩阵和lncRNA 表达矩阵中的A 组和B 组分别进行组内差异lncRNA 分析，其分组为A1 和A2、A3 和A4、B1 和B2、B3 和B4。由于组内样本绝大多数基因无表达差异，P值无限接近于1，因此我们在对lncRNA 的P值作密度分布图时，把范围缩小在0～0.2之间。当把筛选差异lncRNA的标准定为P＜0.05 时，在A 组中lncRNA 表达矩阵和总体RNA 表达矩阵的差异lncRNA 数目分别为9 个（占比2.57%）和7 个（占比2.45%），见图4A；B 组中分别15 个（占比3.46%）和17 个（占比3.65%），见图4B。使用不同表达矩阵筛选差异lncRNA 的FPR并没有太大差别，可能与本研究组内样本之间差异基因数目过少有关。

图4 组内差异表达lncRNAs的P值密度曲线Fig 4 Density curve of P-values of differential expressed lncRNAs within groups

2.5 mRNA和lncRNA的表达水平分布

为了进一步探究使用总体RNA 法差异分析的FPR 更高的原因，我们对mRNA 和lncRNA 的表达水平分布进行了分析。鉴于不同基因之间的原始表达水平（raw counts）差距过大，我们对其作对数化处理后进行展示。可以看到在全部样本中mRNA 和lncRNA 各自表达水平的分布基本一致（图5A），而mRNA和lncRNA之间则明显不同（图5B）。

图5 mRNA和lncRNA表达水平的分布图Fig 5 Distribution of mRNA and lncRNA expression levels

3 讨论

近些年来，第二代测序技术以其不断提高的技术水平和持续降低的成本快速发展，以之为基础的转录组测序也成为大部分研究人员进行基因表达研究的主要选择之一［12］。在鉴定不同环境下或不同组织之间样本的差异表达基因方面，相比于之前发展起来的基因芯片等技术，转录组测序有着明显的优势，例如可以对整个基因组实现更高的覆盖率以及可以容易地检测到低表达基因等［13］。但同时转录组测序数据的分析难度也大大提高了，在进行差异表达分析时，除了要考虑基因长度和测序深度外，其他基因的表达水平也是一个重要因素，部分过高表达的基因会影响低表达基因差异程度的评估［5］。而在全转录组测序中大部分lncRNA 表达水平往往明显低于mRNA，因此在筛选差异lncRNA 时，难免会受到其他高表达mRNA的影响。因此，选择合适的筛选方法对于lncRNA 的研究具有重要意义。

自转录组测序诞生以来，围绕差异表达基因分析方法的研究不断涌现［14］，这些研究往往是对差异分析之前的一个关键步骤——数据标准化方法的改进和创新。尽管最初有人认为转录组测序并不需要复杂的标准化过程［15］，但实际分析时总是需要在样本之间进行对比，而原始reads数目受到的非研究因素影响过多，显然不能直接用来比较。最初得到广泛应用的转录组标准化数值为MORTAZAVI 等［16］开发的RPKM。RPKM 值对基因长度和测序深度作了标准化处理，但缺点也很明显，每个基因RPKM值大小的计算很大程度上受到其他基因表达水平的影响，这可能产生实际并不存在的表达差异。随后出现的其他相对表达量均有类似的问题，但这并不影响它们替代原始表达数据被广泛使用。后来人们开发的标准化算法往往与差异表达算法一起整合在软件包中，它们基于大部分基因的表达水平是不具有差异的这一假设，通过各自的方式计算出归一化因子来对数据进行标准化处理。代表性的有DESeq2所使用的RLE（relative log expression）算法［7］和edgeR所使用的TMM（trimmed mean of Mvalues）算法［8］。理论上这2种标准化方式能够较好地减少高表达基因或高变基因的影响。因此，本研究使用DESeq2和edgeR软件包执行差异分析过程，在另一方面也可以来检验这2种标准化方式的实际效果。

本研究从NCBI_GEO数据库中下载了2组来源于测序质量控制联盟（Sequencing Quality Control Consortium，SEQC）的人类参考RNA测序数据［17-18］。由于这2 组数据中包含不同掺入比例的已知浓度的ERCC spike-in对照，以及近1 000个已在qPCR实验中验证表达的基因，常常也被用于转录组测序分析方法的研究。在获得了包含spike-in表达信息的3个表达矩阵后，过滤掉在所有样本中均不表达的基因，我们发现表达mRNA 的基因比表达lncRNA 的基因数目仅多出3 000～4 000 个，而mRNA 在全基因组水平的表达占比上却远远大于lncRNA，说明绝大多数lncRNA的表达水平普遍低于mRNA。对分开的表达矩阵层级聚类分析结果表明，仅有lncRNA 的表达矩阵也能对不同样本作出良好区分，这与已知许多证明lncRNA 更具有组织特异性的研究基本一致［19-20］。

在差异分析后基于每个spike-in 的P值作出的ROC 曲线来看，不管使用哪个软件包执行差异分析，在lncRNA 表达矩阵中的AUC 值都明显大于另外2 个表达矩阵，而由于mRNA 在所有RNA 中表达占比极高，它们的AUC 值大小不相上下。说明在筛选差异lncRNA时，选择使用lncRNA法进行差异分析更接近实际的差异表达情况。我们通过对使用总体RNA 法多出的假阳性spike-in 进一步分析发现，它们的表达数据在总体RNA 表达矩阵背景下的标准化数值组间差异更大，这意味着即使使用DESeq2 或edgeR 等这类在标准化处理时考虑到高表达基因和高变基因的软件包时，也难以完全避免它们的影响。本研究对组内总体RNA 表达矩阵和lncRNA 表达矩阵的差异lncRNA 分析结果表明，在P＜0.05 的筛选条件下，使用2 种方法分析的FPR 非常接近，这与前面筛选差异spike-in 的研究结果不太一致，但仍能说明直接使用lncRNA表达矩阵筛选差异基因的方法是可行的。

值得注意的是，本研究中筛选差异基因的标准均设为P＜0.05，这是因为在很多研究中该标准被看作是差异具有显著性的分水岭［21］。该标准最早由FISHER R在20世纪20年代中期提出，用来描述农业田间试验的显著性［22］。后来该标准应用于医学领域之后，更高的犯错成本使得人们不得不重视它所带来的诸多问题。首先，部分研究者不惜采取P值操纵（P-hacking）的方式来使结果差异更大［23］，这引起了更多研究人员的不满。其次，P＜0.05 的统计结果意味着仍有约5%的概率得到假阳性结果，而高通量测序数据中动辄上万的基因就会产生相当数目的假阳性结果，大量的差异基因也使得FPR比预估的更高［24］。本研究也间接证实了这一点，该标准下筛选的差异spike-in FPR 都比较高。为了降低FPR，许多转录组差异分析软件额外采用了padj（adjustP-value）来表示 显 著 水 平［25］。padj 是 经 过FDR （false discovery rate）矫正后的P值。本研究在对组内样本执行差异分析时也发现，所有基因的padj 均无限接近于1，即没有假阳性。尽管如此，在实际研究过程中P＜0.05仍是表明研究结果具有显著意义的必要条件，但也应注意不能滥用或过于偏信该标准。

我们通过以结果为导向的研究反推出了使用lncRNA表达矩阵筛选效果更好的结论，而对于使用总体RNA表达矩阵筛选差异lncRNA的FPR更高，这一现象背后的统计学原理缺乏深入的剖析。但通过我们对mRNA 和lncRNA 表达水平分布的统计可以看出，二者显然具有不同的分布规律，这或许是导致这种现象的关键原因。由于DESeq2软件包采用的标准化模型较为复杂，仍待后续进一步的数理验证。此外，本研究分析所采用的样本量较少，尤其在组内分析时较难得到有意义的结果，之后我们将继续在样本量更大的数据集中进一步研究和验证。综上，本研究通过对含有ERCC spike-in对照的全转录组测序数据进行差异分析，发现直接使用lncRNA表达矩阵筛选差异lncRNA的方法特异性和准确性更好。这为今后探究多样化的差异lncRNA筛选方法提供了一定的理论依据。

利益冲突声明/Conflict of Interests

所有作者声明不存在利益冲突。

All authors disclose no relevant conflict of interests.

作者贡献/Authors'Contributions

魏豪主要完成数据分析处理与文章撰写的工作，邱家俊主要负责数据集检索与数据分析指导，颜景斌主要负责总体研究思路和论文整体构思。所有作者均阅读并同意了最终稿件的提交。

WEI Hao completed the work of data analysis and wrote the manuscript. QIU Jiajun completed the dataset retrieval and guided the data analysis. YAN Jingbin designed and supervised the overall research ideas and the overall conception of the paper.All the authors have read the last version of paper and consented for submission.

·Received：2022-03-24

·Accepted：2022-07-14

·Published online：2022-07-28