使用paml计算kaks值练习

感谢大壮大佬的教程，教程来源 https://zhuanlan.zhihu.com/p/105159767

软件安装

软件链接：http://abacus.gene.ucl.ac.uk/software/paml.html
安装方法很简单，win的直接下win，linux的下载win的之后，按下面操作即可安装

wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgz
tar xf paml4.9j.tgz
cd paml4.9j
rm bin/*.exe 
cd src 
make -f Makefile 
rm *.o 
cp baseml basemlg codeml pamp evolver yn00 chi2 ../bin 

1、数据准备

初步组装的转录组数据去冗余后使用trandecoder进行ORF预测，会获得cds和pep两个文件，pep为预测到的氨基酸序列，cds为相对应的核苷酸序列。
文件中的ID序列会有注释信息，如下图

image.png

可以先将准备的cds和pep文件的注释信息去掉
改名的change_name.py脚本内容如下

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())

脚本运行方式为

python change_name.py filename > newfilename

我准备了016和017两个转录组样本的cds和pep文件，更改后的ID如下

image.png

准备好pep和cds文件之后，开始进行同源直系基因的筛选，标准按最宽松的来，evalue设置阈值1e-5，然后取双向最优比对序列即可（如果是发文章的话，请最好按前人的标准去筛选）

###先使用makeblastdb建库，如果是老版本的blast的话，则使用formatdb
makeblastdb -in 016_new.pep -dbtype prot
###使用blast进行alignment，得到一个表格，输出格式为m6
blastp -query 017_new.pep -db 016_new.pep -evalue 1e-5 -max_target_seqs 1 -num_threads 2 -out 1_2_blastp_out.m6 -outfmt 6
###而后使用shell指令进行简单的排序跟取值，如果存在多个物种比对，则需要写更为复杂的脚本去取双向最优序列
cut -f1-2 1_2_blastp_out.m6 | sort | uniq > 1_2.homolog

现在我们就得到了双向最优配对序列的序列名，之后需要使用ParaAT进行后续操作，ParaAT的下载链接为https://ngdc.cncb.ac.cn/tools/paraat

首先将各个样品的cds文件和pep文件分别cat起来和成一个总文件。
同时注意，直接用transdecoder输出的cds跟pep因为序列名后半部分会不相同，所以需要先更改成统一的ID格式。
ParaAT需要全局调用从能正常使用，同时要确保clustalw已经安装，不然会报错。
-p是指定线程数的文件，直接生成一个文件里面填个数字即可
之后运行程序

ParaAT.pl -h 1_2.homolog -n all.cds -a all.pep -m clustalw2 -p proc -f paml -o yeast_paml

运行发现报错，需要把安装路径添加到环境变量中

image.png

运行成功会生成一个文件夹yeast_paml，里面就是paml的输入格式，文件内容如下，数字2表示这里面有两个比对序列，1008是序列长度

image.png

如果引入外类群之类的或者多个物种相互比对，则会出现下面的情况

image.png

然后我们需要把yeast_paml里面的文件总数统计一下

ls -l yeast_paml/|wc -l
30830
###即我们里面有30830个配对文件，在后面的计算中需要将这个数字填进paml的配置文件种

###之后将所有的文件cat到一起
cat yeast_paml/*.paml > test.codon

paml的输入文件有三个，seq文件，即我们上面生成得到的test.codon，tree文件，即各物种间的关系，另外一个是配置文件，里面就是调整我们输入数据，输出数据，使用模型，进化速率等等的调整。

第一个，seq文件，同源直系基因整合文件即上面paml文件整合的test.codon

第二个，tree文件，即各物种间的相互关系

第三个，codeml.ctl，配置文件

可以直接从安装目录拷贝一份并更参数，文件内容如下

     seqfile = stewart.aa * sequence data filename
     treefile = stewart.trees      * tree structure file name
      outfile = mlc           * main result file name

        noisy = 9  * 0,1,2,3,9: how much rubbish on the screen
      verbose = 1  * 0: concise; 1: detailed, 2: too much
      runmode = 0  * 0: user tree;  1: semi-automatic;  2: automatic
                   * 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise

      seqtype = 2  * 1:codons; 2:AAs; 3:codons-->AAs
    CodonFreq = 2  * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table

*        ndata = 10
        clock = 0  * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
       aaDist = 0  * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
   aaRatefile = dat/jones.dat  * only used for aa seqs with model=empirical(_F)
                   * dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own

        model = 2
                   * models for codons:
                       * 0:one, 1:b, 2:2 or more dN/dS ratios for branches
                   * models for AAs or codon-translated AAs:
                       * 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
                       * 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)

      NSsites = 0  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
                   * 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
                   * 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                   * 13:3normal>0

        icode = 0  * 0:universal code; 1:mammalian mt; 2-10:see below
        Mgene = 0
                   * codon: 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff
                   * AA: 0:rates, 1:separate

    fix_kappa = 0  * 1: kappa fixed, 0: kappa to be estimated
        kappa = 2  * initial or fixed kappa
    fix_omega = 0  * 1: omega or omega_1 fixed, 0: estimate 
        omega = .4 * initial or fixed omega, for codons or codon-based AAs

    fix_alpha = 1  * 0: estimate gamma shape parameter; 1: fix it at alpha
        alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
       Malpha = 0  * different alphas for genes
        ncatG = 8  * # of categories in dG of NSsites models

        getSE = 0  * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 1  * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)

   Small_Diff = .5e-6
    cleandata = 1  * remove sites with ambiguity data (1:yes, 0:no)?
*  fix_blength = 1  * 0: ignore, -1: random, 1: initial, 2: fixed, 3: proportional
       method = 0  * Optimization method 0: simultaneous; 1: one branch a time

* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt., 
* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt., 
* 10: blepharisma nu.
* These codes correspond to transl_table 1 to 11 of GENEBANK.

将* ndata = 10 这一行中的 星号去掉 后面的数字改为刚刚得到的基因对数
之后就可以运行paml了

codeml codeml.ctl

之后得到文件

image.png

其中dN,dS,rst,t这些为结果文件，根据需要去画图。