魏 豪，邱家俊，顏景斌

上海市兒童醫(yī)院，上海交通大學(xué)醫(yī)學(xué)院附屬兒童醫(yī)院醫(yī)學(xué)遺傳研究所，上海市胚胎與生殖工程重點實驗室，上海 200040

長 鏈 非 編 碼RNA （long non-coding RNA，lncRNA）是指一類核苷酸長度大于200nt 的非編碼RNA，與其他的非編碼RNA 曾一度被認(rèn)為是“轉(zhuǎn)錄噪聲”［1］。隨著第二代測序技術(shù)的發(fā)展，人們發(fā)現(xiàn)在許多模式生物以及人類的基因組中，lncRNA 均有廣泛的表達(dá)。同時越來越多的功能研究認(rèn)為lncRNA 可以參與調(diào)節(jié)多種重要的細(xì)胞活動，如基因表達(dá)、招募染色質(zhì)修飾物、調(diào)節(jié)X染色體失活、基因組印記、蛋白質(zhì)折疊和蛋白質(zhì)活性等［2］。

以第二代測序技術(shù)為基礎(chǔ)的全轉(zhuǎn)錄組測序逐漸成為研究lncRNA 表達(dá)和功能的重要手段。通過對轉(zhuǎn)錄組測序數(shù)據(jù)的處理和分析，人們可以更加具體地在轉(zhuǎn)錄本異構(gòu)體、核苷酸變異、轉(zhuǎn)錄后堿基修飾等方面進(jìn)行研究。其中的一個主要環(huán)節(jié)就是對不同組別之間的測序數(shù)據(jù)進(jìn)行差異表達(dá)分析，表達(dá)水平的準(zhǔn)確估計對這個過程至關(guān)重要［3］。

在轉(zhuǎn)錄組測序中，計算一個轉(zhuǎn)錄本的表達(dá)水平往往要考慮到測序深度、基因長度等因素。起初人們選擇使用每百萬reads 中每1 000 堿基長度的reads 數(shù)（reads per kilobase per million mapped reads，RPKM）等標(biāo)準(zhǔn)化數(shù)值來表示相對表達(dá)水平，以減少非研究因素的影響，但仍有一些不足，例如在差異分析時部分基因表達(dá)水平過高就會產(chǎn)生大量的假陽性結(jié)果［4-6］。此外人們還開發(fā)出許多R 語言軟件包，例如DESeq2［7］、edgeR［8］、PoissonSeq［9］、CuffDiff［10］等，這些軟件往往采用自研的一套算法直接對原始表達(dá)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，同時分析差異基因表達(dá)水平。

在分析全轉(zhuǎn)錄組測序中l(wèi)ncRNA 的差異表達(dá)水平時，人們往往直接對全轉(zhuǎn)錄組表達(dá)矩陣進(jìn)行差異分析，再從差異分析結(jié)果中將lncRNA 分選出來。但許多研究已經(jīng)表明lncRNA 的表達(dá)水平普遍低于編碼蛋白質(zhì)的mRNA［1］，這就意味著lncRNA 和mRNA 的表達(dá)水平分布或許也不相同，以此為背景對全部RNA的表達(dá)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化和差異分析時，就可能導(dǎo)致一部分沒有表達(dá)差異的lncRNA 被認(rèn)為具有差異，產(chǎn)生假陽性結(jié)果。為了減少差異分析時各種背景因素的影響，人們選擇使用DESeq2 或edgeR 中的標(biāo)準(zhǔn)化模型進(jìn)行分析。已有不少研究表示，DESeq2 和edgeR 差異表達(dá)分析的假陽性率（false positive rate，F(xiàn)PR）更低［11］。因此，選擇合適的分析方法對于篩選差異的lncRNA至關(guān)重要。

本研究對全轉(zhuǎn)錄組測序中僅具有l(wèi)ncRNA 的表達(dá)矩陣進(jìn)行差異表達(dá)分析，并將之與總體RNA 表達(dá)矩陣的分析結(jié)果作比較，從而探究出在篩選差異lncRNA方面誤差更小的方法。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源

從NCBI_GEO 數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo）中下載2個不同來源的人類全轉(zhuǎn)錄組測序數(shù)據(jù)（GSE49712），共10個樣本。A組為5個人類通用參考RNA（universal human reference RNA，UHRR）樣本，每個樣本摻入了2%體積的外源RNA對照物聯(lián)盟（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。

從賽默飛網(wǎng)站（https：//www.thermofisher.cn）中下載ERCC spike-in 對照的fastq 文件和GTF 文件。從GENCODE 網(wǎng)站（https：//www.gencodegenes.org）中下載人類GRCh37 完全注釋基因組以及相應(yīng)的lncRNA注釋基因組。

1.2 獲取表達(dá)矩陣

在獲得原始SRA 格式的數(shù)據(jù)后，對其進(jìn)行格式轉(zhuǎn)換、質(zhì)控（去除低質(zhì)量reads和接頭）、比對和計數(shù)等處理，具體步驟如下。

（1）使用sratoolkit（version 2.10.8）程序?qū)⒃紨?shù)據(jù)轉(zhuǎn)換為fastq格式。

（2）使用trim_galore（version 0.6.6）程序過濾掉低質(zhì)量reads；使用FastQC程序查看質(zhì)控結(jié)果。

（3）使用hisat2（version 2.2.1）程序中的hisat2-build 命令構(gòu)建含有spike-in 序列的參考基因組，然后使用hisat2命令對reads進(jìn)行比對。

（4）使用samtools （version 1.7）程序?qū)Λ@得的二進(jìn)制sam文件進(jìn)行排序、建立索引等處理。

（5）使用awk 命令從人類GRCh37 完全注釋基因組中提取出僅含有編碼蛋白基因的注釋基因組，再將spike-in的注釋信息分別加入各個注釋基因組中。

（6）使用featureCounts（version 2.0.1）程序?qū)Ρ葘蟮钠芜M(jìn)行計數(shù)和注釋，分別獲取總體RNA 表達(dá)矩陣、mRNA表達(dá)矩陣和lncRNA表達(dá)矩陣。

1.3 表達(dá)水平差異分析

分別使用DESeq2 和edgeR 軟件包對A 組和B 組之間的差異表達(dá)水平進(jìn)行分析，獲取以3 種表達(dá)矩陣為背景的差異spik-in 信息。使用DESeq2 分別對A 組和B組組內(nèi)的差異表達(dá)水平進(jìn)行分析，獲取總體RNA表達(dá)矩陣和lncRNA表達(dá)矩陣的差異lncRNA信息。

本研究要比較的2 種分析方法分別為使用總體RNA 表達(dá)矩陣篩選差異lncRNA 的方法（簡稱總體RNA 法）和使用lncRNA 表達(dá)矩陣篩選差異lncRNA的方法（簡稱lncRNA法）。

1.4 樣本層級聚類分析

分別過濾掉每個表達(dá)矩陣中沒有表達(dá)的基因，使用R 語言中的hclust 功能對各個表達(dá)矩陣中所有樣本進(jìn)行層級聚類分析。

1.5 受試者操作特征曲線（ROC）分析

使用pROC 軟件包對3個表達(dá)矩陣的差異spike-in信息進(jìn)行分析，獲取在不同P值下的FPR和真陽性率（true positive ratio，TPR），并計算其曲線下面積（area under the curve，AUC）值。其中在A 組和B 組之間摻入比例為1 的可認(rèn)定為沒有差異，其余3 種比例則認(rèn)定為有差異。

1.6 統(tǒng)計學(xué)方法

對不同方法獲取差異spike-in 信息效果的評價采用ROC 曲線，使用AUC 值預(yù)測特異性和準(zhǔn)確性。部分所得數(shù)據(jù)使用Origin 2022 軟件進(jìn)行統(tǒng)計學(xué)分析并作圖。P＜0.05表示差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 表達(dá)矩陣的獲取

本研究對2個組別共10個樣本的測序數(shù)據(jù)進(jìn)行處理之后，使用不同的注釋基因組對其計數(shù)，獲得了3個表達(dá)矩陣（總體RNA、mRNA、lncRNA）。經(jīng)統(tǒng)計，A 組樣本和B 組樣本中分別有47 583 個和47 177 個已知基因得到了表達(dá)，其中A 組樣本包括19 061 個mRNA 基因和15 711 個lncRNA 基因，B 組樣本包括18 884個mRNA基因和15 652個lncRNA基因。A組樣本中mRNA 和lncRNA 總體表達(dá)水平的平均占比分別為95.58%和3.12%，B 組樣本中mRNA 和lncRNA 總體表達(dá)水平的平均占比分別為93.87%和3.71%。我們對各個表達(dá)矩陣中的樣本進(jìn)行了層級聚類分析，結(jié)果表明所有表達(dá)矩陣均能很好地分辨樣本差異（圖1）。

圖1 A組和B組樣本表達(dá)矩陣的層級聚類圖Fig 1 Hierarchical clustering of the expression profiles from sample A and B

2.2 Spike-in表達(dá)水平差異分析

使用DESeq2軟件包對3個表達(dá)矩陣進(jìn)行差異表達(dá)分析，從中提取出spike-in在A組和B組之間的差異表達(dá)信息。為了確認(rèn)在廣泛采用的差異基因篩選標(biāo)準(zhǔn)（P＜0.05） 下，總體RNA 法和lncRNA 法篩選差異lncRNA的效果，我們統(tǒng)計了spike-in在總體RNA表達(dá)矩陣和lncRNA 表達(dá)矩陣中假陽性和假陰性的數(shù)目（表1），并計算了FPR 和假陰性率（false negative rate，F(xiàn)NR）（圖2A）。結(jié)果發(fā)現(xiàn)，使用總體RNA法分析的spike-in FPR 為0.52（12/23），而使用lncRNA 法分析的spike-in FPR 為0.30（7/23），顯然后者篩選差異spike-in的特異性更好。我們對總體RNA 表達(dá)矩陣中多出的5個假陽性spike-in進(jìn)行了進(jìn)一步分析，除了其中一個spike-in 表達(dá)水平極低以外，其余4 個spikein的表達(dá)數(shù)據(jù)在經(jīng)過標(biāo)準(zhǔn)化處理之后，總體RNA表達(dá)矩陣比lncRNA表達(dá)矩陣的組間差異更大（圖2B）。

圖2 Spike-in RNAs差異表達(dá)分析的結(jié)果Fig 2 Differential expressions of spike-in RNAs

表1 2個表達(dá)矩陣預(yù)測的差異表達(dá)spike-in RNAs數(shù)目Tab 1 Number of differential expression spike-in RNAs predicted by two different expression profiles

2.3 ROC分析

為了進(jìn)一步觀察其他P值條件下的差異基因篩選效果，本研究采用spike-in 在不同表達(dá)矩陣背景下的P值作ROC 曲線。同時為了排除算法因素的影響，我們分別使用DESeq2 和edgeR 進(jìn)行ROC 分析（圖3），其中橫坐標(biāo)和縱坐標(biāo)分別用特異性和準(zhǔn)確性表示，它們在數(shù)值上分別與1-FPR 和TPR 相等。通過每條ROC 曲線下的面積大小，即AUC 值來量化不同分析方法篩選差異spike-in 的效果。在DESeq2 的分析結(jié)果中，Spike-in 在3 個表達(dá)矩陣中的AUC 值大小關(guān)系為AUC（lncRNA）=0.852＞AUC（allRNA）=0.768＞AUC（mRNA）=0.750，而在edgeR 的分析結(jié)果中為AUC（lncRNA）=0.878＞AUC（mRNA）=0.798＞AUC（allRNA）=0.787。顯然，spike-in 以lncRNA 表達(dá)矩陣為背景的AUC值顯著高于其他2個表達(dá)矩陣。

圖3 使用DESeq2和edgeR對spike-in RNAs的ROC分析結(jié)果Fig 3 ROC curve of ERCC spike-in RNAs analyzed by DESeq2 and edgeR

2.4 差異表達(dá)lncRNA的篩選

由于組內(nèi)各個樣本的測序數(shù)據(jù)來源完全相同，理論上可以認(rèn)為組內(nèi)樣本之間基因表達(dá)水平?jīng)]有差異，因此可以用來評估差異分析方法的FPR大小。本研究使用DESeq2 對總體RNA 表達(dá)矩陣和lncRNA 表達(dá)矩陣中的A 組和B 組分別進(jìn)行組內(nèi)差異lncRNA 分析，其分組為A1 和A2、A3 和A4、B1 和B2、B3 和B4。由于組內(nèi)樣本絕大多數(shù)基因無表達(dá)差異，P值無限接近于1，因此我們在對lncRNA 的P值作密度分布圖時，把范圍縮小在0～0.2之間。當(dāng)把篩選差異lncRNA的標(biāo)準(zhǔn)定為P＜0.05 時，在A 組中l(wèi)ncRNA 表達(dá)矩陣和總體RNA 表達(dá)矩陣的差異lncRNA 數(shù)目分別為9 個（占比2.57%）和7 個（占比2.45%），見圖4A；B 組中分別15 個（占比3.46%）和17 個（占比3.65%），見圖4B。使用不同表達(dá)矩陣篩選差異lncRNA 的FPR并沒有太大差別，可能與本研究組內(nèi)樣本之間差異基因數(shù)目過少有關(guān)。

圖4 組內(nèi)差異表達(dá)lncRNAs的P值密度曲線Fig 4 Density curve of P-values of differential expressed lncRNAs within groups

2.5 mRNA和lncRNA的表達(dá)水平分布

為了進(jìn)一步探究使用總體RNA 法差異分析的FPR 更高的原因，我們對mRNA 和lncRNA 的表達(dá)水平分布進(jìn)行了分析。鑒于不同基因之間的原始表達(dá)水平（raw counts）差距過大，我們對其作對數(shù)化處理后進(jìn)行展示。可以看到在全部樣本中mRNA 和lncRNA 各自表達(dá)水平的分布基本一致（圖5A），而mRNA和lncRNA之間則明顯不同（圖5B）。

圖5 mRNA和lncRNA表達(dá)水平的分布圖Fig 5 Distribution of mRNA and lncRNA expression levels

3 討論

近些年來，第二代測序技術(shù)以其不斷提高的技術(shù)水平和持續(xù)降低的成本快速發(fā)展，以之為基礎(chǔ)的轉(zhuǎn)錄組測序也成為大部分研究人員進(jìn)行基因表達(dá)研究的主要選擇之一［12］。在鑒定不同環(huán)境下或不同組織之間樣本的差異表達(dá)基因方面，相比于之前發(fā)展起來的基因芯片等技術(shù)，轉(zhuǎn)錄組測序有著明顯的優(yōu)勢，例如可以對整個基因組實現(xiàn)更高的覆蓋率以及可以容易地檢測到低表達(dá)基因等［13］。但同時轉(zhuǎn)錄組測序數(shù)據(jù)的分析難度也大大提高了，在進(jìn)行差異表達(dá)分析時，除了要考慮基因長度和測序深度外，其他基因的表達(dá)水平也是一個重要因素，部分過高表達(dá)的基因會影響低表達(dá)基因差異程度的評估［5］。而在全轉(zhuǎn)錄組測序中大部分lncRNA 表達(dá)水平往往明顯低于mRNA，因此在篩選差異lncRNA 時，難免會受到其他高表達(dá)mRNA的影響。因此，選擇合適的篩選方法對于lncRNA 的研究具有重要意義。

自轉(zhuǎn)錄組測序誕生以來，圍繞差異表達(dá)基因分析方法的研究不斷涌現(xiàn)［14］，這些研究往往是對差異分析之前的一個關(guān)鍵步驟——數(shù)據(jù)標(biāo)準(zhǔn)化方法的改進(jìn)和創(chuàng)新。盡管最初有人認(rèn)為轉(zhuǎn)錄組測序并不需要復(fù)雜的標(biāo)準(zhǔn)化過程［15］，但實際分析時總是需要在樣本之間進(jìn)行對比，而原始reads數(shù)目受到的非研究因素影響過多，顯然不能直接用來比較。最初得到廣泛應(yīng)用的轉(zhuǎn)錄組標(biāo)準(zhǔn)化數(shù)值為MORTAZAVI 等［16］開發(fā)的RPKM。RPKM 值對基因長度和測序深度作了標(biāo)準(zhǔn)化處理，但缺點也很明顯，每個基因RPKM值大小的計算很大程度上受到其他基因表達(dá)水平的影響，這可能產(chǎn)生實際并不存在的表達(dá)差異。隨后出現(xiàn)的其他相對表達(dá)量均有類似的問題，但這并不影響它們替代原始表達(dá)數(shù)據(jù)被廣泛使用。后來人們開發(fā)的標(biāo)準(zhǔn)化算法往往與差異表達(dá)算法一起整合在軟件包中，它們基于大部分基因的表達(dá)水平是不具有差異的這一假設(shè)，通過各自的方式計算出歸一化因子來對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理。代表性的有DESeq2所使用的RLE（relative log expression）算法［7］和edgeR所使用的TMM（trimmed mean of Mvalues）算法［8］。理論上這2種標(biāo)準(zhǔn)化方式能夠較好地減少高表達(dá)基因或高變基因的影響。因此，本研究使用DESeq2和edgeR軟件包執(zhí)行差異分析過程，在另一方面也可以來檢驗這2種標(biāo)準(zhǔn)化方式的實際效果。

本研究從NCBI_GEO數(shù)據(jù)庫中下載了2組來源于測序質(zhì)量控制聯(lián)盟（Sequencing Quality Control Consortium，SEQC）的人類參考RNA測序數(shù)據(jù)［17-18］。由于這2 組數(shù)據(jù)中包含不同摻入比例的已知濃度的ERCC spike-in對照，以及近1 000個已在qPCR實驗中驗證表達(dá)的基因，常常也被用于轉(zhuǎn)錄組測序分析方法的研究。在獲得了包含spike-in表達(dá)信息的3個表達(dá)矩陣后，過濾掉在所有樣本中均不表達(dá)的基因，我們發(fā)現(xiàn)表達(dá)mRNA 的基因比表達(dá)lncRNA 的基因數(shù)目僅多出3 000～4 000 個，而mRNA 在全基因組水平的表達(dá)占比上卻遠(yuǎn)遠(yuǎn)大于lncRNA，說明絕大多數(shù)lncRNA的表達(dá)水平普遍低于mRNA。對分開的表達(dá)矩陣層級聚類分析結(jié)果表明，僅有l(wèi)ncRNA 的表達(dá)矩陣也能對不同樣本作出良好區(qū)分，這與已知許多證明lncRNA 更具有組織特異性的研究基本一致［19-20］。

在差異分析后基于每個spike-in 的P值作出的ROC 曲線來看，不管使用哪個軟件包執(zhí)行差異分析，在lncRNA 表達(dá)矩陣中的AUC 值都明顯大于另外2 個表達(dá)矩陣，而由于mRNA 在所有RNA 中表達(dá)占比極高，它們的AUC 值大小不相上下。說明在篩選差異lncRNA時，選擇使用lncRNA法進(jìn)行差異分析更接近實際的差異表達(dá)情況。我們通過對使用總體RNA 法多出的假陽性spike-in 進(jìn)一步分析發(fā)現(xiàn)，它們的表達(dá)數(shù)據(jù)在總體RNA 表達(dá)矩陣背景下的標(biāo)準(zhǔn)化數(shù)值組間差異更大，這意味著即使使用DESeq2 或edgeR 等這類在標(biāo)準(zhǔn)化處理時考慮到高表達(dá)基因和高變基因的軟件包時，也難以完全避免它們的影響。本研究對組內(nèi)總體RNA 表達(dá)矩陣和lncRNA 表達(dá)矩陣的差異lncRNA 分析結(jié)果表明，在P＜0.05 的篩選條件下，使用2 種方法分析的FPR 非常接近，這與前面篩選差異spike-in 的研究結(jié)果不太一致，但仍能說明直接使用lncRNA表達(dá)矩陣篩選差異基因的方法是可行的。

值得注意的是，本研究中篩選差異基因的標(biāo)準(zhǔn)均設(shè)為P＜0.05，這是因為在很多研究中該標(biāo)準(zhǔn)被看作是差異具有顯著性的分水嶺［21］。該標(biāo)準(zhǔn)最早由FISHER R在20世紀(jì)20年代中期提出，用來描述農(nóng)業(yè)田間試驗的顯著性［22］。后來該標(biāo)準(zhǔn)應(yīng)用于醫(yī)學(xué)領(lǐng)域之后，更高的犯錯成本使得人們不得不重視它所帶來的諸多問題。首先，部分研究者不惜采取P值操縱（P-hacking）的方式來使結(jié)果差異更大［23］，這引起了更多研究人員的不滿。其次，P＜0.05 的統(tǒng)計結(jié)果意味著仍有約5%的概率得到假陽性結(jié)果，而高通量測序數(shù)據(jù)中動輒上萬的基因就會產(chǎn)生相當(dāng)數(shù)目的假陽性結(jié)果，大量的差異基因也使得FPR比預(yù)估的更高［24］。本研究也間接證實了這一點，該標(biāo)準(zhǔn)下篩選的差異spike-in FPR 都比較高。為了降低FPR，許多轉(zhuǎn)錄組差異分析軟件額外采用了padj（adjustP-value）來表示 顯 著 水 平［25］。padj 是 經(jīng) 過FDR （false discovery rate）矯正后的P值。本研究在對組內(nèi)樣本執(zhí)行差異分析時也發(fā)現(xiàn)，所有基因的padj 均無限接近于1，即沒有假陽性。盡管如此，在實際研究過程中P＜0.05仍是表明研究結(jié)果具有顯著意義的必要條件，但也應(yīng)注意不能濫用或過于偏信該標(biāo)準(zhǔn)。

我們通過以結(jié)果為導(dǎo)向的研究反推出了使用lncRNA表達(dá)矩陣篩選效果更好的結(jié)論，而對于使用總體RNA表達(dá)矩陣篩選差異lncRNA的FPR更高，這一現(xiàn)象背后的統(tǒng)計學(xué)原理缺乏深入的剖析。但通過我們對mRNA 和lncRNA 表達(dá)水平分布的統(tǒng)計可以看出，二者顯然具有不同的分布規(guī)律，這或許是導(dǎo)致這種現(xiàn)象的關(guān)鍵原因。由于DESeq2軟件包采用的標(biāo)準(zhǔn)化模型較為復(fù)雜，仍待后續(xù)進(jìn)一步的數(shù)理驗證。此外，本研究分析所采用的樣本量較少，尤其在組內(nèi)分析時較難得到有意義的結(jié)果，之后我們將繼續(xù)在樣本量更大的數(shù)據(jù)集中進(jìn)一步研究和驗證。綜上，本研究通過對含有ERCC spike-in對照的全轉(zhuǎn)錄組測序數(shù)據(jù)進(jìn)行差異分析，發(fā)現(xiàn)直接使用lncRNA表達(dá)矩陣篩選差異lncRNA的方法特異性和準(zhǔn)確性更好。這為今后探究多樣化的差異lncRNA篩選方法提供了一定的理論依據(jù)。

利益沖突聲明/Conflict of Interests

所有作者聲明不存在利益沖突。

All authors disclose no relevant conflict of interests.

作者貢獻(xiàn)/Authors'Contributions

魏豪主要完成數(shù)據(jù)分析處理與文章撰寫的工作，邱家俊主要負(fù)責(zé)數(shù)據(jù)集檢索與數(shù)據(jù)分析指導(dǎo)，顏景斌主要負(fù)責(zé)總體研究思路和論文整體構(gòu)思。所有作者均閱讀并同意了最終稿件的提交。

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