set.seed(1);
data<-read.table(file="99.02hotmut.txt",sep="\t",head=F);#name name leading/lagging 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[,4]),tem5);
write.table(result,"99.03simulation.txt",append=F,quote=F,row.names=F,col.names=F,sep="\t");

tiff(filename = paste("99.03simulation.tif",sep=""),width=180, height=180, units="mm", res=300);
layout(matrix(c(1:4),2,2,byrow = T));
picname<-c("a","b");
mutposi<-c("leading","lagging");
for(i in c(1:length(mutposi))){
	simdata<-result[which(result[,1]==mutposi[i]),];#
	obsdata<-data[which(data[,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<-50;
	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=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.9);
	text(sumobs,ymax*0.9,"P                              ",font=3);#
	text(sumobs,ymax,"Observed",adj=0.9);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0.2);		
	mtext("Number of genes", 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");
mutposi<-c("leading","lagging");
for(i in c(1:length(mutposi))){
	simdata<-result[which(result[,1]==mutposi[i]),];#
	obsdata<-data[which(data[,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,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=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.9);
	text(sumobs,ymax*0.9,"P                          ",font=3);#
	text(sumobs,ymax,"Observed",adj=0.9);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0.2);		
	mtext("Number of ", side=1, las=1, line=2, col="black",cex=1);
	mtext("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();
}
dev.off();
###



set.seed(1);
data<-read.table(file="99.02hotmut.txt",sep="\t",head=F);#CAB12455 guaA leading 511  9  0 7
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[,4]),tem5);
write.table(result,"99.03simulation_extrem.txt",append=F,quote=F,row.names=F,col.names=F,sep="\t");

tiff(filename = paste("99.03simulation_extrem.tif",sep=""),width=180, height=180, units="mm", res=300);
layout(matrix(c(1:4),2,2,byrow = T));
picname<-c("a","b");
mutposi<-c("leading","lagging");
for(i in c(1:length(mutposi))){
	simdata<-result[which(result[,1]==mutposi[i]),];#
	obsdata<-data[which(data[,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<-0;
	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.9);
	text(sumobs,ymax*0.9,"P                              ",font=3);#
	text(sumobs,ymax,"Observed",adj=0.9);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0.2);		
	mtext("Number of genes", 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");
mutposi<-c("leading","lagging");
for(i in c(1:length(mutposi))){
	simdata<-result[which(result[,1]==mutposi[i]),];#
	obsdata<-data[which(data[,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<-0;
	xmax<-1000;
	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,200), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,xmax,200),at=seq(0,xmax,200), 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.9);
	text(sumobs,ymax*0.9,"P                          ",font=3);#
	text(sumobs,ymax,"Observed",adj=0.9);
	text(xmax*0.8,ymax*0.95,mutposi[i],adj=0.2);		
	mtext("Number of ", side=1, las=1, line=2, col="black",cex=1);
	mtext("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();
}
dev.off();


