modscore_analyis=function(ms2,id){
m=ms2
m[m>1]=1 # clip modscore to 1
#o=0.0005
o=0.001 #shift interval increment
st=0.02 # interval length
#st=0.01
l8=NULL;
lp8=NULL;
for (x in seq(0.5,1-st,by=o)){
y=m[(m>x)&(m<=x+st)];
l8=c(l8,x+0.5*st)
#print (c(x,x+st))
te=ad.test(y,"punif",x,x+st);
lp8=c(lp8,te$p.value);
}

topuniform=max(which(lp8>0.05)); #highest modscore that is not uniformly distributed
print (paste("modscore_cutoff for",id,l8[topuniform],lp8[topuniform],sep=" "));
x=cbind(l8,lp8)
colnames(x)=c("modscore","pvalue_AD_uniform_test");
write.table(x,paste("pvalue_AD_uniform_test_",id,".txt",sep=""),row.names=F)
png(paste("modscore-hist-",id,".png",sep=""))
hist(m,breaks=70,main=id,xlab="modscore")
dev.off()

png(paste("pvalue_AD_uniform_test_",id,".png",sep=""))
plot(l8,-log10(lp8),col="red",pch=20,cex=1,main=id,xlab="modscore",ylab="-log10(pvalue_AD_uniform_test)",xaxt="n")
axis(1,at=seq(0.5,1,by=0.05))
abline(h=-log10(0.05))
text(0.52,-log10(0.05)+0.14,"-log10(0.05)",cex=0.6)
dev.off()
}

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/res/recommended_settings2/sh-clusters-IFN.txt2.csv2.bed.gt5tc.gt0.25TtoC2.plus.noIGG.anno.csv",sep="\t")
ms=unlist(strsplit(as.character(t[,7]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,7]),"_"))[seq(8,length(ms),by=10)])
summary((ms))
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#  0.5000   0.8259   0.9443   0.9060   0.9979 116.8000
modscore_analyis(ms,"IFN")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/res/recommended_settings2/sh-clusters-IL4.txt2.csv2.bed.gt5tc.gt0.25TtoC2.plus.noIGG.anno.csv",sep="\t")
ms=unlist(strsplit(as.character(t[,7]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,7]),"_"))[seq(8,length(ms),by=10)])
summary((ms))
modscore_analyis(ms,"IL4")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-0hrep1.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"0hrep1")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-0hrep2.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"0hrep2")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-0hrep3.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"0hrep3")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-2hrep1.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"2hrep1")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-2hrep2.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"2hrep2")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-2hrep3.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"2hrep3")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-6hrep1.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"6hrep1")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-6hrep2.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"6hrep2")

t=read.table("/data/images/proton/DKlab/mr/parclip/paralyzer/PARalyzer_v1_1b/sh-clusters-6hrep3.txt2.csv3.bed.gt5tc.gt0.25TtoC2.plus.noIGG",sep="\t")
ms=unlist(strsplit(as.character(t[,4]),"_"))
ms=as.numeric(unlist(strsplit(as.character(t[,4]),"_"))[seq(8,length(ms),by=11)])
summary((ms))
modscore_analyis(ms,"6hrep3")