codml 位點模型 (Site model) 假設蛋白質氨基酸序列或密碼子序列中不同位點受到的選擇壓力不同，可以通過 model = 0和多個Nssites 參數(shù)來實現(xiàn)這種情況下dn ds 的計算，進行正選擇分析。

相應家族CDS序列提取和對齊

OrthoFinder 已經對不同物種的氨基酸序列進行Orthologs 和Paralogs 的分類，利用Orthogroups/Orthogroups.txt對目標家族和相應基因的 cds 序列文件進行索引并抽取構成ogs.fa 核苷酸序列文件，每個ogs.fa文件只包含該家族的cds序列。

通過prank 進行比對，必要時可以用gblocks 進行過濾，PAL2NAL.pl將蛋白質的多序列聯(lián)配結果轉換為編碼序列的多序列聯(lián)配，并去除多序列聯(lián)配中產生的空位(gap)，將fasta 格式轉換為phily 格式，同時用raxml 進行進化樹重構，至此codml 所需文件就構建完成。

#  extracting cds fasta of each OG id
multfa2OG.py --OGs OGS --multfa cds.fa --path results

# alignmented nucleotide sequence and makeing tree
mkdir -p PAML/Prank
mkdir -p PAML/Pal2Nal
mkdir -p PAML/NucTree

prank -d="${ogs}.fa" -translate -F -o=PAML/Prank/${ogs} 
pal2nal.pl PAML/Prank/${ogs}.best.pep.fas PAML/Prank/${ogs}.best.nuc.fas -output fasta -nogap -nomismatch > PAML/Pal2Nal/${ogs}.nuc.fa

cd PAML/NucTree 
raxmlHPC-PTHREADS -f a -m GTRGAMMA -T 2 -x $RANDOM -N 100 -n ${ogs}.tree -s ../../raxmlPAML/Pal2Nal/${ogs}.nuc.fa -p $RANDOM

# makeing Phylip fasta
mkdir -p PAML/Phylip
python3 fasta2phylip.py --input Pal2Nal/${ogs}.nuc.fa --output Phylip/${ogs}.nuc.fa.paml

# codeml.ctl
noisy = 1
verbose = 0
runmode = 0
seqtype = 1
CodonFreq = 2
model = 0
NSsites = 0 1 2 7 8
fix_omega = 0
omega = 1

結果處理

正選擇的相關信息：

The first compares M1a (NearlyNeutral) and M2a (PositiveSelection), while the second compares M7 (beta) and M8 (beta&ω). M1a (NearlyNeutral) and M2a (PositiveSelection) are slight modifications of models M1 (neutral) and M2 (selection) in (Nielsen and Yang 1998).

一半來說，提取結果文件m1, m2, m7, m8 模型的lnL數(shù)據(jù)，通過自由度為2的卡方檢驗是更為常見的做法。

names(codeml_out) <- c("Ogs", "m1_ln1", "m2_lnl", "m7_lnl", "m8_lnl")

compare_models <- function(lnl1, lnl2){
  likelihood = 2*(abs(lnl2 - lnl1))
  p = 1 - pchisq(likelihood, df = 2)
  return(p)
}

# M1 (neutral) vs M2 (selection)
## alternative hypothesis: true selection happened.
# M7 (beta) vs M8 (beta& w)
# The M1a-M2a comparison appears to be more stringent than the M7-M8 comparison.
m2_p_pos = compare_models(codeml_out$m1_ln1, codeml_out$m2_lnl)
m8_p_pos = compare_models(codeml_out$m7_lnl, codeml_out$m8_lnl)

padj_M2vsM1 = p.adjust(m2_p_pos, "BH", length(m2_p_pos))
padj_M8vsM7 = p.adjust(m8_p_pos, "BH", length(m8_p_pos))