Codeml（PAML）を用いた dN dS の算出

インストール（Linux）

参考：​http://abacus.gene.ucl.ac.uk/software/#paml-for-unixlinux

# ダウンロード,インストール
git clone https://github.com/abacus-gene/paml.git
tar -xvf paml4.9j.tar
cd paml4.9j
rm bin/*.exe
cd src
make -f Makefile
ls -lF
rm *.o
mkdir ../bin
mv baseml basemlg chi2 codeml evolver infinitesites mcmctree pamp yn00 ../bin

# テストラン
cd ../
bin/codeml
Bash

※ なぜか pamlのフォルダの中からしか動かせなかった（おそらくPATHの問題だが解決できなかった）

ファイルの準備｜アライメント済みのcds配列データ（phylip形式）

Fastaファイルは使用できない。​改行とスペースに関してやや不思議？なルールがある（下例）。

> in.phy
2 12
#ここに改行必須
Sp1  ATGAGGCCTTTC
Sp2  ATGGGAT-TTCT
#配列名と配列の間に２つ以上のスペース必須

cds配列だけでアライメントを取ると読み枠が考慮されず、dNdSの計算が正しくされない。
そのためペプチド配列を参考にしてアライメントを取る必要がある。
以下、その方法の参考。​ついでにfastaからphylipにする（力技）。
依存：tranalign, mafft

# 環境作る
conda create -n alignforcodeml
conda install -c bioconda seqkit emboss mafft -y

# ペプチド配列参考にアライメント
## NやXがあると動かない
mafft --thread 10 pep.fa > pep.a.fa
tranalign cds.fa pep.a.fa cds_basedpep_a.fa

# fastaからphylipにする（cds_basedpep_a.fa > in.phy）
## 改行なしfasta
awk '{if ($1 ~ /^>/) print "\n"$1; else printf $1}' cds_basedpep_a.fa | \
    sed -e '1d' | grep -v ">" > nofold.tmp
## 配列数・配列長のカウント
seqkit stat cds_basedpep_a.fa  > stat.tmp
t=$(awk '{print $4}' stat.tmp | tail -n 1)
l=$(awk '{print $7}' stat.tmp | tail -n 1 | sed 's/,//g')
echo -e "$t $l\n" > header
## 整形
seqkit fx2tab -n cds_basedpep_a.fa > Name
paste -d " " Name nofold.tmp | sed -e 's/ /    /g' > phylip.tmp
cat header phylip.tmp > in.phy
rm *tmp
Bash

基本の実行

bin/codeml codeml.ctl
Bash

※ codeml. ctlの中身を解析に合わせて書き換えて実行
※ 参考：​https://academic.oup.com/mbe/article/40/4/msad041/7140562?login=true
※ 参考：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf​

> codeml.ctl
seqfile = in.phy                   * sequence data file name
outfile = mlc                      * main result file name
treefile = stewart.trees           * tree structure file name

noisy = 9                          * 0,1,2,3,9: how much rubbish on the screen
verbose = 0                        * 1: detailed output, 0: concise output
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
aaDist = 0                         * 0: equal, +: geometric; -: linear, 1-6: G1974, Miyata, c, p, v, a
                                   * 7: AAClasses
aaRatefile = wag.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
                                   * 5: FromCodon0, 6: FromCodon1, 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-11: see below

Mgene = 0                          * 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 = 0.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 = 3                          * # of categories in dG of NSsites models
fix_rho = 1                        * 0: estimate rho; 1: fix it at rho
rho = 0.                           * initial or fixed rho, 0: no correlation

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

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

ペアワイズでdNdSの算出

> codeml.ctl
seqfile = in.phy                 * 配列データ
seqtype = 1                      * cdsだから１
outfile = results.txt            * 出力ファイル名
noisy = 9                        * 詳細な出力を表示
verbose = 1                      * 詳細レベル
runmode = -2                     * ペアワイズの比較モード
model = 2                        * ペアワイズの時ここ意味あるのかわからん
Bash

系統樹を用いたdNdSの算出

> codeml.ctl
seqfile = in.phy                 * 配列データ
seqtype = 1                      * cdsだから１
treefile = in.nwk                * 系統樹データ
outfile = results.txt            * 出力ファイル名
noisy = 9                        * 詳細な出力を表示
verbose = 1                      * 詳細レベル
runmode = 0                      * 系統樹を使用するモード
model = 1                        * すべての枝で独立なomega
NSsites = 0                      * サイトごとに同一なomega
Bash

コメント