kaks calculator批量计算多个基因的选择压力kaks值

文章来源:http://www.cnblogs.com/chenwenyan/

1 文件准备

首先准备两个文件,一个是基因的cds序列,一个是蛋白质序列。

cds序列和蛋白质可以在ensembl网站找到:http://ftp.ensembl.org/pub/current_fasta/

这两个文件的示例如下:

cds序列文件cds.fa

>gene1

ATGGAGGTTGCAATGGTGAGTGCGGAGAGCTCAGGATGCAACAGTCACATGCCTTACGGT

TATGCTGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTTGCTCACTCCAGGGCAGCTGCG

GCAGCTGCCGTTGCAGCGGCCACAGCTGCCGTCGAAGGAAGTGGGGGTTCTGGTGGGGGG

>gene2

ATGGAGGTTGCAATGGTGAGTGCGGAGAGCTCAGGGTGCAACAGTCACATGCCTTATGGT

TATGCTGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTTGCTCACTCCAGGGCAGCTGCA

GCAGCTGCTGTTGCAGCGGCCACAGCTGCTGTCGAAGGTAGCGGGGGTTCTGGTGGGGGC

TCCCAC

>gene3

ATGGAGGTGGCGATGGTGAGTGCGGAGAGCTCAGGGTGCAACAGTCACATGCCTTACGGG

TACGCGGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTGGCTCACTCCCGGGCGGCGGCG

GCCGCCGCCGTCGCGGCTGCCACGGCTGCCGTGGAAGGCAGTGGGGGGCCTGG

蛋白质序列pep.fa

>gene1

MEVAMVSAESSGCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAATAAVEGSGGSGGG

>gene2

MEVAMVSAESSRCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAATAAVEGSGSSGGGSH

>gene3

MEVAMVSAESSGCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAAKAAVEGSGGP

注意:cds.fa和pep.fa文件的序列ID号(gene1,2,3)要一致。

2 对蛋白质序列pep.fa进行比对

进行蛋白质序列比对前,需要先安装mafft软件。

下载mafft软件:

wget https://mafft.cbrc.jp/alignment/software/mafft-7.453-with-extensions-src.tgz

tar -zxvf mafft-7.453-with-extensions-src.tgz

cd mafft-7.453-with-extensions/core

安装:

1)有root权限用户安装法:

make clean

make

su

make install

2)无root权限用户安装法:

vi Makefile

进入makefile文件后,修改第一行和第三行,如下所示两张图,分别为修改前和修改后(请务必修改你有权限的路径):


安装成功后,输入命令mafft --auto pep.fa > alig_pep.fa

生成的alig_pep.fa文件如下:


3 将比对好的蛋白质序列与cds序列比对

这一步需要下载pal2nal.pl文件

wget http://www.bork.embl.de/pal2nal/distribution/pal2nal.v14.tar.gz

tar -zxvf pal2nal.v14.tar.gz

cd pal2nal.v14/

下载后就能看见pal2nal.pl这个文件

随后将蛋白质序列与cds序列比对

pal2nal.pl alig_pep.fa cds.fa -output fasta > cds_pep_aln.fa

比对好的cds_pep_aln.fa文件如下所示:


4 生成计算kaks值时的基因对

计算kaks值前需要先将cds_pep_aln.fa文件拆分:

csplit cds_pep_aln.fa /\>/ -n2 -s {*} -f gene -b "%1d.fa" ; rm gene0.fa

拆分后,会生成gene1.fa 、gene2.fa、 gene3.fa三个文件。

接下来,将生成的gene1.fa、 gene2.fa、 gene3.fa组成新的基因对gene1-gene2、gene1-gene3、gene2-gene3,见如下命令:

fori in $(seq13)docatgene1.fa gene${i}.fa > gene1_${i}.fadonecatgene2.fa gene3.fa > gene2_3.fa

生成如下几个文件:

gene1_1.fa

gene1_2.fa

gene1_3.fa

gene2_3.fa

其中,gene1_2.fa、gene1_3.fa、gene2_3.fa为我们所需的基因对。

这里将他们提取成基因对,而不是多条序列进行计算的原因是,KaKs_Calculator这个软件在处理多序列的输入文件时,会报错:Error. The size of two sequences in 'gene1&gene2' is not equal。我尝试了很多遍,发现只能提取成基因对才不会报这种错误。当然,偶尔运气好的时候,KaKs_Calculator是能处理多序列的kaks值的,为了防止出错,我们还是将他们拆开计算好一点。

5 将gene1_2.fa、gene1_3.fa、gene2_3.fa文件转化为axt文件

转化为axt文件需要下载parseFastaIntoAXT.pl文件,文件地址:https://gitee.com/liaochenlanruo/kaks_pupline/blob/master/parseFastaIntoAXT.pl

下载后,输入如下命令:

foriin*.fadoecho$iperl parseFastaIntoAXT.pl$idone

生成三个文件:

gene1_2.fa.axt

gene1_3.fa.axt

gene2_3.fa.axt

6 计算kaks值

下载安装kakscalculator

下载链接:https://sourceforge.net/projects/kakscalculator2/

输入以下命令:

foriin*.fa.axtdoecho$iKaKs_Calculator -i$i-o${i%%.*}.kaksdone

生成三个文件:gene1_2.kaks、gene1_3.kaks、gene2_3.kaks

到这一步,批量计算kaks值的工作就已经搞定。

这里附上安装kaks_calculator软件可能会遇到报错:

g++ -O -o AXTConvertor AXTConvertor.cpp -lstdc++ -lm

AXTConvertor.cpp: In function ?.ool readPhylip()?.

AXTConvertor.cpp:210:22: error: ?.toi?.was not declared in this scope

if (atoi(num.c_str())!=sequence.size()){

AXTConvertor.cpp: In function ?.ool readNexus()?.

AXTConvertor.cpp:338:39: error: ?.toi?.was not declared in this scope

if (sequence.size()!=atoi(num.c_str())) >{

make: *** [AXTConvertor] Error 1

解决方法在这里:

编辑KaKs.cpp文件,加上include "string.h"

编辑AXTConvertor.cpp文件,加上include "stdlib.h"

编辑GY94.cpp文件,加上 include "string.h"

如无报错请忽略上述内容。

©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容