1.计算RSCU值
可以直接用pip安装CAI
pip install CAI
自己写个RSCU.py计算RSCU值,当然你也可以通过在线版EMBOSS来计算密码子偏性
from CAI import RSCU
from Bio import SeqIO
c2aa = {
'TGT':'Cys',
'UGU':'Cys',
'TGC':'Cys',
'UGC':'Cys',
'GAT':'Asp',
'GAU':'Asp',
'GAC':'Asp',
'GAC':'Asp',
'TCT':'Ser',
'UCU':'Ser',
'TCG':'Ser',
'UCG':'Ser',
'TCA':'Ser',
'UCA':'Ser',
'TCC':'Ser',
'UCC':'Ser',
'AGC':'Ser',
'AGC':'Ser',
'AGT':'Ser',
'AGU':'Ser',
'CAA':'Gln',
'CAA':'Gln',
'CAG':'Gln',
'CAG':'Gln',
'ATG':'Met',
'AUG':'Met',
'AAC':'Asn',
'AAC':'Asn',
'AAT':'Asn',
'AAU':'Asn',
'CCT':'Pro',
'CCU':'Pro',
'CCG':'Pro',
'CCG':'Pro',
'CCA':'Pro',
'CCA':'Pro',
'CCC':'Pro',
'CCC':'Pro',
'AAG':'Lys',
'AAG':'Lys',
'AAA':'Lys',
'AAA':'Lys',
'ACC':'Thr',
'ACC':'Thr',
'ACA':'Thr',
'ACA':'Thr',
'ACG':'Thr',
'ACG':'Thr',
'ACT':'Thr',
'ACU':'Thr',
'TTT':'Phe',
'UUU':'Phe',
'TTC':'Phe',
'UUC':'Phe',
'GCA':'Ala',
'GCA':'Ala',
'GCC':'Ala',
'GCC':'Ala',
'GCG':'Ala',
'GCG':'Ala',
'GCT':'Ala',
'GCU':'Ala',
'GGT':'Gly',
'GGU':'Gly',
'GGG':'Gly',
'GGG':'Gly',
'GGA':'Gly',
'GGA':'Gly',
'GGC':'Gly',
'GGC':'Gly',
'ATC':'Ile',
'AUC':'Ile',
'ATA':'Ile',
'AUA':'Ile',
'ATT':'Ile',
'AUU':'Ile',
'TTA':'Leu',
'UUA':'Leu',
'TTG':'Leu',
'UUG':'Leu',
'CTC':'Leu',
'CUC':'Leu',
'CTT':'Leu',
'CUU':'Leu',
'CTG':'Leu',
'CUG':'Leu',
'CTA':'Leu',
'CUA':'Leu',
'CAT':'HIS',
'CAU':'HIS',
'CAC':'HIS',
'CAC':'HIS',
'CGA':'Arg',
'CGA':'Arg',
'CGC':'Arg',
'CGC':'Arg',
'CGG':'Arg',
'CGG':'Arg',
'CGT':'Arg',
'CGU':'Arg',
'AGG':'Arg',
'AGG':'Arg',
'AGA':'Arg',
'AGA':'Arg',
'TGG':'Trp',
'UGG':'Trp',
'GTA':'Val',
'GUA':'Val',
'GTC':'Val',
'GUC':'Val',
'GTG':'Val',
'GUG':'Val',
'GTT':'Val',
'GUU':'Val',
'GAG':'Glu',
'GAG':'Glu',
'GAA':'Glu',
'GAA':'Glu',
'TAT':'Tyr',
'UAU':'Tyr',
'TAC':'Tyr',
'UAC':'Tyr',
}
seqs = [rec.seq for rec in SeqIO.parse('codon_usage_example.fasta','fasta')]
rscu = RSCU(seqs)
fw = open('rscu.txt','w')
amino_acid = {}
for aa,bb in rscu.items():
if c2aa[aa] not in amino_acid:
amino_acid[c2aa[aa]] = 6
else:
amino_acid[c2aa[aa]] -= 0.5
print(aa,c2aa[aa],round(bb,3),amino_acid[c2aa[aa]])
fw.write(aa+","+c2aa[aa]+","+str(round(bb,3))+","+str(amino_acid[c2aa[aa]])+"\n")
fw.close()
输出结果为下列所示的一个4列文件:
2.R绘堆积图
setwd('C:\\Users\\liu\\Desktop\\1118\\RSCU')
install.packages('aplot')
library(ggplot2)
library(aplot)
df<-read.csv('Flag.rscu.txt',header=F,stringsAsFactors = F)
p1<-ggplot(df,aes(fill=as.character(V4),x=V2,y=V3))+
geom_bar(position = "stack",stat="identity")+
theme_bw()+scale_y_continuous(expand=c(0,0),
limits = c(0,6.2))+
theme(legend.position = "none")+labs(y="RSCU",x="")
#geom_text(aes(label=V1),position=position_stack(vjust=0.5))
pdf('Flag.pdf')
p1
p2<-ggplot(df,aes(x=V2,y=V4))+
geom_label(aes(label=V1,fill=as.character(V4)),
size=2)+
labs(x="",y="")+ylim(3.4,6.1)+
theme_minimal()+
theme(legend.position = "none",
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
p1%>%
insert_bottom(p2,height = 0.3)
dev.off()