Thursday, July 3, 2014

robust.fdr on p-values



rm(list=ls())
source("robust-fdr.R")

setwd("~... ... CNV/pvalues,qvalues/RTFET")
infiles = list.files(, pattern="txt")
  
for( i in 1:length(infiles)){ 
  tb = read.table(infiles[i], sep="\t", header=F)
  
  tmp1 = tb[tb$V4<1.0, ];
  if( length(tmp1[,1]) >2  ) {
    hist(tmp1$V2)
    ret = robust.fdr(tmp1$V4, sides=1, discrete=T )
    tmp1$q1 = ret$q
    tb$q1 = tmp1$q1[match(tb$V1, tmp1$V1)]
  } else {
    tb$q1 = NA;
  }
  summary(tb)
  tb[tb$q1<0.15 & (!is.na(tb$q1)), ]
  
  outfile = paste("output_robustFDR/", infiles[i], ".csv", sep='')
  write.csv(tb, outfile,, row.names=F, quote=F)
}

No comments:

Post a Comment