KaKs_Calculator2.0+ParaAT计算两个物种间

1.软件安装

⭐KaKs_Calculator2.0
⭐ParaAT
KaKs_Calculator2.0

cd KaKs_Calculator2.0/bin/Linux
chmod 744 KaKs_Calculator

echo 'PATH=$PATH:/KaKs_Calculator2.0/bin/Linux' >> ~/.bashrc
source ~/.bashrc

ParaAT:http://download.big.ac.cn/bigd/tools/ParaAT2.0.tar.gz

#解压即可使用
tar xf ParaAT2.0.tar.gz

2.准备两个物种的cds及pep文件

去掉序列名后面的orf注释

#!/usr/bin/env python3
import sys
x = sys.argv[1]
file = open(x, "r")
lines = file.readlines()
for line in lines:
    if ">" in line:
       strings = line.split(" ")[0]
       print(strings)
    else:
       print(line.rstrip())

###保存为change_name.py
python change_name.py 1.cds > 1_new.cds

3.筛选直系同源基因

#建库
 makeblastdb -in ref.pep.fa -dbtype prot
#blast比对
 blastp -query pep.fa -db ref.pep.fa -evalue 1e-5 -max_target_seqs 1 -num_threads 2 -out 1_2_blastp_out.m6 -outfmt 6
#取双向最优配对序列的序列名
 cut -f1-2 1_2_blastp_out.m6 | sort | uniq > 1_2.homolog

4.ParaAT将蛋白序列比对结果转化为cds序列比对结果

此步需注意，file.homolog、file.cds、file.pep三个文件的基因ID需保持一致，且file.cds及file.pep为fasta格式，即“>”后面只接基因ID号，否则会报错

eg：

file.cds	file.homolog

此处cds文件比homolog文件的基因ID多了“_mRNA”，所以要删掉

#合并两个物种的cds及pep序列
cat cds.fa ref.cds.fa >total.cds.fa
cat pep.fa ref.pep.fa >total.pep.fa
#将clustalW和ParaAT写进.bashrc
export PATH=/path/clustalw-2.1/bin/:$PATH
export PATH=/path/ParaAT2.0/:$PATH
# 新建proc文件, 设定12个线程
echo "12" >proc
#运行ParaAT，此步会在环境中搜索KaKs_Calculator的路径
ParaAT.pl -h 1_2.homolog -n total.cds.fa -a total.pep.fa -p proc -m clustalw2 -f axt -g -k -o yeast
-h, 同源基因名称文件
-n, 指定核酸序列文件
-a, 指定蛋白序列文件
-p, 指定多线程文件
-m, 指定比对工具
-g, 去除比对有gap的密码子
-k, 用KaKs_Calculator 计算kaks值
-o, 输出结果的目录
-f, 输出比对文件的格式
*** 也可通过-f参数得到其他软件分析ka/k所需的格式

得到每一对同源基因的Ka/Ks值

*.kaks文件中包括：
● Sequence: Name of Pairwise sequence
● Method: Name of method for calculation of Ka and Ks
● Ka: Nonsynonymous substitution rate
● Ks: Synonymous substitution rate
● Ka/Ks: Selective strength
● P-Value (Fisher): The value computed by Fisher exact test
● Length: Sequence length (after removing gaps and stop codon(s))
● S-Sites: Synonymous sites
● N-Sites: Nonsynonymous sites
● Fold-Sites (0:2:4): 0,2,4-fold degenerate sites
● Substitutions: Substitutions between sequences
● S-Substitutions: Synonymous substitutions
● N-Substitutions: Nonsynonymous substitutions
● Fold-S-Substitutions (0:2:4): Synonymous substitutions at 0,2,4-fold
● Fold-N-Substitutions (0:2:4): Nonsynonymous substitutions at 0,2,4-fold
● Divergence-Time: Divergence time
● Substitution-Rate-Ratio (rTC:rAG:rTA:rCG:rTG:rCA/rCA): Ratios of six substitution
rates to the substitution rate between C and A
● GC(1:2:3): GC content of entire sequences and of three codon positions
● ML-Score: Maximum likelihood score
● AICc: Value of AICc
● Akaike-Weight: Value of Akaike weight for model selection
● Model: Selected model for the method of MS
9
合并所有同源基因对的Ka/Ks值，提取前5列，并排序去重

for i in `ls *.kaks`;do awk 'NR>1{print $1"\t"$3"\t"$4"\t"$5}' $i >>all-kaks.txt;done
#按科学计数法对第四列降序排列
sort -k 4 -g -r all-kaks.txt|uniq >all-kaks.results
#加标题
sed -i '1i\Seq\tKa\tKs\tKa/Ks' all-kaks.results

参考文章：
【比较基因组】如何利用paml计算kaks - 知乎 (zhihu.com)
WGD（全基因组复制）分析——Ka/Ks及4Dtv值计算 - 简书 (jianshu.com)
Kaks_calculator计算ka/ks 值 - 简书 (jianshu.com)
https://www.sciencedirect.com/science/article/pii/S1672022910600083?via%3Dihub