edgeR and DESeq2, part2-獨(dú)立過濾

? 過濾低read計(jì)數(shù)的基因蔗崎，又名獨(dú)立過濾(Independent Filtering)庙睡。

? 每次我們做統(tǒng)計(jì)檢驗(yàn)诀艰，都有可能得出錯誤的結(jié)論讶踪。簡而言之，當(dāng)我們說p值< 0.05顯著躏将，那么我們也可以說5%的情況下我們會報(bào)假陽性硬毕。

? 當(dāng)我們處理1或2個(gè)基因的差異表達(dá)時(shí)呻引，這不是什么大問題，因?yàn)?次測試的5%是很小的吐咳，我們不太可能報(bào)告假陽性(false-positive)逻悠。然而，當(dāng)我們檢查基因組中的每一個(gè)基因(大約20000個(gè))韭脊，看看哪些基因在癌細(xì)胞中被調(diào)控失調(diào)時(shí)童谒，5%*20000=1000 false positive。好消息是沪羔，F(xiàn)DR和Benjamini-Hochberg方法彌補(bǔ)了這個(gè)問題,但是還存在問題饥伊，我們先看一個(gè)例子：

FDR過濾

假設(shè)我們有兩個(gè)獨(dú)立的分布，紅色曲線表示小鼠品系X的體重蔫饰，藍(lán)色曲線表示小鼠品系Y的體重琅豆。如果我們測量3只X品系的老鼠，那些值很大可能接近紅色分布中間篓吁。如果我們測量3只Y品系的老鼠趋距，那些值很大可能接近藍(lán)色分布中間。對這些體重進(jìn)行t檢驗(yàn)將導(dǎo)致p值< 0.05越除，我們將正確地得出測量結(jié)果來自不同的分布。

image-20210106110808050.png

然而，我們時(shí)不時(shí)會得到重疊的值摘盆。在這種情況下翼雀，p值將會很大。這是一個(gè)“假陰性”孩擂。

image-20210106110848581.png

現(xiàn)在狼渊，讓電腦從這些分布中抽取1000個(gè)樣本(3個(gè)X，3個(gè)Y)做1000個(gè)t檢驗(yàn)类垦。畫1000個(gè)p值的直方圖狈邑，我們得到949個(gè)真陽性(p值<0.05),5個(gè)假陽性(p值>0.05),我們做的每一個(gè)測試都應(yīng)該是"真陽性"并且p值< 0.05，這是因?yàn)槊恳粋€(gè)測試都使用兩個(gè)不同分布的樣本蚤认。

image-20210106111606005.png

讓我們通過添加沒有區(qū)別的test來使它更真實(shí)米苹。我們將從已經(jīng)使用兩個(gè)不同的分布計(jì)算的1000個(gè)p值開始。然后再加上1000個(gè)來自相同分布的樣本值砰琢，這些p值應(yīng)該是> 0.05蘸嘶，但偶爾(5%的時(shí)間)，我們會得到p值< 0.05陪汽。

image-20210106111840928.png

一共有993個(gè)p值小于0.05训唱，949個(gè)真陽性來自第一組p值，44個(gè)假陽性來自第二組p值挚冤，因?yàn)閮H僅只有4%的p值<0.05是假陽性况增，我們不需要使用FDR。但這只是因?yàn)槲覀兙幵炝藬?shù)據(jù)训挡，如果是真的數(shù)據(jù)澳骤，我們不知道百分比，所以我們要用FDR舍哄。

image-20210106112434174.png

在使用FDR校正后宴凉，846個(gè)FDR調(diào)整后的p值仍然< 0.05。827個(gè)正陽性值剩下表悬，為原949的89%弥锄，19個(gè)假陽性值剩下，占846個(gè)的2%蟆沫。

image-20210106112928967.png

現(xiàn)在籽暇，讓它更像RNA-Seq，p值的數(shù)目增加到6000饭庞。意味著戒悠，1000個(gè)樣本來自不同的分布，5000個(gè)樣本來自相同的分布舟山，在大多數(shù)情況下來自相同樣本應(yīng)該會給出較大的p值绸狐。

image-20210106113246830.png

有1215個(gè)p值<0.05,949個(gè)p值是真陽性值卤恳，266個(gè)p值是假陽性值，22%的p值是假陽性值寒矿。

image-20210106113657633.png

FDR校正后突琳，僅僅剩下256個(gè)p值<0.05,250個(gè)真陽性值剩下，占949的26%符相。6個(gè)假陽性值剩下拆融，占256的2%。FDR在限制“顯著結(jié)果”中的假陽性數(shù)量方面做得很好啊终，但在保持真實(shí)陽性方面做得并不出色镜豹。

image-20210106114134742.png

現(xiàn)在讓我們把樣本量增加到11000，1000個(gè)樣本來自不同的分布蓝牲，10000個(gè)樣本來自相同的分布趟脂，然后做11000次t檢驗(yàn)，得到10000個(gè)p值搞旭。有1430個(gè)p值<0.05,949個(gè)p值是真陽性值散怖，481個(gè)p值是假陽性值，34%的p值是假陽性值肄渗。FDR校正后镇眷，僅僅剩下56個(gè)p值<0.05,54個(gè)真陽性值剩下，占949的6%翎嫡。2個(gè)假陽性值剩下欠动，每當(dāng)我們增加樣本的檢驗(yàn)數(shù)目時(shí)，通過FDR校正的真陽性的（p小于0.05）數(shù)目都會減少：

image-20210106115039256.png

下圖中惑申，綠線代表通過FDR校正后具伍，真陽性p值<0.05的數(shù)量。橙線代表通過FDR校正后圈驼，假陽性p值<0.05的數(shù)量人芽。這張圖顯示，盡管FDR可以控制假陽性的比率绩脆，隨著檢驗(yàn)數(shù)目的增多萤厅，真陽性的比例卻在下降，這也表明Benjamini Hochberg方法還有改進(jìn)的空間!

image-20210106115508646.png

edgeR和DESeq2過濾

? 一般的想法是靴迫，具有超低read計(jì)數(shù)的基因不能提供有用信息惕味，因此，它們可以從數(shù)據(jù)集中刪除玉锌。換句話說名挥，即使這些基因在生物學(xué)上是有趣的，如果在一種樣本類型中只有1或2個(gè)轉(zhuǎn)錄本主守，而在另一種樣本類型中只有3或4個(gè)轉(zhuǎn)錄本禀倔，就很難得到準(zhǔn)確的read計(jì)數(shù)榄融。

edgeR過濾

在做任何事情之前，edgeR建議去除所有的基因蹋艺，除了那些在兩個(gè)樣品或更多樣本中CPM>1的基因剃袍。CPM=Counts Per Million,它彌補(bǔ)了文庫之間read深度的差異。

計(jì)算CPM:

計(jì)算公式：

image-20210106120653891.png

image-20210106120831864.png

現(xiàn)在我們有了所有樣本中所有基因的CPM值捎谨，讓我們?nèi)コ械幕颍四切┰趦蓚€(gè)樣品或更多樣本中CPM>1的基因憔维。

image-20210106121209831.png

edgR的方法很簡單涛救，但是你應(yīng)該意識到測序深度會影響它。例如业扒，如果一個(gè)樣本有500萬reads检吆，CPM標(biāo)準(zhǔn)化因子=5000000/1000000=5,如果有5reads比對到一個(gè)基因，這個(gè)基因的CPM=5/5=1 CPM.如果一個(gè)樣本有8000萬reads程储，CPM標(biāo)準(zhǔn)化因子=80000000/1000000=80,此時(shí)1CPM=80reads蹭沛。

image-20210106121555961.png

另一方面，有時(shí)你需要一個(gè)更大的CPM閾值章鲤，例如摊灭，你有50000reads比對到一個(gè)樣本，標(biāo)準(zhǔn)化因子CPM=50000/1000000=0.05,如果你有一個(gè)read比對到一個(gè)基因上败徊，它將變成1/0.05 = 20 CPM帚呼。即使這個(gè)基因是在生物學(xué)相關(guān)的水平上轉(zhuǎn)錄，因?yàn)槟阒荒茏x到一個(gè)read皱蹦，它仍然存在很大的噪音煤杀。

edgeR中CPM閾值

我們得到一個(gè)很好的閾值(cut off)？我們通過真實(shí)的數(shù)據(jù)集來說明

• 我從一位同事那里得到一個(gè)數(shù)據(jù)集沪哺，每個(gè)樣本平均有2200萬reads沈自。(4個(gè)“野生型”及4個(gè)“敲除型”樣本)
• 我在沒有過濾單個(gè)基因的情況下對它進(jìn)行了edgeR，這將生成原始的p值
• 使用不同的CPM閾值過濾掉基因辜妓，然后矯正p值

我們使用不同的CPM閾值過濾掉基因枯途，然后矯正p值，繪制成圖嫌拣，其中x軸為最小CPM閾值柔袁，y軸表示經(jīng)過FDR校正后，p值<0.05的基因數(shù)量异逐。當(dāng)x=0時(shí)捶索，沒有過濾掉任何基因，當(dāng)x=1(推薦的閾值)時(shí)灰瞻，此時(shí)基因數(shù)量為3400腥例，因?yàn)橛泻芏鄏eads辅甥，建議的閾值太嚴(yán)格了。使用較低的閾值可以鑒別出大約200個(gè)更重要的基因燎竖。

image-20210106124901166.png

edgeR的主旨是：要小心璃弄，在計(jì)算p值后，嘗試不同的CPM閾值(cutoff)构回。

DESeq和edgeR的區(qū)別

1. edgeR查看單個(gè)樣本夏块，并確保至少有2個(gè)CMP大于閾值。相比之下纤掸，DESeq2查看某個(gè)基因所有樣本均一化reads的平均值脐供，>閾值便保留。此時(shí)你可能會想:“酷!!我用DESeq2的方法借跪，但是如果測量的基因有異常值呢?”

• DESeq2有一個(gè)異常值值檢測方法(我們將在另一個(gè)StatQuest中討論)政己，但它只在每個(gè)類別有兩個(gè)以上的樣本時(shí)才生效。

下圖是分別使用DESeq和edgeR對同一批數(shù)據(jù)進(jìn)行處理掏愁，它們都在同一區(qū)域達(dá)到峰值歇由，這兩種方法都會產(chǎn)生相似的閾值。現(xiàn)在讓我們看看另一個(gè)不同之處果港。

2.另一個(gè)不同是改變了x軸沦泌，DESeq2繪制了分位數(shù)與顯著基因數(shù)量的圖，而不是最小CPM閾值京腥。0%的基因低于閾值赦肃，20%的基因低于閾值，40%的基因低于閾值公浪，分位數(shù)是有用的他宛，因?yàn)檎缥覀兯吹降模珻PM依賴于測序深度欠气，但無論如何厅各，分位數(shù)總是分位數(shù)。無論庫中有800萬或8000萬reads预柒，10%的基因總是小于0.1分位數(shù)队塘。

image-20210106131846588.png

我們可以選擇分位數(shù)和最小CPM，

3.DESeq2對這些點(diǎn)擬合一條曲線宜鸯，DESeq2然后在擬合曲線上找到最大位置憔古。閾值是曲線上的最大位置，減去擬合曲線與原始值之間的標(biāo)準(zhǔn)差淋袖。換句話說鸿市，在峰值噪聲范圍內(nèi)的第一個(gè)分位數(shù)是CPM閾值。如果沒有原始值超過閾值，則不進(jìn)行過濾焰情。

image-20210106132824031.png

現(xiàn)在我們知道了edgeR和DESeq2是如何過濾基因的陌凳。

• edgeR是保留那些在2個(gè)或2個(gè)以上樣本中，CPM大于最小閾值的基因
• DESeq2保留那些平均CPM大于最小CPM的基因内舟，然后繪制顯著基因與分位數(shù)的散點(diǎn)圖合敦，找到擬合曲線，再用最大值減去噪聲验游，即是閾值

image-20210106133105442.png

建議

• 如果使用edgeR充岛，在計(jì)算p值后計(jì)算CPM閾值
• 應(yīng)用DEseq2的方法很容易找到最優(yōu)CPM，以edgeR的基因選擇標(biāo)準(zhǔn)
• 如果你這樣做了批狱，請確保你引用了這兩個(gè)出版物!

• 如果你使用DESeq2裸准，當(dāng)每個(gè)分類只有2個(gè)樣本時(shí)要小心異常值。

  
  
  image-20210106133521094.png