R二項分布檢驗與FDR校正

R語言二項分布檢驗與FDR校正

二項分布是重復n次的實驗，且每次實驗都是獨立的，只有兩種結果，并且相互對立的，生活中最常見的是投硬幣~~~在生物領域內也有很多符合此類分布的，如二倍體動物等位基因，來源于父本和母本的重組等。具體公式什么的就不寫了，寫個關于ASE的例子吧。

[plain] view plain copy
    # cat binom.r | R --slave --args <file>  
    args <- commandArgs()  
    fa <- read.table(args[4], header=FALSE, sep="\t")  
    n1 = fa$V8  
    n2 = fa$V12  
    len = length(n1)  
    pv = numeric(len)  
    for(i in 1:len){  
    pv[i] = 0  
    if(n1[i] > n2[i]) pv[i]=binom.test(n1[i], n1[i]+n2[i], p=1/2, alternative="greater")$p.value else pv[i]=binom.test(n2[i], n1[i]+n2[i], p=1/2, alternative="greater")$p.value  
    }  
    qv <- p.adjust(pv, method="fdr")#fdr校正  
    fa$pv = formatC(pv, digits=4)  
    fa$fdr = formatC(qv, digits=4)  
    write.table(fa, file=paste(args[4],".out",sep=""), sep="\t", col.names=FALSE, row.names=FALSE, quote=FALSE)  

其實也不算是什么例子，就是把過程寫下來熟悉一下寫法而已。

通過控制FDR(False Discovery Rate)來決定P值的域值。假設你挑選了R個差異表達的基因，其中有S個是真正有差異表達的，另外有V個其實是沒有差異表達的，是假陽性的。實踐中希望錯誤比例Q＝V/R平均而言不能超過某個預先設定的值（比如0.05），在統計學上，這也就等價于控制FDR不能超過5％。對所有候選基因的p值進行從小到大排序，則若想控制fdr不能超過q，則只需找到最大的正整數i，使得 p(i)<= (i*q)/m。然后，挑選對應p(1)，p(2)，……，p(i)的基因做為差異表達基因，這樣就能從統計學上保證fdr不超過q。