密码子偏向性分析的论文中通常会计算每个蛋白编码基因的总的GC含量和GC1/GC2/GC3的含量。找到一个现成的工具是emboss工具集中的cusp命令，在线版链接是 http://www.bioinformatics.nl/cgi-bin/emboss/cusp
但是这个工具好像只能一次计算一条序列的GC含量，如果一次性上传多条蛋白编码基因序列，计算结果是所有序列加到一起的值，而不是每条序列一个值。
试着自己写了一个python脚本
import sys
from Bio import SeqIO
from Bio.SeqUtils import Seq
from Bio.SeqUtils import GC
in_fasta = sys.argv[1]
out_results = sys.argv[2]
def calculate_GC123(seq):
GC1 = []
GC2 = []
GC3 = []
for i in range(1,len(seq)+1):
if i%3 == 1:
GC1.append(str(seq[i-1]))
elif i%3 == 2:
GC2.append(str(seq[i-1]))
else:
GC3.append(str(seq[i-1]))
GC1_seq = Seq(''.join(GC1))
GC2_seq = Seq(''.join(GC2))
GC3_seq = Seq(''.join(GC3))
return [round(GC(GC1_seq),2),round(GC(GC2_seq),2),round(GC(GC3_seq),2)]
fw = open(out_results,'w')
fw.write("%s\t%s\t%s\t%s\t%s\n"%("gene_name","GC_total","GC_first","GC_second","GC_third"))
for rec in SeqIO.parse(in_fasta,'fasta'):
GC_total = round(GC(rec.seq),2)
GC123 = calculate_GC123(rec.seq)
GC_first = GC123[0]
GC_second = GC123[1]
GC_third = GC123[2]
fw.write("%s\t%s\t%s\t%s\t%s\n"%(rec.id,str(GC_total),str(GC_first),str(GC_second),str(GC_third)))
fw.close()
print("Well done!")
使用需要安装biopython
使用方式
python calculate_GC_at_1_2_3_codon_position.py cds.fasta output.txt
第一个输入文件是自己的cds序列，fasta格式。第二个是输出文件，自己指定名字。
另外使用codonW时遇到的一个问题 使用如下命令运行codonw
codonW/codonw cds.fasta -all_indices -nomenu
得到一个结果文件cds.out
这里面有一些指标
比如T3s C3s GC3s是什么意思？我最开始以为是密码子第三位的T含量和C含量和GC含量，但是这个数值和 cusp计算的结果不一样。
大家有人知道这个指标是什么意思吗？
欢迎大家关注我的公众号
小明的数据分析笔记本
