5/30 MIKUセミナー

の塩基配列
21個のアミノ酸
4*4*4*21の分割表
を作るのに、次元をどうこうしたいという話。
コドンを一生懸命打ったのでメモ。

#適当にmRNA配列を作成する
mRNA <- function(n=100, prob = rep(1, 4)){
	bases <- c("A", "T", "G", "C")
	res <- sample(bases, size=n, prob=prob, replace=TRUE)
	res1 <- NULL
	for(i in 1:length(res)){
		res1 <- paste(res1, a[i], sep="")
	}
	return(res1)
}
mRNA <- function(n=100, prob = rep(1, 4)){
	bases <- c("A", "T", "G", "C")
	res <- sample(bases, size=n, prob=prob, replace=TRUE)
	return(res)
}
w <- 100
mRNAseq <- lapply(mapply(rep, 0, w), unique)
for(e in 1:w){
	mRNAseq[[e]] <- mRNA(n = sample(seq(99, 300, by=3), size=1))
}

amino <- c("A", "C", "D", "E", "F", "G", "H", "I", "K", "L", "M", "N", "P", "Q", "R", "S", "T", "V", "W", "Y")
end <- c("T", "C", "A", "G")
codon <- list(paste("G", "C", end, sep=""),
			paste("T", "G", end[1:2], sep=""),
			paste("G", "A", end[1:2], sep=""),
			paste("G", "A", end[3:4], sep=""),
			paste("T", "T", end[1:2], sep=""),
			paste("G", "G", end, sep=""),
			paste("C", "A", end[1:2], sep=""),
			paste("A", "T", end[1:3], sep=""),
			paste("A", "A", end[3:4], sep=""),
			paste("C", "T", end, sep=""),
			paste("A", "T", end[4], sep=""),
			paste("A", "A", end[1:2], sep=""),
			paste("C", "C", end, sep=""),
			paste("C", "A", end[3:4], sep=""),
			paste("C", "G", end, sep=""),
			paste("T", "C", end, sep=""),
			paste("A", "C", end, sep=""),
			paste("G", "T", end, sep=""),
			paste("T", "G", end[4], sep=""),
			paste("T", "A", end[1:2], sep=""))
			
			#c(paste("T", "A", end[3:4], sep=""),paste("T", "G", "G", sep="")))
names(codon) <- amino
library(MCMCpack)
n_codon <- sapply(codon, length)
codon_prob <-lapply(mapply(rep, 0, length(codon)), unique) 
for(i in 1:length(n_codon)){
	codon_prob[[i]] <- rdirichlet(1, rep(1, n_codon[i]))
}

#アミノ酸の使われ方も変える
amino_prob <- rdirichlet(1, rep(1, length(amino)))

seqs <- 100
result <- lapply(mapply(rep, 0, seqs), unique)
for(resseq in 1:seqs){
	use_amino <- sample(30:300, size=1)
	temp_res <- rep(0, use_amino)
	for(ua in 1:use_amino){
		temp_use <- sample(1:length(amino), size=1)
		temp_codon <- sample(codon[[temp_use]], prob=codon_prob[[temp_use]], size=1)
		temp_res[ua] <- temp_codon
	}
	result[[resseq]] <- temp_res
}

　
20120610追記
簡潔に書くスクリプトを教えていただいた。

f1 <- c("G", "T", "G", "G", "T", "G", "C", "A", "A", "C", "A", "A", "C", "C", "C", "T", "A", "G", "T", "T")
s1 <- c("C", "G", "A", "A", "T", "G", "A", "T", "A", "T", "T", "A", "C", "A", "G", "C", "C", "T", "G", "A")
f2 <- c(1,1,1,3,1,1,1,1,3,1,4,1,1,3,1,1,1,1,4,1)
s2 <- c(4,2,2,4,2,4,2,3,4,4,4,2,4,4,4,4,4,4,4,2)
codon <- unname(mapply(function(x, y, i, j) paste(x, y, end[i:j], sep=""), f1, s1, f2, s2))