魏 豪，邱家俊，顏景斌

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

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

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

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

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

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

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

1.1 數(shù)據(jù)來(lái)源

從NCBI_GEO 數(shù)據(jù)庫(kù)（https：//www.ncbi.nlm.nih.gov/geo）中下載2個(gè)不同來(lái)源的人類(lèi)全轉(zhuǎn)錄組測(cè)序數(shù)據(jù)（GSE49712），共10個(gè)樣本。A組為5個(gè)人類(lèi)通用參考RNA（universal human reference RNA，UHRR）樣本，每個(gè)樣本摻入了2%體積的外源RNA對(duì)照物聯(lián)盟（external RNA control consortium，ERCC）混合物1；B組為5個(gè)人腦參考RNA（human brain reference RNA，HBRR）樣本，每個(gè)樣本摻入了2%體積的ERCC混合物2。ERCC spike-in 對(duì)照是由92 種含polyA 序列的人工合成的外源寡核苷酸組成，它們?cè)贏組和B組之間的摻入比例有以下4種：1∶2、2∶3、1∶1和4∶1。

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

1.2 獲取表達(dá)矩陣

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

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

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

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

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

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

（6）使用featureCounts（version 2.0.1）程序?qū)Ρ葘?duì)后的片段進(jìn)行計(jì)數(shù)和注釋?zhuān)謩e獲取總體RNA 表達(dá)矩陣、mRNA表達(dá)矩陣和lncRNA表達(dá)矩陣。

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

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

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

1.4 樣本層級(jí)聚類(lèi)分析

分別過(guò)濾掉每個(gè)表達(dá)矩陣中沒(méi)有表達(dá)的基因，使用R 語(yǔ)言中的hclust 功能對(duì)各個(gè)表達(dá)矩陣中所有樣本進(jìn)行層級(jí)聚類(lèi)分析。

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

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

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

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

2 結(jié)果

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

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

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

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

使用DESeq2軟件包對(duì)3個(gè)表達(dá)矩陣進(jìn)行差異表達(dá)分析，從中提取出spike-in在A組和B組之間的差異表達(dá)信息。為了確認(rèn)在廣泛采用的差異基因篩選標(biāo)準(zhǔn)（P＜0.05） 下，總體RNA 法和lncRNA 法篩選差異lncRNA的效果，我們統(tǒng)計(jì)了spike-in在總體RNA表達(dá)矩陣和lncRNA 表達(dá)矩陣中假陽(yáng)性和假陰性的數(shù)目（表1），并計(jì)算了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的特異性更好。我們對(duì)總體RNA 表達(dá)矩陣中多出的5個(gè)假陽(yáng)性spike-in進(jìn)行了進(jìn)一步分析，除了其中一個(gè)spike-in 表達(dá)水平極低以外，其余4 個(gè)spikein的表達(dá)數(shù)據(jù)在經(jīng)過(guò)標(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個(gè)表達(dá)矩陣預(yù)測(cè)的差異表達(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 曲線。同時(shí)為了排除算法因素的影響，我們分別使用DESeq2 和edgeR 進(jìn)行ROC 分析（圖3），其中橫坐標(biāo)和縱坐標(biāo)分別用特異性和準(zhǔn)確性表示，它們?cè)跀?shù)值上分別與1-FPR 和TPR 相等。通過(guò)每條ROC 曲線下的面積大小，即AUC 值來(lái)量化不同分析方法篩選差異spike-in 的效果。在DESeq2 的分析結(jié)果中，Spike-in 在3 個(gè)表達(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個(gè)表達(dá)矩陣。

圖3 使用DESeq2和edgeR對(duì)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)各個(gè)樣本的測(cè)序數(shù)據(jù)來(lái)源完全相同，理論上可以認(rèn)為組內(nèi)樣本之間基因表達(dá)水平?jīng)]有差異，因此可以用來(lái)評(píng)估差異分析方法的FPR大小。本研究使用DESeq2 對(duì)總體RNA 表達(dá)矩陣和lncRNA 表達(dá)矩陣中的A 組和B 組分別進(jìn)行組內(nèi)差異lncRNA 分析，其分組為A1 和A2、A3 和A4、B1 和B2、B3 和B4。由于組內(nèi)樣本絕大多數(shù)基因無(wú)表達(dá)差異，P值無(wú)限接近于1，因此我們?cè)趯?duì)lncRNA 的P值作密度分布圖時(shí)，把范圍縮小在0～0.2之間。當(dāng)把篩選差異lncRNA的標(biāo)準(zhǔn)定為P＜0.05 時(shí)，在A 組中l(wèi)ncRNA 表達(dá)矩陣和總體RNA 表達(dá)矩陣的差異lncRNA 數(shù)目分別為9 個(gè)（占比2.57%）和7 個(gè)（占比2.45%），見(jiàn)圖4A；B 組中分別15 個(gè)（占比3.46%）和17 個(gè)（占比3.65%），見(jiàn)圖4B。使用不同表達(dá)矩陣篩選差異lncRNA 的FPR并沒(méi)有太大差別，可能與本研究組內(nèi)樣本之間差異基因數(shù)目過(guò)少有關(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 更高的原因，我們對(duì)mRNA 和lncRNA 的表達(dá)水平分布進(jìn)行了分析。鑒于不同基因之間的原始表達(dá)水平（raw counts）差距過(guò)大，我們對(duì)其作對(duì)數(shù)化處理后進(jìn)行展示?？梢钥吹皆谌繕颖局衜RNA 和lncRNA 各自表達(dá)水平的分布基本一致（圖5A），而mRNA和lncRNA之間則明顯不同（圖5B）。

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

3 討論

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

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

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

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

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

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