無生物學(xué)重復(fù)RNA-seq分析
CORNAS: coverage-dependent RNA-Seq analysis of gene expression data without biological replicates
BMC Bioinformatics 的一篇文章中提出了一種新的差異基因分析方法欢摄。
這篇文章提出了CORNAS(COverage-dependent RNA-Seq) 方法瞬沦，利用貝葉斯方法來推斷真實(shí)基因表達(dá)數(shù)的 后驗(yàn)分布。
其創(chuàng)新型之一該方法包括了由RNA樣品濃度決定的覆蓋度參數(shù)坏平，之二是真實(shí)基因表達(dá)量后驗(yàn)分布的比較為尋找差異表達(dá)基因提供了一個(gè)參照揍拆。
這種方法針對(duì)無重復(fù)樣本的數(shù)據(jù)是有一定優(yōu)勢(shì)的
目前現(xiàn)有的差異基因計(jì)算方法主要思路是**通過某種方法矯正后的觀測(cè)基因counts代表了實(shí)際為止的真實(shí)基因counts渠概。
**常用的矯正思路是首先考慮所含基因個(gè)數(shù)，基因表達(dá)量和樣本library大小等因素進(jìn)行樣本間的標(biāo)準(zhǔn)化嫂拴，然后使用校正后的基因counts數(shù)進(jìn)行有參或者無參差異表達(dá)檢驗(yàn)播揪。
這些軟件都有一個(gè)假設(shè)是觀測(cè)到的基因counts數(shù)充分代表了真實(shí)的基因表達(dá)情，但是這一假設(shè)并不是完全沒有問題筒狠。
在基因組測(cè)序中猪狈，覆蓋深度是一個(gè)非常重要的概念，但是在轉(zhuǎn)錄組測(cè)序中因?yàn)檗D(zhuǎn)錄組大小在不同組織或者不同細(xì)胞中都有差別辩恼，而且細(xì)胞間的mRNA種類也是高度變化的雇庙，因此覆蓋深度的概念并沒有很好的延伸到轉(zhuǎn)錄組中。依據(jù)2009年的一篇報(bào)道灶伊，想要準(zhǔn)確覆蓋人細(xì)胞系的95%轉(zhuǎn)錄本大致需要700M reads (但是這個(gè)說法個(gè)人認(rèn)為也不能全信)疆前，因?yàn)閷?shí)際測(cè)序深度遠(yuǎn)遠(yuǎn)比這個(gè)要少，所以測(cè)序帶來的強(qiáng)大的隨機(jī)效應(yīng)就對(duì)統(tǒng)計(jì)結(jié)果有非常大的影響聘萨，目前的通用做法是盡量增加樣本的重復(fù)竹椒。除此之外，現(xiàn)在已知 RNAseq 文庫制備和測(cè)序都存在一些誤差和偏好性匈挖。
從以上問題和角度出發(fā)碾牌，這篇文章提出了CORNAS(COverage-dependent RNA-Seq) 方法，利用貝葉斯方法來推斷真實(shí)基因表達(dá)數(shù)的后驗(yàn)分布儡循。
其創(chuàng)新型之一該方法包括了由RNA樣品濃度決定的覆蓋度參數(shù)，之二是真實(shí)基因表達(dá)量后驗(yàn)分布的比較為尋找差異表達(dá)基因提供了一個(gè)參照征冷。
這種方法針對(duì)無重復(fù)樣本的數(shù)據(jù)是有一定優(yōu)勢(shì)的择膝，在文章中，作者也是使用的真實(shí)無重復(fù)數(shù)據(jù)進(jìn)行了測(cè)試检激，另外比較對(duì)象也并非是DEseq和EdgeR這類有參檢驗(yàn)的軟件肴捉，而是NOISeq和GFOLD兩款軟件，至于使用這兩款軟件的原因作者提到是因?yàn)楹统Ｓ玫腄Eseq以及EdgeR相比叔收，這兩款軟件返回的假陽性數(shù)目更少齿穗。
關(guān)于如何計(jì)算樣本覆蓋度的問題可以參考原文，另外作者發(fā)現(xiàn)在給定真實(shí)基因表達(dá)數(shù)時(shí)描述觀測(cè)基因表達(dá)數(shù)最合適的模型是廣義泊松分布（而不是過去常用的負(fù)二項(xiàng)分布）饺律，將真實(shí)基因數(shù)和廣義泊松分布參數(shù)以及通過RNAseq樣品濃度得到測(cè)序覆蓋度相關(guān)聯(lián)就可以確定真實(shí)基因表達(dá)量的后驗(yàn)分布窃页，進(jìn)而用這個(gè)分布作為無重復(fù)RNA-seq試驗(yàn)進(jìn)行差異基因分析的基礎(chǔ)。
綜上，如果手上恰好有無重復(fù)試驗(yàn)的RNA-seq數(shù)據(jù)要做差異表達(dá)分析脖卖，不妨試試NOISeq乒省，GFOLD和這個(gè)新出的CORNAS 這幾個(gè)軟件。
摘要
背景:
目前在RNA-Seq實(shí)驗(yàn)中畦木，在調(diào)用差異表達(dá)基因的統(tǒng)計(jì)方法中袖扛，假設(shè)一個(gè)被調(diào)整的觀察到的基因計(jì)數(shù)代表一個(gè)未知的真實(shí)基因計(jì)數(shù)。
這種調(diào)整通常包括一個(gè)歸一化步驟來解釋異質(zhì)樣本庫的大小十籍，然后將結(jié)果歸一化的基因計(jì)數(shù)作為參數(shù)或非參數(shù)差異基因表達(dá)測(cè)試的輸入蛆封。
一個(gè)真實(shí)基因的分布，每一個(gè)都有不同的概率勾栗，可以導(dǎo)致相同的觀察到的基因計(jì)數(shù)娶吞。
重要的是，序列覆蓋信息目前沒有明確納入用于RNA-Seq分析的任何統(tǒng)計(jì)模型中械姻。
結(jié)果:
我們開發(fā)了一種快速貝葉斯方法妒蛇，該方法利用RNA樣本濃度確定的測(cè)序覆蓋信息來估計(jì)真實(shí)基因計(jì)數(shù)的后驗(yàn)分布。
與NOISeq和GFOLD相比楷拳，我們的方法有更好的或比較好的性能绣夺，根據(jù)模擬實(shí)驗(yàn)和真實(shí)的未復(fù)制數(shù)據(jù)的實(shí)驗(yàn)結(jié)果。我們將先前未使用的測(cè)序覆蓋參數(shù)納入RNA-Seq數(shù)據(jù)的差異基因表達(dá)分析的過程中欢揖。
結(jié)論:我們的研究結(jié)果表明陶耍，我們的方法可用于在有限數(shù)量的復(fù)制和低測(cè)序覆蓋率的實(shí)驗(yàn)中克服分析瓶頸。
**背景 **
與特定表型類型顯著相關(guān)的大規(guī)模的基因簽名挖掘是轉(zhuǎn)錄組分析通常期望的結(jié)果她混。
RNA測(cè)序（RNA-Seq）已成為基因表達(dá)譜分析的首選工具烈钞，它在幾個(gè)重要方面補(bǔ)充了傳統(tǒng)的微陣列：
它更徹底地采樣轉(zhuǎn)錄組，檢測(cè)異構(gòu)體坤按，并且在沒有對(duì)靶轉(zhuǎn)錄組的預(yù)先知識(shí)的情況下工作[1,2]毯欣。
自從發(fā)表第一篇RNA-Seq論文[3]以來，對(duì)RNA-Seq的廣泛興趣導(dǎo)致了測(cè)序平臺(tái)如454臭脓，Illumina和Solexa的快速開發(fā)和部署酗钞。這些平臺(tái)自然促成數(shù)據(jù)處理和分析方法的并行開發(fā)，以從RNA-Seq數(shù)據(jù)中提取生物學(xué)意義来累。
典型的RNA-Seq數(shù)據(jù)分析首先選擇通過質(zhì)量控制標(biāo)準(zhǔn)的讀數(shù)砚作，將它們映射到參考基因組，然后量化基因計(jì)數(shù)嘹锁。標(biāo)準(zhǔn)化后葫录，所得到的數(shù)據(jù)矩陣已準(zhǔn)備好進(jìn)行統(tǒng)計(jì)分析，對(duì)此可以使用令人眼花繚亂的數(shù)字替代方法[4,5]领猾。不管這些參數(shù)是參數(shù)（例如DESeq [6]米同，EdgeR [7]骇扇，DEGSeq [8]，BaySeq [9]）還是非參數(shù)（例如NOISeq [10]窍霞，SAMSeq [11]）匠题，觀察到的基因計(jì)數(shù)是實(shí)際真實(shí)基因計(jì)數(shù)的充分表示。
在基因組測(cè)序中但金，總閱讀長度與基因組大小的比率提供了覆蓋度量韭山，這對(duì)于評(píng)估組裝基因組的完整性很重要。然而冷溃，擴(kuò)展轉(zhuǎn)錄組大小的覆蓋概念并不簡(jiǎn)單钱磅。首先，同一生物體內(nèi)不同組織類型的轉(zhuǎn)錄組大小不同似枕，甚至相同組織類型的細(xì)胞之間轉(zhuǎn)錄組大小也不同[12]盖淡。接下來，細(xì)胞之間mRNA相對(duì)比例可以變化很大[13]凿歼。例如褪迟，在基因相同的酵母細(xì)胞中，已經(jīng)觀察到每個(gè)細(xì)胞超過800個(gè)mRNA拷貝的變異[14]答憔。
為了精確定量人細(xì)胞系中95％的轉(zhuǎn)錄物味赃，需要高達(dá)7億（M）的讀數(shù)[15]。相比之下虐拓，RNA-Seq實(shí)驗(yàn)通常產(chǎn)生的讀數(shù)遠(yuǎn)小于700M [16]心俗，足夠低以使隨機(jī)效應(yīng)對(duì)統(tǒng)計(jì)分析結(jié)果的解釋產(chǎn)生重大影響。如果不大幅度降低測(cè)序成本蓉驹，研究人員往往不得不優(yōu)先重復(fù)讀數(shù)總數(shù)的增加城榛，因?yàn)檫@是提高差異基因表達(dá)分析統(tǒng)計(jì)效能的最佳策略[17]。然后出現(xiàn)幾個(gè)挑戰(zhàn)态兴。假設(shè)有完美的閱讀圖譜和定量分析募判，目前還不清楚觀察到的基因計(jì)數(shù)是否代表真正的基因計(jì)數(shù)弧满，因?yàn)榇罅康恼鎸?shí)基因計(jì)數(shù)可能由于強(qiáng)大的隨機(jī)效應(yīng)而產(chǎn)生特定的觀察基因計(jì)數(shù)荸实。使問題復(fù)雜化是技術(shù)性RNA-Seq文庫制備和測(cè)序中固有的偏見[18]蜗帜，最近才引起嚴(yán)重關(guān)注的問題[19]。
我們開發(fā)了CORNAS（依賴于平衡的RNA-Seq）敢订，貝葉斯方法來推斷真實(shí)基因計(jì)數(shù)的后驗(yàn)分布。這種方法的新穎之處在于它結(jié)合了由RNA樣品濃度確定的覆蓋參數(shù)罢吃。隨后楚午，真實(shí)基因計(jì)數(shù)后驗(yàn)分布的比較為調(diào)用差異表達(dá)基因（DEG）提供了基礎(chǔ)。我們報(bào)告了CORNAS在未復(fù)制的RNA-Seq實(shí)驗(yàn)中的應(yīng)用尿招，并討論了其用于克服這些實(shí)驗(yàn)的分析局限性的可能性矾柜。
結(jié)果
真實(shí)基因計(jì)數(shù)和樣本覆蓋率的定義
我們首先將真基因計(jì)數(shù)定義為為測(cè)序運(yùn)行準(zhǔn)備的樣品中基因的mRNA拷貝總數(shù)阱驾。該定義適用于包含單個(gè)或多個(gè)單元格的樣本。由于后者原則上可以來源于多種不同的真實(shí)基因計(jì)數(shù)怪蔑，因此不能僅憑觀察到的基因計(jì)數(shù)確定該值里覆。然而，關(guān)于樣本覆蓋率的信息可以改進(jìn)估計(jì)真實(shí)基因計(jì)數(shù)的過程缆瓣。
樣品（b）的覆蓋度被定義為測(cè)序的cDNA片段（S）除以總cDNA片段群體大行稀（N）。單端測(cè)序產(chǎn)生一個(gè)讀數(shù)來代表一個(gè)cDNA測(cè)序弓坞，而配對(duì)末端測(cè)序產(chǎn)生兩個(gè)讀數(shù)代表一個(gè)測(cè)序的cDNA隧甚。
在Illumina測(cè)序方案中，樣品覆蓋率的計(jì)算可以基于mRNA樣品濃度渡冻。我們推斷在PCR之前的步驟中產(chǎn)生的cDNA的量為合理估計(jì)樣本覆蓋率提供了關(guān)鍵因素戚扳，因?yàn)椋?/p>
1）樣品文庫制備過程中的片段化步驟導(dǎo)致cDNA分子大小（500bp）的同質(zhì)性;
2）PCR后的體積和濃度已知（40μL的200nM cDNA）和;
3）PCR循環(huán)的數(shù)量是已知的（14個(gè)循環(huán)）族吻。cDNA片段進(jìn)行PCR以提高獲得至少一個(gè)測(cè)序覆蓋率的機(jī)會(huì)帽借。
假設(shè)完美的擴(kuò)增效率，每個(gè)cDNA片段在PCR過程中擴(kuò)增214次超歌。因此砍艾，PCR之前的cDNA片段數(shù)估計(jì)為4.818×1012 /214≈300M（詳見附加文件1）。我們使用這個(gè)數(shù)量作為估計(jì)的片段總體大小來確定覆蓋范圍握础，因?yàn)樗c我們開始使用的mRNA數(shù)量非常相似辐董。
機(jī)會(huì)機(jī)制為觀察基因計(jì)數(shù)產(chǎn)生廣義泊松分布
當(dāng)將cDNA片段加載到測(cè)序運(yùn)行中時(shí)，假定從加載的cDNA片段隨機(jī)產(chǎn)生短讀數(shù)禀综。因此简烘，真正的基因計(jì)數(shù)會(huì)誘導(dǎo)觀察到的基因計(jì)數(shù)的概率分布。為了找到最能描述后者的概率模型定枷，
我們進(jìn)行了一系列模擬來確定均方差觀察到的基因計(jì)數(shù)在6個(gè)覆蓋值下的關(guān)系：0.5,0.4,0.25,0.1,0.01和0.001孤澎。假設(shè)分別為150M，120M欠窒，75M覆旭，30M，3M和0.3M讀數(shù)計(jì)算這些覆蓋率從總片段群體大小300M開始測(cè)序岖妄。對(duì)于每個(gè)覆蓋范圍型将，我們生成了1到100,000的實(shí)際計(jì)數(shù)值的觀察計(jì)數(shù)的經(jīng)驗(yàn)分布（詳見“方法”部分）。
模擬結(jié)果提供了三個(gè)重要觀察結(jié)果：觀察計(jì)數(shù)的平均值與覆蓋率成正比荐虐，隨著覆蓋率的增加七兜，發(fā)生不均勻分布（即方差小于平均值）（附加文件1：圖S1），線性模型充分描述了平均方差比和覆蓋率之間的關(guān)系（等式6）福扬。這些結(jié)果表明廣義泊松（GP）分布[20]適用于給定真實(shí)基因計(jì)數(shù)（T）的觀測(cè)基因計(jì)數(shù)（X）的分布腕铸。帶GP的GP的概率質(zhì)量函數(shù)參數(shù)λ1和λ2由下式給出
意味著λ2= 1-√m惜犀，其中方差比（0 <m <4）。觀察到的基因數(shù)的平均值,如果真實(shí)基因數(shù)與覆蓋率b和真實(shí)計(jì)數(shù)k的乘積成比例狠裹，則給出λ1=bk√m,具有平均值λ1的泊松分布是當(dāng)m = 1時(shí)GP的特例虽界。
用于估計(jì)真實(shí)基因計(jì)數(shù)的貝葉斯模型觀察到的基因數(shù)和測(cè)序覆蓋率
公式中GP模型的重要性。 1源于反向條件使我們能夠考慮的事實(shí)給出觀察基因計(jì)數(shù)和測(cè)序覆蓋率的真實(shí)基因計(jì)數(shù)（T）的概率分布（即真基因計(jì)數(shù)的后驗(yàn)分布）涛菠。讓我們假設(shè)一個(gè)統(tǒng)一的先驗(yàn)分布的T值超過1莉御，2，...碗暗。 颈将。貝葉斯定理的應(yīng)用產(chǎn)生：
其中k≥x。請(qǐng)注意言疗，盡管我們之前使用的是不正確的晴圾，但后驗(yàn)分布是正確的。有趣的是噪奄，伽馬分布提供了一個(gè)很好的近似值等式2（數(shù)學(xué)證明見[21]）死姚。在這里，我們發(fā)現(xiàn)如果伽馬分布的平均值μ和方差σ2與覆蓋度b和觀測(cè)基因數(shù)x相關(guān)勤篮，則近似值非常好（附加文件1：圖S2）：
因此都毒，近似伽瑪分布的概率密度由下式給出
近似提供計(jì)算有效的手段來計(jì)算真實(shí)基因計(jì)數(shù)的后驗(yàn)分布的累積分布函數(shù)。在未復(fù)制的RNA-Seq實(shí)驗(yàn)中調(diào)用差異表達(dá)基因的統(tǒng)計(jì)檢驗(yàn)可以基于真實(shí)基因計(jì)數(shù)的后驗(yàn)分布
（等式2）如下碰缔。對(duì)于單個(gè)對(duì)照和單個(gè)治療樣本账劲，如果我們有關(guān)于測(cè)序的信息
覆蓋對(duì)照樣品（b0）和處理樣品（b1），然后金抡，給定觀察到的基因計(jì)數(shù)
對(duì)照（x0）和治療（x1）組瀑焦，其真實(shí)基因計(jì)數(shù)的后驗(yàn)分布近似
伽馬（公式5）。如果后者的后驗(yàn)平均值較大梗肝，我們宣布一個(gè)基因在治療組中差異上調(diào)榛瓮，其第0.5個(gè)百分點(diǎn)至少比對(duì)照組的第99.5個(gè)百分點(diǎn)高1.5倍（默認(rèn)）。相反巫击，如果后者具有較小的后部禀晓，則基因在治療組中被差異下調(diào),平均而言，其第99.5百分位至少是對(duì)照組的第0.5百分位的1.5倍（默認(rèn)）（圖坝锰。1）粹懒。這個(gè)過程很快，因?yàn)橘が敺植嫉陌俜治粩?shù)很容易計(jì)算出來顷级。此外崎淳，使用這個(gè)程序聲明基因差異表達(dá)意味著兩個(gè)樣本中真實(shí)基因計(jì)數(shù)的差異為0.9952≈0.99至少1.5折。
CORNAS****的性能評(píng)估
我們使用模擬和實(shí)際數(shù)據(jù)集進(jìn)行了一系列比較CORNAS與NOISeq [10]和GFOLD [22]性能的測(cè)試愕把。我們選擇了GFOLD和NOISeq拣凹，因?yàn)閮烧叨加袌?bào)道在應(yīng)用于未復(fù)制的RNA-Seq數(shù)據(jù)集時(shí)，與其他流行方法（如DESeq2和edgeR [5]）相比恨豁，在標(biāo)記為差異表達(dá)的基因中返回相對(duì)較少的假陽性數(shù)[5]嚣镜。
測(cè)試1：在模擬的真實(shí)基因計(jì)數(shù)數(shù)據(jù)中檢測(cè)差異表達(dá)的基因
我們使用四個(gè)覆蓋范圍測(cè)試了CORNAS：0.5,0.25,0.1和0.01，模擬真實(shí)基因計(jì)數(shù)范圍從1到10,000橘蜜。記錄調(diào)用差異表達(dá)基因（DEG）的相對(duì)頻率為100個(gè)獨(dú)立樣本
對(duì)照組與對(duì)照組之間無倍數(shù)變化（無效）菊匿，1.5倍變化（弱效應(yīng)）和2倍變化（更強(qiáng)效）的試驗(yàn)。在沒有倍數(shù)變化的情況下计福，假陽性率（FPR）被估計(jì)為DEG呼叫率跌捆。真實(shí)陽性率（TPR）或靈敏度是弱和強(qiáng)效應(yīng)情景下的DEG呼叫率。
一般來說象颖，我們觀察到隨著覆蓋率和覆蓋率的增加佩厚，誤報(bào)率降低和DEG呼叫率增加
越來越多的真實(shí)基因計(jì)數(shù)（圖2）。與GFOLD和CORNAS默認(rèn)值相比说订，當(dāng)真實(shí)基因計(jì)數(shù)較低時(shí)抄瓦，NOISeq產(chǎn)生最大的FPR。 NOISeq的靈敏度通常很好陶冷，除非0.01的低覆蓋率;
當(dāng)真實(shí)數(shù)量超過1,000時(shí)钙姊，其DEG呼叫率開始下降。 GFOLD表現(xiàn)出非常低的靈敏度埂伦，其中
與[5]中報(bào)道的保守行為是一致的煞额。CORNAS顯示出對(duì)FPR的極好控制，并且依賴于檢測(cè)的倍數(shù)變化閾值DEG在弱信號(hào)和強(qiáng)信號(hào)情況下沾谜。例如膊毁，CORNAS默認(rèn)值（φ= 1.5）在弱信號(hào)情況下表現(xiàn)非常差，因此如果檢測(cè)到這些基因是有意義的类早，那么應(yīng)將φ調(diào)整為較低的值媚媒，如1（CORNAS set1）。一般而言涩僻，CORNAS的靈敏度隨著真實(shí)計(jì)數(shù)的增加而增加缭召，并且覆蓋值為0.1或更大時(shí)迅速收斂到1。
測(cè)試2：compcodeR模擬
觀察到的基因計(jì)數(shù)的分布通常使用負(fù)二項(xiàng)式分布建模逆日，并且compcodeR R軟件包[23]提供了一個(gè)模擬器來模擬基于這種分布的RNA-Seq計(jì)數(shù)數(shù)據(jù)嵌巷。假定基因長度相等并設(shè)定在1000個(gè)堿基處。我們使用[23]中提供的示例來創(chuàng)建對(duì)照組治療比較（每組5個(gè)重復(fù)）室抽，其中對(duì)照組中有624個(gè)上調(diào)基因和625個(gè)下調(diào)基因搪哪，用于模擬12,498個(gè)基因的轉(zhuǎn)錄組。從這個(gè)數(shù)據(jù)矩陣中坪圾，共構(gòu)建了25個(gè)未復(fù)制的數(shù)據(jù)集晓折。對(duì)于CORNAS惑朦，我們?cè)u(píng)估了兩種不同覆蓋率的結(jié)果
對(duì)樣本進(jìn)行比較;估計(jì)比compcodeR覆蓋率（CORNAS_10xless）小10倍，另一個(gè)小于100倍（CORNAS_100xless）（Supporting Data）漓概。我們做了兩個(gè)單獨(dú)的NOISeq運(yùn)行漾月，一個(gè)沒有長度歸一化（NOISeq_nln），另一個(gè)使用M值歸一化的修整均值（NOISeq_tmmnl）胃珍。
所有方法的陽性預(yù)測(cè)值（PPV）和敏感度均較低;盡管如此梁肿，CORNAS顯示出比其他方法更高的靈敏度，而GFOLD具有相對(duì)更好的PPV（圖3a）觅彰。所有方法的F分?jǐn)?shù)非常相似（表1）吩蔑。與其他方法相比，CORNAS稱為更大的DEG集大小填抬。與NOISeq_nln不同烛芬，由CORNAS調(diào)用的更大的DEG集合大小并未顯著降低其PPV。CORNAS_100xless和CORNAS_10xless表現(xiàn)出類似的表現(xiàn)痴奏。
比較的平均運(yùn)行時(shí)間約為NOISeq_nln和NOISeq_tmmnl的三分鐘蛀骇，GFOLD為一分鐘，CORNAS_10xless和CORNAS_100xless為三秒读拆。
測(cè)試3：人類性別特異性基因表達(dá)
CORNAS對(duì)真實(shí)數(shù)據(jù)適用性的評(píng)估基于人類淋巴母細(xì)胞RNA-Seq來自Pickrell的研究數(shù)據(jù)[24]擅憔。在這個(gè)數(shù)據(jù)集中，男性和女性的性別構(gòu)成了兩個(gè)表型類別檐晕，所以真正的DEG可以純粹使用生物推理來確定暑诸。差異表達(dá)的基因被鑒定為具有Y染色體相關(guān)表達(dá)的19個(gè)基因[5]。在生物學(xué)基礎(chǔ)上沒有差異表達(dá)的基因包括61個(gè)X失活的(XiE)基因[25,26]和11個(gè)看家基因[27]辟灰。
我們從725對(duì)可能的組合中隨機(jī)選取100對(duì)單身女性 - 單身男性對(duì)（29名女性个榕，25名男性）（支持?jǐn)?shù)據(jù)），并比較GFOLD芥喇，NOISeq和CORNAS西采。我們的研究結(jié)果表明，NOISeq與CORNAS和GFOLD相比表現(xiàn)不佳继控，而GFOLD表現(xiàn)略好于CORNAS（圖3b械馆，表1）。但是武通，類似于compcodeR仿真結(jié)果霹崎，CORNAS稱為較大的DEG集。
NOISeq的平均運(yùn)行時(shí)間約為2分鐘冶忱，GFOLD為30秒尾菇，CORNAS為10秒。
測(cè)試4：在組織特異性基因表達(dá)數(shù)據(jù)中的覆蓋效應(yīng)
Marioni數(shù)據(jù)集[28]由來自人類肝臟和腎臟的RNA-Seq數(shù)據(jù)在兩種不同的序列中組成加載濃度，3 pM（高）和1.5 pM（低）派诬。我們調(diào)查了CORNAS是否會(huì)被簡(jiǎn)單地根據(jù)不同濃度進(jìn)行DEG調(diào)用時(shí)被誤導(dǎo)劳淆，當(dāng)兩個(gè)樣品都來自同一組織時(shí)。在CORNAS中假陽性率低千埃，在相同濃度的相同組織樣品中沒有DEG呼叫進(jìn)行比較（附加文件1：表S1）憔儿。但是，對(duì)于不同的樣品GFOLD顯示出比CORNAS更少的誤報(bào)放可。在所有情況下，NOISeq都返回了最高FPR朝刊。
一組4863個(gè)基因被鑒定為在組織表達(dá)數(shù)據(jù)庫TISSUES中編目的人肝臟或腎臟組織中唯一表達(dá)[29]耀里。同樣，與CORNAS和GFOLD相比拾氓，NOISeq表現(xiàn)不佳冯挎，而CORNAS表現(xiàn)最好（圖3c，表1）咙鞍。對(duì)于不同組織類型之間的所有12種比較房官，最大的DEG組被CORNAS調(diào)用，并且NOISeq最小的那個(gè)续滋。
通常對(duì)于不同的組織類型翰守，由NOISeq和GFOLD調(diào)用的DEG集顯示差的重疊，與GFOLD和CORNAS之間的重疊相比疲酌，和NOISeq和CORNAS之間（附加文件1：圖S3）蜡峰。CORNAS表示針對(duì)不同組織類型的更獨(dú)特的DEG呼叫。與此同時(shí)朗恳，來自GFOLD或NOISeq的大部分DEG呼叫也被CORNAS調(diào)用湿颅。
NOISeq的平均運(yùn)行時(shí)間約為5分鐘，GFOLD約為30秒粥诫，CORNAS為5秒油航。
PCR擴(kuò)增效率對(duì)靈敏度的影響
雖然我們?cè)跇?gòu)建模型時(shí)假設(shè)了完美的PCR擴(kuò)增效率，但我們?nèi)栽谠u(píng)估可能性影響95％怀浆，90％谊囚，85％和80％PCR效率對(duì)CORNAS靈敏度和FPR的影響（詳見“方法”一節(jié)）。正如我們沒有的那樣揉稚，CORNAS似乎對(duì)小的違反完美PCR擴(kuò)增效率的魯棒性,即使PCR效率達(dá)到80％秒啦，也會(huì)發(fā)現(xiàn)靈敏度和FPR的實(shí)質(zhì)性變化。曲線下面積（AUC）為所有四種測(cè)試預(yù)期覆蓋率的接收者操作特征（ROC）圖均小于5％（圖4和附加文件1：圖S4）搀玖。
討論
CORNAS作為估計(jì)真實(shí)基因數(shù)量的框架
作為RNA-Seq計(jì)數(shù)數(shù)據(jù)建模中負(fù)二項(xiàng)分布的替代方案余境，GP模型正在被越來越多地研究[30-33]。在這里，我們展示了一種機(jī)會(huì)機(jī)制芳来，它自然使GP成為觀察基因計(jì)數(shù)數(shù)據(jù)的模型含末。通過將GP的參數(shù)與真實(shí)基因計(jì)數(shù)和測(cè)序聯(lián)系起來使用RNA樣本濃度的覆蓋率，我們因此能夠確定真實(shí)基因計(jì)數(shù)的后驗(yàn)分布即舌。這種分布形成了在未復(fù)制的RNA-Seq實(shí)驗(yàn)中進(jìn)行DEG調(diào)用的基礎(chǔ)佣盒。
目前，在生物體的基因模型上映射的讀取深度用于估計(jì)RNA-Seq實(shí)驗(yàn)中的覆蓋度顽聂。我們知道樣品中mRNA的總量在Illumina測(cè)序儀中未被捕獲肥惭，其具有固定的有限飽和量，其可能超過或低于樣品濃度紊搪。覆蓋面通常被認(rèn)為是代表性不足蜜葱，這是一種限制這通常被認(rèn)為是可以通過深度測(cè)序糾正的，該測(cè)序用于檢測(cè)mRNA表達(dá)非常低的基因[15,17,34]耀石。我們覆蓋參數(shù)的范圍（0到1之間）應(yīng)該覆蓋最多沒有完成深度排序的實(shí)際案例牵囤。如果估計(jì)的覆蓋范圍不止一個(gè)，我們不建議使用CORNAS滞伟。
我們的研究做了幾個(gè)假設(shè)來簡(jiǎn)化模型揭鳞，其中之一是100％有效的PCR擴(kuò)增。我們模擬的PCR擴(kuò)增效率的影響確實(shí)表明梆奈，當(dāng)我們高估覆蓋率時(shí)野崇，靈敏度和FPR增加，但是這種變化沒有不利影響鉴裹。為了保持觀察計(jì)數(shù)模型舞骆，在當(dāng)前工作中假設(shè)理想的隨機(jī)cDNA片段采樣（因此后驗(yàn)分布）我們可以非常簡(jiǎn)單地研究將覆蓋參數(shù)引入到DEG呼叫過程中的效果。由于真正的RNA-Seq實(shí)驗(yàn)含有文庫制備偏差径荔，通過完整的測(cè)序過程模擬器督禽，如rlsim [35]，可以更好地探索這種偏差的影響总处。
在我們的模擬中未明確處理觀察到的基因計(jì)數(shù)的潛在變異來源,涉及不同算法將短讀數(shù)映射到參考基因組的方式（例如使用BWA [36]狈惫，OSA [37]，TopHat [38]和Bowtie [39]）鹦马，以及如何量化這些映射讀數(shù)（例如使用HTSeq [40]和Cufflinks [38]）胧谈。我們建議由于這種變異來源觀察到的基因計(jì)數(shù)的變化是相對(duì)的不重要，因此不會(huì)嚴(yán)重影響真實(shí)基因計(jì)數(shù)的后驗(yàn)分布荸频。首先菱肖，可以提高讀取對(duì)齊質(zhì)量的算法[41]，從而最大限度地減少計(jì)數(shù)錯(cuò)誤旭从。此外稳强，read-mapper和基因計(jì)數(shù)量化的組合已經(jīng)經(jīng)過實(shí)驗(yàn)研究场仲，并且可獲得最可靠的最佳建議以獲得最可靠的觀察到的基因計(jì)數(shù)（如[42]所建議的OSA + HTSeq）。
CORNAS的穩(wěn)健性
盡管如此退疫，CORNAS在compcodeR仿真中表現(xiàn)出與GFOLD和NOISeq相當(dāng)?shù)男阅?基于觀察到的基因計(jì)數(shù)的不同數(shù)據(jù)模型（即廣義泊松比負(fù)二項(xiàng)式）渠缕。該發(fā)現(xiàn)提供了將CORNAS框架整合到當(dāng)前RNA-Seq數(shù)據(jù)分析方案中的信心。此外褒繁，盡管覆蓋范圍已被估計(jì)亦鳞，并因此出現(xiàn)錯(cuò)誤，但兩個(gè)CORNAS設(shè)置（10xless和100xless）平均顯示類似的性能棒坏。CORNAS達(dá)成了一個(gè)很好的妥協(xié),與GFOLD和NOISeq相比燕差，靈敏度，PPV和DEG設(shè)置大小之間俊抵。在現(xiàn)實(shí)世界的實(shí)驗(yàn)中谁不，CORNAS如果能夠更可靠地確定覆蓋范圍，則可以超越競(jìng)爭(zhēng)方法徽诲，例如來自Marioni數(shù)據(jù)集在測(cè)試4中。
如果不包含coverage參數(shù)中的信息吵血，諸如GFOLD和NOISeq等傳統(tǒng)方法用于分析未復(fù)制的RNA-Seq計(jì)數(shù)數(shù)據(jù)要么過于保守谎替，要么很少打電話，但其中大多數(shù)是真正的陽性（GFOLD）蹋辅，或者在非常低的覆蓋范圍情況下（例如b = 0.01）做出相對(duì)更多的誤報(bào)（NOISeq）（圖2）钱贯。另一方面，我們發(fā)現(xiàn)CORNAS很好地控制了FPR并且當(dāng)覆蓋率不太小時(shí)（例如b≥0.1）具有高TPR侦另。此外秩命，如果檢測(cè)到較弱的倍數(shù)變化差異是令人感興趣的，
那么可以將折變參數(shù)（φ）從1.5減小到1.0（詳見“方法”部分）褒傅。對(duì)于弱信號(hào)和強(qiáng)信號(hào)弃锐，CORNAS在倍數(shù)變化參數(shù)為1.0時(shí)的TPR分布與NOISeq類似，除非覆蓋率非常低殿托。隨著真實(shí)基因數(shù)量的增加霹菊，CORNAS繼續(xù)顯示TPR普遍上升，而NOISeq則下降支竹。
現(xiàn)在旋廷，大多數(shù)RNA-Seq實(shí)驗(yàn)在測(cè)序前沒有報(bào)告起始材料中實(shí)際RNA量的估計(jì)值。因此礼搁，我們只能通過模擬研究使用后驗(yàn)均值校正觀察基因數(shù)的效果饶碘。鑒于令人鼓舞的結(jié)果，研究人員可能希望收集關(guān)于未來覆蓋參數(shù)的信息馒吴，以利用CORNAS分析真實(shí)RNA-Seq數(shù)據(jù)集扎运。
分析未復(fù)制的RNA-Seq計(jì)數(shù)數(shù)據(jù)時(shí)的一個(gè)主要問題是在缺少生物學(xué)重復(fù)時(shí)缺乏有效的歸一化方法瑟曲。在這里，我們已經(jīng)展示了CORNAS所基于的貝葉斯框架绪囱，通過處理基因真實(shí)計(jì)數(shù)的后驗(yàn)分布來避免歸一化問題测蹲。結(jié)果，轉(zhuǎn)錄物長度信息不是必需的鬼吵。這使得CORNAS適用于具有不完整或不斷演變的轉(zhuǎn)錄組參考數(shù)據(jù)的生物體扣甲，因?yàn)樾碌霓D(zhuǎn)錄信息不會(huì)改變真實(shí)計(jì)數(shù)隨著時(shí)間的推移如何估計(jì)。我們的研究結(jié)果表明齿椅，CORNAS可以用作克服分析瓶頸的手段琉挖，有限的重復(fù)和低測(cè)序覆蓋率導(dǎo)致DEGs，使用定量PCR和NanoString nCounter [43]等平臺(tái)進(jìn)行下游驗(yàn)證具有更好的前景涣脚。將CORNAS擴(kuò)展到多重復(fù)制的結(jié)果將在其他地方發(fā)布示辈。
結(jié)論
我們開發(fā)了CORNAS（COverage-dependent RNA-Seq），一種快速貝葉斯方法遣蚀，它結(jié)合了一個(gè)新的覆蓋參數(shù)來估計(jì)真實(shí)基因計(jì)數(shù)的后驗(yàn)分布矾麻。在CORNAS框架下，可以使用由RNA樣品濃度確定的來自序列覆蓋的正交信息來提高調(diào)用DEG的準(zhǔn)確性芭梯。通過對(duì)實(shí)際數(shù)據(jù)集的模擬和分析险耀，我們發(fā)現(xiàn)在未復(fù)制的RNA-Seq實(shí)驗(yàn)的情況下，CORNAS的性能與GFOLD和NOIseq相當(dāng)或更好玖喘。
方法
CORNAS作為R程序?qū)嵤┧ξ晒┫螺d（https://github.com/joel-lzb/CORNAS）。用于模擬和數(shù)據(jù)分析工作的Perl和R腳本可在https://github.com/joel-lzb/CORNAS_上找到Supporting_Data累奈。我們?cè)谶\(yùn)行RedHat 6操作系統(tǒng)的96 GB RAM的IBM System x3650 M3（2x6核心Xeon 5600）機(jī)器上進(jìn)行了計(jì)算機(jī)實(shí)驗(yàn)贬派。圖形使用R [44]中的ggplot繪制。模擬片段抽樣過程以及覆蓋率和觀測(cè)計(jì)數(shù)的均值/方差之間的關(guān)系,考慮具有相同長度的N個(gè)cDNA片段的群體澎媒。在這項(xiàng)研究中搞乏，我們?cè)O(shè)置N = 300×106（300M）。我們使用以下數(shù)目的測(cè)序讀數(shù)（S）：150M旱幼，120M查描，75M，30M柏卤，3M和0.3M冬三，用于模擬0.5,0.4,0.25,0.1,0.01和0.001的覆蓋率（S / N）分別。為了模擬來自cDNA片段群體的取樣過程缘缚，我們首先對(duì)每個(gè)N cDNA進(jìn)行索引分子從1到N.接下來勾笆，我們使用Fisher-Yates shuffle算法來洗牌指數(shù)，創(chuàng)建置換π=（π1桥滨，π2窝爪，...弛车，πN）。π的前S個(gè)元素表示測(cè)序片段的索引蒲每。對(duì)于從1到100,000的每個(gè)真實(shí)基因計(jì)數(shù)k纷跛，我們確定相應(yīng)的觀察基因數(shù)為</pre>
其中IA是事件A為真時(shí)的值為1的指示器函數(shù)，否則為0邀杏∑兜欤總共進(jìn)行2000次迭代，并從中估計(jì)觀察到的基因計(jì)數(shù)的平均值和方差望蜡。理論上唤崭，從這個(gè)過程中產(chǎn)生的觀察計(jì)數(shù)遵循超幾何分布。因此我們能夠計(jì)算給定覆蓋率b的超幾何分布的均值和方差（m）的比率為：
如果N非常大（300M）脖律，k遠(yuǎn)小于N（≤100K）谢肾，b = S / N。對(duì)于足夠小的b小泉，m≈1+ b芦疏。
抽樣過程自然導(dǎo)致觀測(cè)計(jì)數(shù)的超幾何分布，因?yàn)镹是有限的微姊。然而眯分，實(shí)際上N很大并且
未知，因此需要一個(gè)沒有上限的近似分布（參見GP模型中的“結(jié)果”部分）柒桑。
將后驗(yàn)均值和后驗(yàn)變化建模為覆蓋的函數(shù)
觀察計(jì)數(shù)x的真實(shí)計(jì)數(shù)k的后驗(yàn)分布如公式2。然后我們確定了后驗(yàn)分布的均值和方差與覆蓋參數(shù)的關(guān)系噪舀。我們首先將均值和方差建模為線性函數(shù)魁淳，使得：
其中參數(shù)Gm和Gs是梯度，并且Im和Is分別是截距与倡。然后界逛，我們?yōu)槊總€(gè)參數(shù)擬合模型
在各種覆蓋范圍進(jìn)行模擬（附加文件1：圖S2）。等式3和4是模擬均值的最終近似值和作為觀察基因計(jì)數(shù)（x）和測(cè)序范圍（b）的函數(shù)的后驗(yàn)分布的變化纺座。
CORNAS評(píng)估
程序設(shè)置
需要在CORNAS中設(shè)置兩個(gè)參數(shù)息拜。第一個(gè)是α，其用于確定較低（1-α）/ 2×100百分位數(shù)（p（1-α）/ 2）和上（1 +α）/ 2×100百分位數(shù)（p（1 +α）/ 2）净响。第二個(gè)參數(shù)是折變截止φ少欺。為了進(jìn)行DEG調(diào)用，我們需要p +（1-α）/ 2 / p-（1 +α）/2≥φ馋贤，其中上標(biāo)+和 - 分別表示后驗(yàn)分布具有較高和較低的平均值赞别。默認(rèn)設(shè)置是α= 0.99和φ= 1.5。這些值可以改變CORNAS更保守（例如增加α和/或φ）或更自由（例如降低α和/或φ）配乓。NOISeq在q = 0.9截止時(shí)運(yùn)行仿滔。 GFOLD運(yùn)行的折疊變化為0.01惠毁。如果基因的表達(dá)被認(rèn)為是上調(diào)的GFOLD值為1或更大并且如果GFOLD值為-1或更小則下調(diào)
性能指標(biāo)
真陽性（TP）是已知在兩個(gè)樣品之間差異表達(dá)的基因，并且通過評(píng)估的方法檢測(cè)為DEG崎页。假DEG電話是誤報(bào)（FP）鞠绰，而假陰性（FN）則錯(cuò)過了真正的DEG呼叫。對(duì)于DEG調(diào)用方法飒焦，其正向預(yù)測(cè)值（PPV）是正確調(diào)用的比例DEG（TP /（TP + FP））蜈膨，其靈敏度是稱為真實(shí)DEG的比例（TP /（TP + FN））。我們共同考慮了測(cè)試2,3和4的每種方法的靈敏度和PPV荒给。計(jì)算每個(gè)比較的F評(píng)分（靈敏度和PPV的調(diào)和平均值）為2×（靈敏度×PPV）/（靈敏度+ PPV）??丈挟。報(bào)道了每種方法的Themean F評(píng)分。
對(duì)于測(cè)試1志电，假陽性率（FPR）根據(jù)無倍數(shù)變化情況確定為真陰性（TN）明確已知
（FP /（FP + TN））曙咽，而靈敏度的計(jì)算方法與弱，強(qiáng)效應(yīng)場(chǎng)景中的測(cè)試2,3和4相似挑辆。
測(cè)試1：在模擬的真實(shí)基因計(jì)數(shù)數(shù)據(jù)中檢測(cè)差異表達(dá)的基因
對(duì)于此模擬例朱，考慮了三種生物效應(yīng)情景：無倍數(shù)變化（無效），1.5倍變化（弱效果）鱼蝉，2倍變化（強(qiáng)效果）洒嗤。在這三種情況下考慮的最大真實(shí)數(shù)分別為10,000,6,666和5,000。
我們假定每個(gè)真實(shí)的基因計(jì)數(shù)都是由一個(gè)基因發(fā)出的魁亦，因此在所有三種情況下所有真實(shí)基因計(jì)數(shù)的集合總共對(duì)應(yīng)于21,666個(gè)基因渔隶。按照模擬片段抽樣過程中描述的程序產(chǎn)生每個(gè)基因的觀察計(jì)數(shù)〗嗄危總共進(jìn)行了100次迭代說明觀察到的基因計(jì)數(shù)的取樣變異性间唉。在特定方法需要基因長度信息的情況下，我們將其設(shè)置為1000個(gè)堿基利术。
測(cè)試2：compcodeR模擬
我們根據(jù)compcodeR論文[23]中描述的方法生成模擬數(shù)據(jù)集B_625_625呈野。DEG的數(shù)量占基因總數(shù)的10％（12,498）。
測(cè)試3：人類性別特異性基因表達(dá)
對(duì)于來自尼日利亞的29名女性和25名男性組成的Pickrell研究印叁，我們使用了總數(shù)
來自已發(fā)表論文[24]的測(cè)序讀數(shù)和來自ReCount數(shù)據(jù)庫的RNA-Seq計(jì)數(shù)數(shù)據(jù)[45]被冒。每個(gè)樣品的測(cè)序覆蓋率計(jì)算為報(bào)告的總讀數(shù)除以標(biāo)準(zhǔn)300M cDNA片段大小。對(duì)于具有多個(gè)測(cè)序運(yùn)行的樣品轮蜕，我們?nèi)∷傻目傋x數(shù)的平均值昨悼。
測(cè)試4：在組織特異性基因表達(dá)數(shù)據(jù)中的覆蓋效應(yīng)
在Marioni數(shù)據(jù)集中，相同的人類肝臟和腎臟樣品在每條七條泳道中測(cè)序肠虽，其中5個(gè)泳道的RNA濃度為3pM幔戏，另外兩個(gè)泳道的濃度為1.5pM。14條泳道分兩次進(jìn)行測(cè)序税课。為了減少技術(shù)變化闲延，我們僅使用來自運(yùn)行2的數(shù)據(jù)痊剖，其中不同濃度的加載在相同條件下運(yùn)行和時(shí)間。我們將代表樣品轉(zhuǎn)錄組的cDNA片段的數(shù)量估計(jì)為負(fù)載濃度的乘積垒玲，加載量（假定標(biāo)準(zhǔn)為120μL）和Avogadro常數(shù)6.022×1023mol-1陆馁。根據(jù)2016年6月14日提取自TISSUES [29]的策劃信息確定了所使用的真實(shí)的DEGs集。我們選擇了737個(gè)人類腎臟基因和4,126個(gè)人類肝臟基因合愈，這些基因支持實(shí)驗(yàn)驗(yàn)證結(jié)果,并且可以用Ensembl基因ID鑒定叮贩。
PCR****擴(kuò)增效率對(duì)靈敏度的影響
使用與測(cè)試1相同的數(shù)據(jù)集進(jìn)行評(píng)估。為了模擬研究中PCR擴(kuò)增效率的影響佛析，
我們重新計(jì)算了每種COSNAS的測(cè)序覆蓋率益老，通過減少由PCR不同擴(kuò)增效率引起的PCR之前假定的片段總數(shù)（分別為95％，90％寸莫，85％捺萌，80％擴(kuò)增效率的總片段分別為70％，49％膘茎，34％桃纯，23％）。有關(guān)如何確定PCR擴(kuò)增效率的詳細(xì)信息可以在附加文件1中找到披坏。
例如态坦，有一個(gè)樣本完全擴(kuò)增，但測(cè)序范圍為0.25棒拂，在PCR之前將有300M片段伞梯，并產(chǎn)生75M讀數(shù)。假設(shè)產(chǎn)生的讀數(shù)保持不變帚屉，但PCR擴(kuò)增效率現(xiàn)在為95％壮锻，則測(cè)序覆蓋率估計(jì)將為0.36（75M /（300M×0.7））。新的覆蓋范圍然后用于0.25覆蓋的CORNAS以95％的PCR擴(kuò)增效率運(yùn)行涮阔。對(duì)于每個(gè)覆蓋范圍，根據(jù)在無效情況下調(diào)用的DEG的數(shù)量計(jì)算FPR灰殴，
并且從強(qiáng)效果方案中調(diào)用DEG計(jì)算靈敏度敬特。我們使用ROCR R軟件包[46]生成了接收器操作特性（ROC）曲線。通過固定α= 0.99然后從1.5變化到0.75獲得用于制作差異表達(dá)呼叫的截止點(diǎn)牺陶，并通過固定φ= 0.75伟阔，然后將α從0.99變?yōu)?.01。
