(12) 塩基置換

#(U,C,A,G)=(1,2,3,4)
b<-1:4
#b<-c("T","C","A","G")
codons<-as.matrix(expand.grid(b,b,b))
aas<-c(
"F","F","L","L","S","S","S","S","Y","Y","X","X","C","C","X","W","L","L","L","L","P","P","P","P","H","H","Q","Q","R","R","R","R","I","I","I","M","T","T","T","T","N","N","K","K","S","S","R","R","V","V","V","V","A","A","A","A","D","D","E","E","G","G","G","G"
)


nonStops<-which(aas!="X")
st<-which(aas=="M")
Stops<-which(aas=="X")

L<-10
peptides<-c(st,sample(nonStops,L-1,replace=TRUE))

bases<-c()
for(i in 1:L){
	tmpaa<-peptides[i]
	candidates<-which(aas==aas[tmpaa])
	tmpcodon<-sample(candidates,1)
	bases<-c(bases,as.vector(codons[tmpcodon,]))
}
bases

base2aa<-function(b,cs){
	len<-length(b)/3
	ret<-c()
	for(i in 1:len){
		triplet<-b[(i*3):((i-1)*3+1)]
		for(j in 1:length(cs[,1])){
			if(prod(triplet==cs[j,])){
				ret<-c(ret,aas[j])
			}
		}
	}
	ret
}
peptides2<-base2aa(bases,codons)

Niter<-1000

M<-matrix(0,Niter+1,length(bases))
M[1,]<-bases

rt<-0.1
for(i in 1:Niter){
	loop<-TRUE
	while(loop){
		rs<-runif(length(bases))
		newbs<-bases
		muts<-which(rs<rt)
		#print(muts)
		if(length(muts)>0){
			tmpb<-sample(b,length(muts),replace=TRUE)
			newbs[muts]<-tmpb
		}
		#print(bases)
		#print(newbs)
		newpeptides<-base2aa(newbs,codons)
		#print(newpeptides)
		#print(peptides)
		if(identical(newpeptides,peptides2)){
			loop<-FALSE
			M[i+1,]<-newbs
		}
		
	}
	
}


image(M)

M2<-t(t(M)-bases%%4)
image(t(sign(M2%%4)))
grid(nx=length(M2[1,]),ny=1)