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/sh-clusters-IFN.txt2.csv2.bed.gt5tc.gt0.25TtoC2.plus.noIGG",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/sh-clusters-IL4.txt2.csv2.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=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)]) ms2=ms; 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)]) ms2=append(ms2,ms); 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)]) ms2=append(ms2,ms); summary((ms2)) modscore_analyis(ms,"0h_all_reps") 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)]) ms2=ms; 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)]) ms2=append(ms2,ms); 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)]) ms2=append(ms2,ms); summary((ms2)) modscore_analyis(ms,"2h_all_reps") 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)]) ms2=ms; 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)]) ms2=append(ms2,ms); 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)]) ms2=append(ms2,ms); summary((ms2)) modscore_analyis(ms,"6h_all_reps")