寫在前面的廢話
最近在處理一批RNA-seq的數(shù)據(jù)，里面混入了spike-in夺巩。利用spike-in矯正之后贞让，樣本A的基因表達(dá)量普遍比樣本B的基因表達(dá)量高3-5倍，這和我所熟知的背景知識(shí)是一致的柳譬。
但是當(dāng)我使用DESeq2對(duì)基因表達(dá)量進(jìn)行差異表達(dá)分析時(shí)喳张，上調(diào)的基因和下調(diào)的基因竟然相差無(wú)幾，都有兩三千個(gè)……這不符合邏輯啊美澳，前前后后思考了一遍销部，發(fā)現(xiàn)是我對(duì)DESeq2的理解太淺的緣故。
太長(zhǎng)不看系列
spike-in是已知的其他物種基因組制跟，可以對(duì)基因表達(dá)數(shù)據(jù)進(jìn)行絕對(duì)值的矯正
個(gè)人認(rèn)為舅桩，F(xiàn)PKM/RPKM/TPM都是基因的相對(duì)表達(dá)值
DESeq2會(huì)對(duì)測(cè)序深度進(jìn)行矯正，將普遍高表達(dá)的樣本A看作是測(cè)序深度過高所導(dǎo)致雨膨，從而影響差異基因的篩選
廢話超多系列
RNA的spike-in是一組序列已知的RNA轉(zhuǎn)錄本擂涛，目前普遍使用的是ERCC( The External RNA Control Consortium)搞出來(lái)的一組RNA序列。當(dāng)然聊记，也可以使用序列相似度較低的物種的序列作為spike-in撒妈。如兩種常見的酵母，S.pombe和S. cerevisiae排监，他們的序列可以作為彼此的spike-in進(jìn)行后續(xù)矯正狰右。
通過在樣品制備過程中，混入指定數(shù)量的spike-in舆床，我們就可以知道不同樣本中的基因絕對(duì)比表達(dá)值棋蚌。
如等細(xì)胞數(shù)的樣本A和樣本B，在每個(gè)樣本中挨队，我加入了等量的spike-in谷暮。最后分析發(fā)現(xiàn)，spike-in占樣本A的1%瞒瘸，但是占樣本B的5%坷备。這表明樣本A的RNA表達(dá)量也許普遍比樣本B的表達(dá)量高五倍左右。
下面是一個(gè)簡(jiǎn)單的草圖情臭，希望可以幫助理解省撑，左邊是細(xì)胞赌蔑，右邊是RNA
說(shuō)完了spike-in矯正的原理，接下來(lái)講講DESeq2文庫(kù)矯正的原理竟秫。
常用的normalization方法有FPKM/RPKM/TPM娃惯，但是這些方法只能對(duì)測(cè)序深度和基因長(zhǎng)度進(jìn)行矯正，沒有考慮樣本的文庫(kù)組成成分可能對(duì)表達(dá)量產(chǎn)生的影響肥败。這里舉一個(gè)例子：
看起來(lái)兩個(gè)不同的樣本之間除了基因A1CF趾浅，其余基因都是差異表達(dá)的。但事實(shí)上馒稍，這是由于樣本A中A2M表達(dá)量過高皿哨，導(dǎo)致樣本中其他基因的相對(duì)表達(dá)量較低。如果我們把兩個(gè)樣本中的A2M基因去除纽谒，重新計(jì)算FPKM证膨，我們會(huì)發(fā)現(xiàn)兩個(gè)樣本中其他五個(gè)基因的表達(dá)量其實(shí)相差無(wú)幾。
因此DESeq2需要解決兩個(gè)問題：
1鼓黔、樣本的測(cè)序深度
2央勒、樣本的文庫(kù)組成
這里以一個(gè)簡(jiǎn)單的例子講解一下DESeq2是如何解決這兩個(gè)問題的。
首先對(duì)樣本量進(jìn)行以自然對(duì)數(shù)為底數(shù)的log轉(zhuǎn)換：
DESeq2除了可以使用澳化，還可以使用和崔步，但因?yàn)镽默認(rèn)的log是以e為底數(shù)，因此這里使用缎谷。我曾經(jīng)畫圖時(shí)為了偷懶使用井濒，導(dǎo)致組會(huì)上被老板批評(píng)，因?yàn)楦闵诺那拜厒兤毡橄矚g使用或者進(jìn)行對(duì)數(shù)轉(zhuǎn)換慎陵。我太難了╥﹏╥...
注：我們可以看到眼虱，表達(dá)量為0的值在去對(duì)數(shù)之后，變成了負(fù)無(wú)窮席纽。
對(duì)每一行的值進(jìn)行平均值計(jì)算
第二列和第三列都比較好理解，第一列因?yàn)闃颖?的gene1的log值是-inf撞蚕，因此gene1這一列的平均值也是-inf润梯。這里還有一點(diǎn)值得關(guān)注，當(dāng)我們將gene3的平均表達(dá)的log轉(zhuǎn)換為正常數(shù)值時(shí)甥厦，≈73.7纺铭，而對(duì)基因表達(dá)值的原始矩陣計(jì)算得到的平均值為(33+55+200)/3=96。
我們可以看到96>73.7刀疙，因此這種取對(duì)數(shù)的方法可以使得基因更不容易受到異常值的影響舶赔。事實(shí)上這種取完log之后做平均的方法，計(jì)算的是幾何平均數(shù)
移除具有-inf值的基因
當(dāng)一個(gè)樣本或者多個(gè)樣本中基因的表達(dá)量為0時(shí)谦秧，這個(gè)基因就會(huì)被移除竟纳。
事實(shí)上撵溃，通過這個(gè)可以讓我們最后的基因表達(dá)矩陣更加集中于管家基因（house-keeping gene）
對(duì)基因的reads count取log并減去基因的log值的平均數(shù)，具體如下：
通過這個(gè)計(jì)算方式锥累，我們將得到每個(gè)樣本中g(shù)eneX的表達(dá)水平與 geneX在所有樣本平均表達(dá)水平的比值
計(jì)算步驟四所得的樣本表達(dá)矩陣中缘挑，每個(gè)樣本中的基因表達(dá)中位數(shù)
這里使用中位數(shù)，可以進(jìn)一步減少異常值對(duì)于scale factor(校正因子)的影響桶略。至于為什么中位數(shù)有這個(gè)好處语淘，我想這里應(yīng)該不需要詳細(xì)闡述了
將步驟5計(jì)算的每個(gè)樣本的中位數(shù)，進(jìn)行指數(shù)計(jì)算
通過該方法际歼，我們最終得到了樣本的校正因子(scale factor)
將原始樣本表達(dá)矩陣除以步驟6所得到的scale factor
以上七步就是DESeq2對(duì)于樣本文庫(kù)矯正的方法惶翻，
為避免文章太長(zhǎng)，對(duì)前文有所遺忘鹅心。我們?cè)谶@里再簡(jiǎn)單的敘述一下吕粗，DESeq2對(duì)于樣本文庫(kù)的矯正做了些什么：
1、使用log轉(zhuǎn)換巴帮，去除那些只在某幾個(gè)樣本中表達(dá)的基因溯泣，也減少了極端值對(duì)于最終結(jié)果的影響
2、取每個(gè)樣本的中位數(shù)榕茧，進(jìn)一步減少極端值的影響垃沦，使得結(jié)果更加關(guān)注于在多個(gè)樣本中穩(wěn)定表達(dá)的基因
通過對(duì)上述DESeq2的文庫(kù)normalization方法的學(xué)習(xí)，我了解到我的RNA-seq數(shù)據(jù)是受到了DESeq2自動(dòng)文庫(kù)矯正的影響用押。這個(gè)是可以關(guān)閉的
這可真是個(gè)bug啊肢簿，DESeq2明明是好心卻辦出錯(cuò)了事，無(wú)奈我只能使用logFoldChange和T檢驗(yàn)的方法篩選差異表達(dá)基因……
寫在后面的話
我在想蜻拨，DESeq2這么牛池充，作者一定很聰明，不可能沒想到這個(gè)問題缎讼。此外收夸，這么多人在使用DESeq2，一定也遇到過類似的問題……
果不其然血崭，DESeq2對(duì)于這個(gè)問題早有設(shè)計(jì)卧惜，我還是太naive了……
具體怎么實(shí)現(xiàn)，咱們下篇文章見(～￣▽￣)～
參考資料
Statquest：DESeq2, part1: Library Normalization
