set.seed(1);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=F);##name name leading/lagging aa_length total_mutaion hostport_posi mutation_posi
times<-10000;
tem5<-NULL;
for(i in c(1:times)){
	tem4<-NULL;
	for(j in c(1:dim(data)[1])){
			tem<-sample(c(1:as.numeric(data[j,4])),as.numeric(data[j,5]),replace=T);
			tem2<-sapply(split(tem,tem),length);
			tem3<-length(which(tem2>1));
			tem4<-rbind(tem4,tem3);
		}
	tem5<-cbind(tem5,tem4);
}
result<-cbind(as.character(data[,3]),as.numeric(data[,9]),tem5);
write.table(result,"99.02simulation.txt",append=F,quote=F,row.names=F,col.names=F,sep="\t");

set.seed(1);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=F);
times<-10000;
tem5<-NULL;
for(i in c(1:times)){
	tem4<-NULL;
	for(j in c(1:dim(data)[1])){
			tem<-sample(c(1:as.numeric(data[j,7])),as.numeric(data[j,5]),replace=T);
			tem2<-sapply(split(tem,tem),length);
			tem3<-length(which(tem2>1));
			tem4<-rbind(tem4,tem3);
		}
	tem5<-cbind(tem5,tem4);
}
result<-cbind(as.character(data[,3]),as.numeric(data[,9]),tem5);
write.table(result,"99.02simulation_extreme.txt",append=F,quote=F,row.names=F,col.names=F,sep="\t");

set.seed(1);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=F);
#data[which(data[,8]>-1 & data[,8]<data[,4]),8]<-data[which(data[,8]>-1 & data[,8]<data[,4]),4];
data[which(data[,8]==-1),8]<-data[which(data[,8]==-1),4];
#data[which(data[,8]==-1),8]<-data[which(data[,8]==-1),7];
#data<-data[which(data[,8]>-1),];
times<-10000;
tem5<-NULL;
for(i in c(1:times)){
	tem4<-NULL;
	for(j in c(1:dim(data)[1])){
			tem<-sample(c(1:as.numeric(data[j,8])),as.numeric(data[j,5]),replace=T);
			tem2<-sapply(split(tem,tem),length);
			tem3<-length(which(tem2>1));
			tem4<-rbind(tem4,tem3);
		}
	tem5<-cbind(tem5,tem4);
}
result<-cbind(as.character(data[,3]),as.numeric(data[,9]),tem5);
write.table(result,"99.02simulation_available.txt",append=F,quote=F,row.names=F,col.names=F,sep="\t");

####
result<-read.table(file="99.02simulation.txt",sep="\t",head=F);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=F);
times<-10000;

tiff(filename = paste("99.02simulation.tif",sep=""),width=178, height=178, units="mm", res=200);
layout(matrix(c(1:4),2,2,byrow = T));
picname<-c("a","b");
mutposi<-c("leading","lagging");
result2<-result[which(as.numeric(result[,2])>=200),];
data2<-data[which(data[,9]>=200),];
for(i in c(1:length(mutposi))){
	simdata<-result2[which(result2[,1]==mutposi[i]),];
	obsdata<-data2[which(data2[,3]==mutposi[i]),];
	sumobs<-dim(obsdata[which(obsdata[,6]>0),])[1];
	sumsim<-apply(simdata[,c(3:dim(simdata)[2])],2,function(x){length(which(as.numeric(x)>0))});
	ymax<-2000;
	xmin<--10;
	xmax<-80;
	hist(sumsim,breaks=30,xlab="",ylab="Frequency",main=paste("median=",median(sumsim),sep=""),xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",cex.lab=1.2,freq=T);
	axis(side=1,seq(0,xmax,20), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,20),at=seq(0,xmax,20), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,length.out=5), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,length.out=5)/times,at=seq(0,ymax,length.out=5), side=2, las=3, line=1, col="black",cex=1);	
	arrows(sumobs,ymax*0.8,sumobs,0,lwd=3);
	text(sumobs,ymax*0.9,paste("   = ",length(which(sumsim>sumobs))/times*100,"%",sep=""),adj=0.1);
	text(sumobs,ymax*0.9,"P",adj=1,font=3);
	text(sumobs,ymax,"Observed",adj=0.15);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0);		
	mtext("Number of genes  (length > 200 aa)", side=1, las=1, line=2, col="black",cex=1);
	mtext("which have hotspot muations' sites", side=1, las=1, line=3, col="black",cex=1);
	mtext(picname[i],side=3, las=1, line=1, col="black",cex=1.5,adj=-0.2);	
	box();
}


picname<-c("c","d");
result2<-result[which(as.numeric(result[,2])>=200),];
data2<-data[which(data[,9]>=200),];
for(i in c(1:length(mutposi))){
	simdata<-result2[which(result2[,1]==mutposi[i]),];
	obsdata<-data2[which(data2[,3]==mutposi[i]),];
	sumobs<-sum(obsdata[,6]);
	sumsim<-apply(simdata[,c(3:dim(simdata)[2])],2,function(x){sum(as.numeric(x))});
	ymax<-2000;
	xmin<--10;
	xmax<-80;
	hist(sumsim,breaks=30,xlab="",ylab="Frequency",main=paste("median=",median(sumsim),sep=""),xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",cex.lab=1.2,freq=T);
	axis(side=1,seq(0,xmax,20), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,20),at=seq(0,xmax,20), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,length.out=5), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,length.out=5)/times,at=seq(0,ymax,length.out=5), side=2, las=3, line=1, col="black",cex=1);	
	arrows(sumobs,ymax*0.8,sumobs,0,lwd=3);
	text(sumobs,ymax*0.9,paste("   = ",length(which(sumsim>sumobs))/times*100,"%",sep=""),adj=0.1);
	text(sumobs,ymax*0.9,"P",adj=1,font=3);
	text(sumobs,ymax,"Observed",adj=0.15);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0);		
	mtext("Number of hotspot muations' sites", side=1, las=1, line=2, col="black",cex=1);
	mtext("in genes  (length > 200 aa)", side=1, las=1, line=3, col="black",cex=1);
	mtext(picname[i],side=3, las=1, line=1, col="black",cex=1.5,adj=-0.2);	
	box();
}

dev.off();

####
result<-read.table(file="99.02simulation_extreme.txt",sep="\t",head=F);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=F);
times<-10000;

tiff(filename = paste("99.02simulation_extreme.tif",sep=""),width=178, height=178, units="mm", res=200);
layout(matrix(c(1:4),2,2,byrow = T));
picname<-c("a","b");
mutposi<-c("leading","lagging");
result2<-result[which(as.numeric(result[,2])>=200),];
data2<-data[which(data[,9]>=200),];
for(i in c(1:length(mutposi))){
	simdata<-result2[which(result2[,1]==mutposi[i]),];
	obsdata<-data2[which(data2[,3]==mutposi[i]),];
	sumobs<-dim(obsdata[which(obsdata[,6]>0),])[1];
	sumsim<-apply(simdata[,c(3:dim(simdata)[2])],2,function(x){length(which(as.numeric(x)>0))});
	ymax<-3000;
	xmin<--10;
	xmax<-250;
	hist(sumsim,breaks=30,xlab="",ylab="Frequency",main=paste("median=",median(sumsim),sep=""),xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",cex.lab=1.2,freq=T);
	axis(side=1,seq(0,xmax,50), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,50),at=seq(0,xmax,50), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,length.out=5), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,length.out=5)/times,at=seq(0,ymax,length.out=5), side=2, las=3, line=1, col="black",cex=1);	
	arrows(sumobs,ymax*0.8,sumobs,0,lwd=3);
	text(sumobs,ymax*0.9,paste("   = ",length(which(sumsim>sumobs))/times*100,"%",sep=""),adj=0.1);
	text(sumobs,ymax*0.9,"P",adj=1,font=3);
	text(sumobs,ymax,"Observed",adj=0.15);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0);		
	mtext("Number of genes  (length > 200 aa)", side=1, las=1, line=2, col="black",cex=1);
	mtext("which have hotspot muations' sites", side=1, las=1, line=3, col="black",cex=1);
	mtext(picname[i],side=3, las=1, line=1, col="black",cex=1.5,adj=-0.2);	
	box();
}


picname<-c("c","d");
result2<-result[which(as.numeric(result[,2])>=200),];
data2<-data[which(data[,9]>=200),];
for(i in c(1:length(mutposi))){
	simdata<-result2[which(result2[,1]==mutposi[i]),];
	obsdata<-data2[which(data2[,3]==mutposi[i]),];
	sumobs<-sum(obsdata[,6]);
	sumsim<-apply(simdata[,c(3:dim(simdata)[2])],2,function(x){sum(as.numeric(x))});
	ymax<-2000;
	xmin<--10;
	xmax<-500;
	hist(sumsim,breaks=30,xlab="",ylab="Frequency",main=paste("median=",median(sumsim),sep=""),xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",cex.lab=1.2,freq=T);
	axis(side=1,seq(0,xmax,100), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,100),at=seq(0,xmax,100), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,length.out=5), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,length.out=5)/times,at=seq(0,ymax,length.out=5), side=2, las=3, line=1, col="black",cex=1);	
	arrows(sumobs,ymax*0.8,sumobs,0,lwd=3);
	text(sumobs,ymax*0.9,paste("   = ",length(which(sumsim>sumobs))/times*100,"%",sep=""),adj=0.1);
	text(sumobs,ymax*0.9,"P",adj=1,font=3);
	text(sumobs,ymax,"Observed",adj=0.15);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0);		
	mtext("Number of hotspot muations' sites", side=1, las=1, line=2, col="black",cex=1);
	mtext("in genes  (length > 200 aa)", side=1, las=1, line=3, col="black",cex=1);
	mtext(picname[i],side=3, las=1, line=1, col="black",cex=1.5,adj=-0.2);	
	box();
}

dev.off();

####
result<-read.table(file="99.02simulation_available.txt",sep="\t",head=F);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=F);
times<-10000;

#data<-data[which(data[,8]>-1),];
#result<-result[which(data[,8]>-1),];

tiff(filename = paste("99.02simulation_available.tif",sep=""),width=178, height=178, units="mm", res=200);
layout(matrix(c(1:4),2,2,byrow = T));
picname<-c("a","b");
mutposi<-c("leading","lagging");
result2<-result[which(as.numeric(result[,2])>=200),];
data2<-data[which(data[,9]>=200),];
for(i in c(1:length(mutposi))){
	simdata<-result2[which(result2[,1]==mutposi[i]),];
	obsdata<-data2[which(data2[,3]==mutposi[i]),];
	sumobs<-dim(obsdata[which(obsdata[,6]>0),])[1];
	sumsim<-apply(simdata[,c(3:dim(simdata)[2])],2,function(x){length(which(as.numeric(x)>0))});
	ymax<-3000;
	xmin<--5;
	xmax<-40;
	hist(sumsim,breaks=30,xlab="",ylab="Frequency",main=paste("median=",median(sumsim),sep=""),xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",cex.lab=1.2,freq=T);
	axis(side=1,seq(0,xmax,10), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,10),at=seq(0,xmax,10), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,length.out=4), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,length.out=4)/times,at=seq(0,ymax,length.out=4), side=2, las=3, line=1, col="black",cex=1);	
	arrows(sumobs,ymax*0.8,sumobs,0,lwd=3);
	text(sumobs,ymax*0.9,paste("   = ",length(which(sumsim>sumobs))/times*100,"%",sep=""),adj=0.1);
	text(sumobs,ymax*0.9,"P",adj=1,font=3);
	text(sumobs,ymax,"Observed",adj=0.15);
	text(xmax*-0.1,ymax*0.95,mutposi[i],adj=0);		
	mtext("Number of genes  (length > 200 aa)", side=1, las=1, line=2, col="black",cex=1);
	mtext("which have hotspot muations' sites", side=1, las=1, line=3, col="black",cex=1);
	mtext(picname[i],side=3, las=1, line=1, col="black",cex=1.5,adj=-0.2);	
	box();
}


picname<-c("c","d");
result2<-result[which(as.numeric(result[,2])>=200),];
data2<-data[which(data[,9]>=200),];
for(i in c(1:length(mutposi))){
	simdata<-result2[which(result2[,1]==mutposi[i]),];
	obsdata<-data2[which(data2[,3]==mutposi[i]),];
	sumobs<-sum(obsdata[,6]);
	sumsim<-apply(simdata[,c(3:dim(simdata)[2])],2,function(x){sum(as.numeric(x))});
	ymax<-2000;
	xmin<--10;
	xmax<-80;
	hist(sumsim,breaks=30,xlab="",ylab="Frequency",main=paste("median=",median(sumsim),sep=""),xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",cex.lab=1.2,freq=T);
	axis(side=1,seq(0,xmax,20), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,20),at=seq(0,xmax,20), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,length.out=5), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,length.out=5)/times,at=seq(0,ymax,length.out=5), side=2, las=3, line=1, col="black",cex=1);	
	arrows(sumobs,ymax*0.8,sumobs,0,lwd=3);
	text(sumobs,ymax*0.9,paste("   = ",length(which(sumsim>sumobs))/times*100,"%",sep=""),adj=0.1);
	text(sumobs,ymax*0.9,"P",adj=1,font=3);
	text(sumobs,ymax,"Observed",adj=0.15);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0);		
	mtext("Number of hotspot muations' sites", side=1, las=1, line=2, col="black",cex=1);
	mtext("in genes  (length > 200 aa)", side=1, las=1, line=3, col="black",cex=1);
	mtext(picname[i],side=3, las=1, line=1, col="black",cex=1.5,adj=-0.2);	
	box();
}

dev.off();

