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。