选择压力分析，从序列下载到结果绘图
选择压力及基因在进化过程中所受到的压力，也称自然选择。
- 在遗传学中， Ka/Ks或者dN/dS表示的是非同义突变频率(Ka)和同义突变频率（Ks）之间的比例。
- 一般认为，同义突变不受自然选择，而非同义突变则受到自然选择作用。由于非同义替换往往对于生命体有害，所以在负选择的作用下，非同义替换常常会在群体中被逐渐消灭。
- 当Ka/Ks > 1时，基因受到正选择（positive selection）
- 受到正选择的基因可能是物种受到自然选择/人工选择的结果。对研究物种环境适应性、人工选育/驯化的影响具有重要意义。
具体知识我不作赘述，感兴趣的可以去看其他文章，本文主要介绍选择压力分析具体的流程及方法。这方面网上的教程很少，我也只是在各种网上的教程里自己摸索整理出的方法，方法不一定对，仅供参考。
1. 物种选择及序列下载
- 物种的选择
在进行选择压力之前，我们首先要确定分析的物种，物种一般选择近缘物种，数量不宜太多，一般5-6个。不然会导致单拷贝直系同源基因太少。
- 序列下载
我们需要下载物种的CDS序列或自己测序获得的CDS序列进行分析，下载的数据库推荐两个：
-
NCBI Genome Database,这个是大家平时可能最常用的基因数据库。
image.png
-
中国国家基因库，这个数据库回定期同步国际数据库的资源，网络加载速度比 NCBI 快很多且中文界面比较友好。
image.png
2. 序列去冗余及格式调整
- 这里我以Ameiurus melas基因组为例，首先下载CDS序列，如下图所示，我们首先需要修改其序列名称，只保留红框中我们需要的部分，这步可以自己编写脚本，或者使用我写的Python脚本。（使用方法： python 1_modify_NCBI_ID.py YourFastaFile.fa）
image.png
然后对序列最长CDS序列进行提取，并将物种名称加到ID的前面，得到如下格式：
image.png
- 另外加上自己组装的Trinity序列，首先从中提取最长转录本，方法见我另一篇文章用Python提取FASTA中最长转录本,或者直接使用我写好的脚本。再使用TransDecoder预测CDS序列。最后提取最长CDS序列，原理和上面类似。
3. 直系单拷贝同源基因筛选及系统发育树构建
这里使用OrthoFinder寻找同源基因并建树,得到直系单拷贝同源基因和物种系统发育树结果。
4. 多序列比对及格式转换
-
首先将所有的直系单拷贝同源基因进行多序列比对，这里我这了很多软件，发现常规的mafft, muscle 等软件比对后会改变密码子序列。于是使用了PhyloSuite中的mafft多序列比对，选择密码子（codon）模式。
image.png
-
最后使用PhyloSuite的格式转换工具，转换比对完成后的序列为Paml所需要的格式。到这里，所有的数据准备就完成了，就可以开始选择压力分析了。
image.png
5. 选择压力分析
目前我们有了选择压力分析所必须的文件：
- 比对好的Paml格式密码子序列。
- 物种的系统发育树。（计算时要在选择的物种名后加上 #1）
选择压力分析使用的是Paml中的codeml软件包。EasyCodeML做了很好的可视化界面，但是不能批量运算，只能一次分析一个，并且容易卡死（Mac 和 windows上都遇到过）。
所以我自己写了一个脚本callCodeml（使用方法：python3 callCodeml.py 序列文件夹 物种发育树文件）来计算前景枝受到的选择压力，调用codeml批量运算枝位点模型（Branch site model），并使用chi2检验，只输出 P < 0.05 的结果，最终将结果自动整理成表格。
image.png
至此，选择压力分析完成，我们只需要根据OrthoFinder的文件即可找到OG00XXXX同源组所对应的基因名称。
6. 绘图
困了，睡觉，改天有缘再写。
