淺談多重檢驗校正FDR

 Posted: 四月 12, 2017  Under: Basic  By Kai  no Comments

例如，在我們對鑑定到的差異蛋白做GO功能註釋後，通常會計算一個p值。當某個蛋白的p值小於0.05（5%）時，我們通常認為這個蛋白在兩個樣本中的表達是有差異的。但是仍舊有5%的概率，這個蛋白並不是差異蛋白。那麼我們就錯誤地否認了原假設（在兩個樣本中沒有差異表達），導致了假陽性的產生（犯錯的概率為5%）。

如果檢驗一次，犯錯的概率是5%；檢測10000次，犯錯的次數就是500次，即額外多出了500次差異的結論（即使實際沒有差異）。 為了控制假陽性的次數，於是我們需要對p值進行多重檢驗校正，提高閾值。

第一種方法Bonferroni，最簡單嚴厲的方法。
例如，如果檢驗1000次，我們就講閾值設定為5% / 1000 = 0.00005；即使檢驗1000次，犯錯誤的概率還是保持在N×1000 = 5%。最終使得預期犯錯誤的次數不到1次，抹殺了一切假陽性的概率。但是該方法雖然簡單，但是檢驗過於嚴格，導致最後找不到顯著表達的蛋白（假陰性）。

第二種方法FDR（False Discovery Rate）
相對Bonferroni來說，FDR用比較溫和的方法對p值進行了校正。其試圖在假陽性和假陰性間達到平衡，將假/真陽性比例控制到一定範圍之內。例如，如果檢驗1000次，我們設定的閾值為0.05（5%），那麼無論我們得到多少個差異蛋白，這些差異蛋白出現假陽性的概率保持在5%之內，這就叫FDR＜5%。

那麼我們怎麼從p value 來估算FDR呢，人們設計了幾種不同的估算模型。其中使用最多的是Benjamini and Hochberg方法，簡稱BH法。雖然這個估算公式並不夠完美，但是也能解決大部分的問題，主要還是簡單好用！

FDR的計算方法
除了可以使用excel的BH計算方法外，對於較大的資料，我們推薦使用R命令p.adjust。 p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods # c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", # "fdr", "none")
 

我們還可以從R命令p.adjust的原始碼，瞭解其執行的機制是什麼。

> p.adjust
function (p, method = p.adjust.methods, n = length(p)){ method <- match.arg(method) if (method == "fdr") method <- "BH" nm <- names(p) p <- as.numeric(p) …… BH = { i <- lp:1L o <- order(p, decreasing = TRUE) ro <- order(o) pmin(1, cummin(n/i * p[o]))[ro] } …… p0 }

其實該函式表達的意思是這樣的：

1. 我們將一系列p值、校正方法（BH）以及所有p值的個數（length(p)）輸入到p.adjust函式中。
2. 將一系列的p值按照從大到小排序，然後利用下述公式計算每個p值所對應的FDR值。 公式：p * (n/i)， p是這一次檢驗的p value，n是檢驗的次數，i是排序後的位置ID（如最大的P值的i值肯定為n，第二大則是n-1，依次至最小為1）。
3. 將計算出來的FDR值賦予給排序後的p值，如果某一個p值所對應的FDR值大於前一位p值（排序的前一位）所對應的FDR值，則放棄公式計算出來的FDR值，選用與它前一位相同的值。因此會產生連續相同FDR值的現象；反之則保留計算的FDR值。
4. 將FDR值按照最初始的p值的順序進行重新排序，返回結果。

參考：http://www.omicshare.com/forum/thread-173-1-1.html