更新20181216
今天又把KaKs_calculator2.0的帮助手册翻出来看了一遍,找到了一句话:Before calculation, gaps and stop codons between compared sequences will be removed.软件应该是自己会自动删除gaps和终止密码子,比对后应该不用自己手动删除终止密码子
参考 组学大讲堂 公众号分享的文章 如何使用 KaKs_Calculator 计算Ka/Ks 重复论文 Chloroplast genome analysis of resurrection tertiary relict Haberlea rhodopensis highlights genes important for desiccation stress response 中ccsA基因的ka/ks的计算过程:
论文中用到了五个物种ccsA基因,根据拉丁名搜索NCBI organelle genome数据库找到对应的CDS序列,使用 mafft 比对
mafft --auto --clustalout ccsA_KaKs_pra.fasta > ccsA_aligned.aln
(小疑问:mafft比对结果是默认把氨基酸大写字母转换成小写字母吗?)
KaKs_Calculator2.0 download | SourceForge.net 由此网址下载KaKs_calculator
解压出来以后分别用到src目录下的 AXTConvertor 程序和bin/Linux目录下的 KaKs_Calculator 程序,使用之前使用 chmod u+x AXTConvertor
chmod u+x KaKs_Calculator 赋予文件可执行权限
之后用 AXTConvertor 程序将比对结果文件转化为axt文件
./src/AXTConvertor ccsA_aligned.aln ccsA.axt
然后用 KaKs_Calculator 计算 Ka/Ks值
./bin/Linux/KaKs_Calculator -i ccsA.axt -o ccsA_2.kaks -m YN -c 11
-o参数指定输出文件, -m参数指定方法,-c参数指定密码子表;具体内容可以使用 -h 参数查看帮助文档
ccsA_fasta_and_kaks_result.rar_免费高速下载|百度网盘-分享无限制
分析所用到的fasta文件和最终输出结果的百度云链接
问题:输入5条序列会输出C52个Ka/Ks值,最终的结果是取平均值吗?
最后用ggplot2简单画一幅箱线图
pal2nal.pl脚本将蛋白比对转换成CDS比对命令
perl pal2nal.pl example_pep_aln.fasta example_CDS.fasta -output.fasa > example_CDS_aln.fasta
下载地址 http://www.bork.embl.de/pal2nal/#Download
更新:20181215
比对格式之间互相转化(使用biopython的AlignIO模块)
from Bio import AlignIO
from Bio import Alphabet
AlignIO.convert("inputfile","fasta","outfile","clustal",alphabet=Alphabet.generic_dna)
参考文献
微信文章: 《 一个计算同义突变(Ks)与非同义突变(Ka)的方法》(其中有1:先比对蛋白序列,然后利用perl脚本转换成CDS比对;2:将核酸序列联配转换成AXT格式的python脚本)