# R code
dir<-"~/R/chongli/DMR_win/single_end/test/"

head<-c("5E.5A","5E.6B","6D.5A","6D.6B")

chr<-paste0(".chr",1:5)

suffix<-c(".-.400",".+.400",".-.800",".+.800")

count<-array(0,c(4,5,4))
dimnames(count)[[1]]<-head
dimnames(count)[[2]]<-chr
dimnames(count)[[3]]<-suffix

for(h in 1:4){
  for(n in 1:5){
    for(k in 1:4){
      filename<-paste0(dir,head[h],chr[n],suffix[k])
      #if(is.na(file.info(filename)[1])==FALSE){  
        bonferroni<-read.table(filename,header=T)
        fdrtest<-bonferroni[order(bonferroni[,4]),]   #sort by the pvalue column
        fdr.pvalue<-fdrtest[,4] 
        q.value<-0.1                             #q.value to be determined

        i<-1
        repeat{
	    if( fdr.pvalue[i] > q.value*i/length(fdr.pvalue) | i==length(fdr.pvalue)) break  # loop stop condition: fdr algorithm or meet the end
           else  i<-i+1
          }
        if(i!=1){
             fdrtest<-fdrtest[1:i,]                                                               #only keep the pre-i lines
             fdrtest<-fdrtest[order(fdrtest[,1]),]
	}else{
             fdrtest<-fdrtest[i,]
        }
        count[h,n,k]<-i
        filenames<-paste0(head[h],chr[n],suffix[k],".fdr0.1")
        write.table(fdrtest,file=filenames,quote=F,row.names=F,sep="\t")  #sort by chr then address
      #}
    }
  }
}
sink("total.info")
count
sink()