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

欢迎来到"bio生物信息"的世界

今天给大家带来“批量计算kaks值”的技能。

关于kaks的背景知识我就不介绍了,感兴趣的自行搜索,这里直接开始讲怎么批量计算kaks值。

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.fapep.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文件后,修改第一行和第三行,如下所示两张图,分别为修改前和修改后(请务必修改你有权限的路径):

image
image

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

生成的alig_pep.fa文件如下:


image

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文件如下所示:

image

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

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

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

拆分后,会生成gene1.fagene2.fagene3.fa三个文件。

接下来,将生成的gene1.fagene2.fagene3.fa组成新的基因对gene1-gene2gene1-gene3gene2-gene3,见如下命令:

for i in $(seq 1 3)

do

cat gene1.fa gene${i}.fa > gene1_${i}.fa

done

cat gene2.fa gene3.fa > gene2_3.fa

生成如下几个文件:

gene1_1.fa

gene1_2.fa

gene1_3.fa

gene2_3.fa

其中,gene1_2.fagene1_3.fagene2_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

下载后,输入如下命令:

for i in *.fa

do

echo $i

perl parseFastaIntoAXT.pl $i

done

生成三个文件:

gene1_2.fa.axt

gene1_3.fa.axt

gene2_3.fa.axt

6 计算kaks值

下载安装kakscalculator

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

输入以下命令:

for i in *.fa.axt

do

echo $i

KaKs_Calculator -i $i -o ${i%%.*}.kaks

done

生成三个文件:gene1_2.kaksgene1_3.kaksgene2_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"

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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,294评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,493评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,790评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,595评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,718评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,906评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,053评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,797评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,250评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,570评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,711评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,388评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,018评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,796评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,023评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,461评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,595评论 2 350

推荐阅读更多精彩内容