最近一直在做lncRNA的分析燕侠，其中的lncRNA的差異表達(dá)分析中镇草，需要對(duì)reads count 進(jìn)行歸一化，之前沒有考慮很多恍飘，就用的通常的流程：

hisat2→stringtie→prepDE.py/featureCount→DESeq2

其中的DESeq2 的歸一化部分榨崩，也是我們通常稱的標(biāo)準(zhǔn)化，是我們關(guān)注的重點(diǎn)常侣，DESeq2主要原理：通過計(jì)算一個(gè)歸一化因子蜡饵，并進(jìn)行變換，進(jìn)而提高中等表達(dá)基因的地位胳施。
歸一化和標(biāo)準(zhǔn)化總是搞不清溯祸，我還特地查了二者的共性和區(qū)別：

共性：歸一化和標(biāo)準(zhǔn)化本質(zhì)上都是一種線性變換。線性變換保持線性組合與線性關(guān)系式不變，這保證了特定模型不會(huì)失效焦辅，歸一化和標(biāo)準(zhǔn)化的本質(zhì)都是縮放和平移博杖。
區(qū)別：他們的區(qū)別直觀的說就是歸一化的縮放是 “拍扁” 統(tǒng)一到區(qū)間（0-1），而標(biāo)準(zhǔn)化的縮放是更加 “彈性” 和 “動(dòng)態(tài)” 的筷登，和整體樣本的分布有很大的關(guān)系剃根。

下文統(tǒng)一使用歸一化。

這篇文獻(xiàn)是我在查詢歸一化方法的時(shí)侯前方，Github上一個(gè)網(wǎng)友推薦《Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data》狈醉，翻譯過來就是：
《Illumina 高通量 RNA-seq 數(shù)據(jù)差異分析的歸一化方法比較 》。

下面開始介紹這篇2015年發(fā)表在 BMC Bioinformatics 文獻(xiàn)惠险。

如果覺得前面內(nèi)容太多苗傅，可以直接跳到最后看總結(jié)。

一班巩、背景和目的

1.1 背景：

RNA-seq 技術(shù)的快速發(fā)展和測(cè)序成本的降低使其成為一種廣泛應(yīng)用的基因表達(dá)定量技術(shù)渣慕。 由于歸一化在RNA-seq 數(shù)據(jù)分析中的重要性，人們提出了各種歸一化方法抱慌。
歸一化方法：
非豐度估計(jì)）的歸一化方法（non-abundance normalization
1. RC（row count）：每個(gè)基因的原始計(jì)數(shù)是所有run基因計(jì)數(shù)的總和逊桦。
2. UQ（upper quartile）：上四分位數(shù)是通過對(duì)所有樣本的基因計(jì)數(shù)應(yīng)用0.75的上四分位數(shù)來計(jì)算的，主要在芯片測(cè)序數(shù)據(jù)中使用抑进。
3. Med（median）：中位數(shù)計(jì)算為所有樣本基因計(jì)數(shù)的中位數(shù)强经。
4. TMM（Trimmed mean of M-values normalization）：M值的 trim 均值是一種用于RNA-seq 數(shù)據(jù)差異表達(dá)分析的標(biāo)度歸一化方法。 這種歸一化方法是在R包 edgeR中實(shí)現(xiàn)的单匣。 使用包中的 CalcNormFactors 函數(shù)計(jì)算縮放因子（scaling factors）夕凝，然后通過將基因計(jì)數(shù)除以每次運(yùn)行的每個(gè)縮放因子來獲得重新縮放的基因計(jì)數(shù)。 TMM是所有樣品 重新縮放(rescaled) 基因計(jì)數(shù)的總和户秤。
5. DESeq：DESeq是一種基于負(fù)二項(xiàng)分布模型的差異基因表達(dá)分析方法码秉， 方差和均值通過局部回歸聯(lián)系起來，并給出了一個(gè)也給出比例因子(scale factors) 的實(shí)現(xiàn)鸡号。 它在DESeq包中转砖，通過 EstimateSizeFactorsFormatrix函數(shù)，可以計(jì)算每次運(yùn)行的縮放因子鲸伴。 將基因計(jì)數(shù)除以每個(gè)標(biāo)度因子后府蔗，DESeq值被計(jì)算為所有樣品重新縮放基因計(jì)數(shù)的總和。
6. Q (quantiles)：分位數(shù)以前被用來歸一化數(shù)組之間的單通道或A-value 芯片數(shù)據(jù)汞窗。 R包limma中的normalizequantiles函數(shù)將矩陣的列歸一化為具有相同的分位數(shù)姓赤。 這里，我們將函數(shù)輸出的總值設(shè)置為分位數(shù)的歸一化值仲吏。
7. RPKM：這種方法通過對(duì)總轉(zhuǎn)錄本長度和測(cè)序 reads 數(shù)進(jìn)行歸一化不铆，從RNA-seq數(shù)據(jù)中量化基因表達(dá)蝌焚。 RPKM值可以使用以下定義輕松計(jì)算：

8. ERPKM:：RPKM的變形體，采用effective transcript length誓斥，但是作用不大只洒。由于reads的長度不為零，而reads 概率取決于有效長度劳坑，我們使用有效reads長度計(jì)算了每千貝每百萬映射reads的有效轉(zhuǎn)錄本毕谴。

豐度估計(jì)的歸一化方法（abundance normalization）：使用機(jī)器學(xué)習(xí)算法進(jìn)行豐度估計(jì)
     RSEM: 不同于以往的歸一化方法。 提出了一種結(jié)合期望最大化算法的有向圖模型來估計(jì)豐度距芬。 RSEM提供了一個(gè)從RNA-SEQ數(shù)據(jù)中量化基因豐度的軟件包涝开，因此我們通過準(zhǔn)備參考轉(zhuǎn)錄本數(shù)據(jù)和輸入RNA-SEQ數(shù)據(jù)計(jì)算RSEM值來生成參考指數(shù)。 
     Sailfish：在豐度估計(jì)中框仔，被介紹為無比對(duì)忠寻。它利用K-MER的概念對(duì)RNA-Seq reads進(jìn)行索引和計(jì)數(shù)。 在這里存和，我們使用偏置后的估計(jì)K-MERS數(shù)作為估計(jì)計(jì)數(shù)。 

1.2 目的：

對(duì)目前存在的歸一化方法進(jìn)行比較衷旅，以便為將來的實(shí)驗(yàn)選擇最合適的方法產(chǎn)生合適的指導(dǎo)方針捐腿。

二、方法和材料

2.1 方法

==通俗解釋：==
用 實(shí)驗(yàn)測(cè)量的 qRT-PCR值 和每種歸一化計(jì)算出的RNA-seq的豐度（或者可以理解為歸一化后的 reads count）預(yù)測(cè)進(jìn)行比較, 用 可以評(píng)估 數(shù)據(jù)相關(guān)性 的統(tǒng)計(jì)方法柿顶，來評(píng)估每種歸一化方法計(jì)算結(jié)果的準(zhǔn)確性茄袖。

【這里我是存在疑惑的，qRT-PCR實(shí)驗(yàn)中也存在誤差嘁锯，其數(shù)據(jù)結(jié)果是否可行宪祥，重復(fù)了幾次，文中作者僅使用原作者數(shù)據(jù)進(jìn)行的評(píng)估家乘，自己沒有實(shí)際做實(shí)驗(yàn)得到 qRT-PCR 數(shù)據(jù)】

==學(xué)術(shù)解釋：==
用Shapiro-Wilk正態(tài)性檢驗(yàn)QRT-PCR值的分布和所有的歸一化結(jié)果蝗羊。
根據(jù)P值<0.05的檢測(cè)結(jié)果，QRTPCR值和歸一化結(jié)果均不呈正態(tài)分布仁锯。
為了對(duì)數(shù)據(jù)進(jìn)行描述耀找，我們使用Spearman的秩相關(guān)系數(shù)，通過計(jì)算每種歸一化方法的RNA-Seq豐度預(yù)測(cè)與測(cè)量的QRT-PCR值之間的相似性來評(píng)估性能业崖。

Spearman相關(guān)系數(shù)作為一種非參數(shù)方法來度量?jī)蓚€(gè)變量之間的線性相關(guān)性野芒。 它被計(jì)算為數(shù)據(jù)秩上的皮爾遜相關(guān)系數(shù)。 對(duì)于一組大小為n的基因和相應(yīng)的n個(gè)原始數(shù)據(jù)双炕，變量 x 為 QRT-PCR 基因表達(dá)值狞悲，變量 y 為歸一化方法的結(jié)果，相關(guān)Rs 使用下面的公式計(jì)算妇斤。

Spearman相關(guān)系數(shù)將產(chǎn)生+1和-1之間的值摇锋，其中+1表示總的正相關(guān)丹拯，0表示不相關(guān)，而-1表示總的負(fù)相關(guān)乱投。 最接近+1或-1的值表示最高的相關(guān)性咽笼，因此也是最好的歸一化結(jié)果。

2.2 材料

高通量 RNA-seq 數(shù)據(jù)收集自NCBI數(shù)據(jù)庫中的SRA (Sequence Read Archive) 數(shù)據(jù)庫戚炫，如果自己的實(shí)驗(yàn)室沒有條件進(jìn)行測(cè)序剑刑，又想做數(shù)據(jù)分析，可以在這個(gè)數(shù)據(jù)庫下載數(shù)據(jù)双肤，自己處理和分析施掏。
腦組織(HBR)和組織類型混合物(UHR)的原始 RNA-seq 數(shù)據(jù)
35 bp reads length SRA010153.1
76 bp reads length SRA039286
一般現(xiàn)在的測(cè)序數(shù)據(jù)都是reads length 或者說插入片段長度都是 150或者151 bp的雙端測(cè)序數(shù)據(jù)。

結(jié)果

3.1 歸一化方法的比較

1. Spearman相關(guān)分析顯示茅糜，無論reads 長度如何七芭，RC，UQ蔑赘，Med狸驳，TMM，DESeq 和Q 都沒有明顯改善基因表達(dá)正乘跞化耙箍。
2. ERPKM并沒有取得比RPKM更好的結(jié)果。 使用有效的轉(zhuǎn)錄本長度顯然不能改善歸一化結(jié)果

• 在實(shí)驗(yàn)結(jié)果中酥馍，RPKM結(jié)合Salifish 的歸一化方法幾乎可以取代qRT-PCR測(cè)量辩昆。 而對(duì)于76 bp 序列數(shù)據(jù)，采用無豐度估計(jì)歸一化方法的RC獲得了最好的結(jié)果旨袒，其次是相關(guān)性相似的RSEM汁针； Salifish方法產(chǎn)生了更差的相關(guān)值。

• 從比較結(jié)果來看砚尽，歸一化方法并不是所有序列數(shù)據(jù)都必須的施无。 TMM、DESeq和Q等樣本間歸一化方法無論 reads 長度如何都不能顯著提高基因表達(dá)必孤，但當(dāng)比對(duì)精度較低時(shí)帆精，RPKM可能更有效。

• 對(duì)于35 bp的reads 數(shù)據(jù)隧魄，在兩種豐度估計(jì)歸一化方法中卓练，Salifish方法與RPKM結(jié)合的歸一化結(jié)果比RSEM更好，因?yàn)樗鼪]有比對(duì)购啄，也是一種相當(dāng)高效的組合襟企。 然而，==當(dāng)比對(duì)精度較高時(shí)狮含，RC似乎足夠用于實(shí)際實(shí)驗(yàn)中的基因表達(dá)計(jì)算顽悼。 ==

• 表1給出了八種非豐度估計(jì)歸一化方法（未應(yīng)用豐度估計(jì)歸一化）的Spearman相關(guān)系數(shù)結(jié)果曼振。

  
  

• 表2給出了結(jié)合RC和RPKM的兩種豐度估計(jì)方法（RSEM和Salifish方法）的Spearman相關(guān)系數(shù)結(jié)果。

  
  

3.2 添加 poly-A tails 的比較

• 換句話說蔚龙，通過增加Poly-A尾長冰评，可以將更多的 reads 映射到一個(gè)參考轉(zhuǎn)錄本。
• 總之木羹，選擇合適的Poly-A尾部長度可以改善差異分析甲雅； 然而，基于本研究中觀察到的最小效應(yīng)坑填，Poly-A尾長可能可以忽略不計(jì)抛人。

四、總結(jié)

1. 對(duì)于reads長度為35 bp的樣本脐瑰，在8種非豐度估計(jì)歸一化方法中妖枚，RPKM的相關(guān)值高于RC、UQ苍在、MED绝页、TMM和Q，證明在歸一化中考慮轉(zhuǎn)錄本長度是相當(dāng)有效的寂恬。

   • 這個(gè)結(jié)果在我們進(jìn)行短序列測(cè)序抒寂，例如miRNA測(cè)序的時(shí)侯就可以作為參考，一般的miRNA定量用的是RPM掠剑，華中農(nóng)大夏瑞老師課題組開發(fā)的sRNAminer 中用的是RP10M （reads per 10 million），也就是說郊愧，RPM主要應(yīng)用于sRNA（sRNA長度變化不大）朴译，來消除測(cè)序深度bias；而RPKM/FPKM應(yīng)用范圍會(huì)更廣泛属铁，用于同時(shí)消除長度及測(cè)序深度這兩種bias眠寿。

2. 通過使用有效轉(zhuǎn)錄本長度，ERPKM與RPKM相比并沒有改善歸一化結(jié)果焦蘑。 在結(jié)合豐度估計(jì)歸一化方法后盯拱，對(duì)歸一化結(jié)果進(jìn)行了改進(jìn)。 特別是RPKM和Salifish結(jié)合例嘱，我們建議研究者在未來的分析中作為一種歸一化方法狡逢，幾乎可以取代QRT-PCR，因?yàn)橛^察到了近0.8的相關(guān)性拼卵。 而且奢浑，Salifish是無比對(duì)的，比RSEM更節(jié)省時(shí)間腋腮。 RSEM也產(chǎn)生了良好的效果雀彼。

   • 這估計(jì)是ERPKM沒有RPKM流行的原因壤蚜，因?yàn)閑ffective transcript length 并沒用什用，現(xiàn)在TPM準(zhǔn)確性更高徊哑，使用的更頻繁袜刷，TPM和RPKM計(jì)算公式相同。
     • TPM與RPKM的區(qū)別：唯一的不同是計(jì)算操作的順序莺丑，TPM是先去除了基因長度的影響著蟹，而RPKM/FPKM是先去除測(cè)序深度的影響，TPM實(shí)際上改進(jìn)了RPKM方法在不同樣品間定量的不準(zhǔn)確性窒盐。

3. 然而草则，對(duì)于reads長度為76bp的樣本，沒有一種歸一化方法改善相關(guān)結(jié)果蟹漓。 因此炕横，我們得出結(jié)論，當(dāng)比對(duì)精度較高時(shí)葡粒，RC足以用于實(shí)際實(shí)驗(yàn)中的基因表達(dá)計(jì)算份殿。 另外，Poly-A尾的影響實(shí)驗(yàn)是通過在轉(zhuǎn)錄本數(shù)據(jù)中添加腺嘌呤（0-25）嗽交，結(jié)果表明卿嘲，選擇適當(dāng)?shù)膒oly-A尾部長度可以改善差異分析，但在本研究中似乎沒有影響夫壁。 因此拾枣，我們建議研究者在基因差異表達(dá)分析的歸一化步驟中不需要考慮ploy-a tails。

   • 現(xiàn)在有參考基因組的物種比對(duì)精度都比較高盒让，因?yàn)镽NA-seq測(cè)序質(zhì)量已經(jīng)不是問題梅肤，因此可以使用 Row count 作為基因表達(dá)的計(jì)算，然后計(jì)算一下Fold Change邑茄，再取對(duì)數(shù)姨蝴，比如我現(xiàn)在的這個(gè)課題就完全適用。
     計(jì)算方法如下：
     E = mean(group1)
     B=mean(group2)
     FC = (E-B) / min(E,B)

參考文獻(xiàn)

• Li P, Piao Y, Shon HS, Ryu KH. Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data. BMC Bioinformatics. 2015 Oct 28;16:347. doi: 10.1186/s12859-015-0778-7. PMID: 26511205; PMCID: PMC4625728.
• deweylab/RSEM: RSEM: accurate quantification of gene and isoform expression from RNA-Seq data (github.com)
• 【生信】【轉(zhuǎn)錄組】RSEM轉(zhuǎn)錄本定量 - 知乎 (zhihu.com)

本文由mdnice多平臺(tái)發(fā)布