研究密碼子偏向性的論文通常都會(huì)分析RSCU值，論文中通常會(huì)用堆積柱形圖來展示RSCU的值，之前在論文里也看到過下面這幅圖的形式展示RSCU分析的結(jié)果
之前也有人在公眾號(hào)留言問過這幅圖改如何實(shí)現(xiàn)，但是自己當(dāng)時(shí)也不知道。今天看Y叔的公眾號(hào)文章aplot包：讓你畫出更復(fù)雜的圖，解決的主要問題是拼圖的時(shí)候坐標(biāo)軸對(duì)齊。
看過這篇文章后突然想到展示密碼子偏向性RSCU值的這幅圖可以借助拼圖來實(shí)現(xiàn)，下面介紹自己的實(shí)現(xiàn)過程。
首先是計(jì)算RSCU值
我借助python中的CAI模塊實(shí)現(xiàn)
https://github.com/Benjamin-Lee/CodonAdaptationIndex
使用pip直接安裝
pip install CAI
計(jì)算RSCU值
from CAI import RSCU
from Bio import SeqIO
seqs = [rec.seq for rec in SeqIO.parse('codon_usage_example.fasta','fasta')]
rscu = RSCU(seqs)
rscu是一個(gè)字典，密碼子是鍵，對(duì)應(yīng)的RSCU是值
寫一個(gè)簡(jiǎn)單的腳本獲得使用R語言的ggplot2作圖的輸入文件
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()
獲得的rscu.txt四列
- 第一列是密碼子
- 第二列是對(duì)應(yīng)的氨基酸
- 第三列是RSCU值
- 第四列是數(shù)字，用來填充顏色和控制位置
接下來是ggplot2作圖代碼
install.packages('aplot')
library(ggplot2)
library(aplot)
df<-read.csv('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))
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)
這樣就得到了開頭提到的圖
還發(fā)現(xiàn)了一個(gè)R包可以分析密碼子偏向性
sscu，具體用法沒看，用到再說
還發(fā)現(xiàn)了一個(gè)網(wǎng)站分析密碼子偏向性
http://www.codons.org/Help.html#CU
還看到了一個(gè)python模塊可以把對(duì)應(yīng)的蛋白質(zhì)序列弄回核苷酸序列，不知道什么情況下會(huì)用到
https://pypi.org/project/codon-harmony/
歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本
需要示例數(shù)據(jù)和代碼的可以關(guān)注公眾號(hào)留言哈！
也歡迎大家到B站 https://space.bilibili.com/355787260
關(guān)注我，有時(shí)間會(huì)把這些文字內(nèi)容錄制成視頻發(fā)布到b站的
