set.seed(1);
data<-read.table(file="99.01hotmut.txt",sep="\t",head=T);
sum(data[,"avisite"])/334
times<-10000;
result<-NULL;
for(k in seq(0,1,1)){
	for(j in c(1:dim(data)[1])){
		tem4<-NULL;
		for(i in c(1:times)){
			tem<-sample(c(1:(data[j,"alainlen"]-(data[j,"alainlen"]-data[j,"avisite"])*k)),data[j,"mutnum"],replace=T);
			tem2<-sapply(split(tem,tem),length);
			tem3<-length(which(tem2>1));
			tem4<-c(tem4,tem3);
		}
		result<-rbind(result,c(k,as.character(data[j,"leadorlag"]),data[j,"genelen"],tem4));
	}
}
write.table(result,"99.02simulation.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=T);
times<-10000;
this1<-function(k,picname,tbreaks,ymax,ygap,xmin,xmax,xgap,xname1,xname2,xfont,xname3,xname4,xfont2){
	leadobs<-sum(data[which(data[,"leadorlag"]=="leading" & data[,"hotnum"]>0),"hotnum"]);
	laggobs<-sum(data[which(data[,"leadorlag"]=="lagging" & data[,"hotnum"]>0),"hotnum"]);
	ratobs<-laggobs/leadobs;
	leadsim<-result[which(result[,1]==k & result[,2]=="leading"),];
	laggsim<-result[which(result[,1]==k & result[,2]=="lagging"),];
	leadsum<-apply(leadsim[,c(4:dim(leadsim)[2])],2,function(x){sum(as.numeric(x))});
	laggsum<-apply(laggsim[,c(4:dim(laggsim)[2])],2,function(x){sum(as.numeric(x))});
	ratsim<-laggsum/leadsum;
	pvalue<-length(which(ratsim>ratobs))/times*2;
	if(pvalue>1){
		pvalue<-2-pvalue;
	}
	hist(ratsim,breaks=tbreaks,xlab="",ylab="",main="",xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",freq=T);
	axis(side=1,seq(xmin,xmax,xgap), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(xmin,xmax,xgap),at=seq(xmin,xmax,xgap), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,ygap), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,ygap)/times,at=seq(0,ymax,ygap), side=2, las=3, line=1, col="black",cex=1);	
	arrows(ratobs,ymax*0.8,ratobs,ymax*0.6,lwd=3,col="black");
	text(xmin+(xmax-xmin)*0.85,ymax*0.95,paste("   = ",sprintf("%0.2f",pvalue),sep=""),adj=0.1,cex=1.2);
	text(xmin+(xmax-xmin)*0.85,ymax*0.95,"P",adj=1,font=3,cex=1.2);
	text(ratobs,ymax*0.9,"Observed",cex=1.2);		
	mtext("Frequency", side=2, las=3, line=3, col="black",cex=1.2);	
	mtext(xname1, side=1, las=1, line=2, col="black",cex=1.2,font=xfont);
	mtext(xname2, side=1, las=1, line=3, col="black",cex=1.2,font=xfont);
	mtext(xname4, side=1, las=1, line=3, col="black",cex=1.2,font=xfont2);
	mtext(xname3, side=1, las=1, line=3.5, col="black",cex=0.8);
	mtext(picname,side=3, las=1, line=1, col="black",cex=2.5,adj=-0.2,font=2);	
	box();
}
this2<-function(k,picname,tbreaks,ymax,ygap,xmin,xmax,xgap,xname1,xname2,xfont,xname3,xname4,xfont2){
	leadobs<-length(which(data[,"leadorlag"]=="leading" & data[,"hotnum"]>0));
	laggobs<-length(which(data[,"leadorlag"]=="lagging" & data[,"hotnum"]>0));
	ratobs<-laggobs/leadobs;
	leadsim<-result[which(result[,1]==k & result[,2]=="leading"),];
	laggsim<-result[which(result[,1]==k & result[,2]=="lagging"),];
	leadsum<-apply(leadsim[,c(4:dim(leadsim)[2])],2,function(x){length(which(x>0))});
	laggsum<-apply(laggsim[,c(4:dim(laggsim)[2])],2,function(x){length(which(x>0))});
	ratsim<-laggsum/leadsum;
	pvalue<-length(which(ratsim>ratobs))/times*2;
	if(pvalue>1){
		pvalue<-2-pvalue;
	}
	hist(ratsim,breaks=tbreaks,xlab="",ylab="",main="",xlim=c(xmin,xmax),ylim=c(0,ymax),xaxt="n",yaxt="n",freq=T);
	axis(side=1,seq(xmin,xmax,xgap), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(xmin,xmax,xgap),at=seq(xmin,xmax,xgap), side=1, las=1, line=0.7, col="black",cex=1);
	axis(side=2,seq(0,ymax,ygap), tcl=-0.5, line=0, labels=FALSE);
	mtext(seq(0,ymax,ygap)/times,at=seq(0,ymax,ygap), side=2, las=3, line=1, col="black",cex=1);	
	arrows(ratobs,ymax*0.8,ratobs,ymax*0.6,lwd=3,col="black");
	text(xmin+(xmax-xmin)*0.85,ymax*0.95,paste("   = ",sprintf("%0.2f",pvalue),sep=""),adj=0.1,cex=1.2);
	text(xmin+(xmax-xmin)*0.85,ymax*0.95,"P",adj=1,font=3,cex=1.2);
	text(ratobs,ymax*0.9,"Observed",cex=1.2);	
	mtext("Frequency", side=2, las=3, line=3, col="black",cex=1.2);	
	mtext(xname1, side=1, las=1, line=2, col="black",cex=1.2,font=xfont);
	mtext(xname2, side=1, las=1, line=3, col="black",cex=1.2,font=xfont);
	mtext(xname4, side=1, las=1, line=3, col="black",cex=1.2,font=xfont2);
	mtext(xname3, side=1, las=1, line=3.5, col="black",cex=0.8);
	mtext(picname,side=3, las=1, line=1, col="black",cex=2.5,adj=-0.2,font=2);	
	box();

}

tiff(filename = paste("99.02simulation.tif",sep=""),width=250, height=200, units="mm", res=300);
layout(matrix(c(1:4),2,2,byrow=T));
par(oma=c(0,0,0,0),mar=c(5,5,4,3));
this1(1,"a",20,times*0.3,times*0.1,0.1,0.35,0.05,"Ratio of # of multi-hit sites","between lagging and leading strands (    )",1,
"                                                                                             N",
"                                                        R",3);
this2(1,"b",20,times*0.3,times*0.1,0.1,0.35,0.05,"Ratio of # of multi-hit-site-containing genes","between lagging and leading strands (    )",1,
"                                                                                             G",
"                                                        R",3);
this1(0,"c",20,times*0.4,times*0.1,0,0.8,0.2,"","R",3,"       N","",1);
this2(0,"d",20,times*0.4,times*0.1,0,0.8,0.2,"","R",3,"       G","",1);
dev.off();


