paml计算dn／ds的小实例
更新：20190319
使用EasyCodeML完成dnds的计算计算，自己尝试了分支模型和位点模型的小实计算，简直不要太容易！计算
更新：20190314 教程 http://wiki.cathdb.info/wiki-beta/doku.php?小实id=tutorials:eccb_t2_codeml
- The selevtive presure in protein coding genes can be detected within the framework of comparative genomics.
- The selective presure is assumed to be defined by the ratio (w) dN/dS.
- dS represents the synonymous rate (changing the amino acid)
- dN represents the non-synonymous rate (keeping the amino acid)
- Under purifying selection, natural selection prevents the replacement of amino acids, so the dN will be lower than the dS.
- Values of w<1, =1, and > 1 means negative purifying selection, neutral evolution, and positive selection.
（这个视频教程不是我的，我也只是计算跟着这个教程来学习的！）
以下内容全部在window系统操作
使用到的小实软件
- MEGA 序列比对；进化树构建
下载地址 https://megasoftware.net/ - EasyCodeML 比对格式转化为PAML要求的格式
下载地址 https://github.com/BioEasy/EasyCodeML - PAML 运用分支模型估计dn/ds
下载地址 http://abacus.gene.ucl.ac.uk/software/paml.html (win系统不用编译，下载解压出来直接可以使用)
第一步：下载序列
根据教程中提到的计算物种拉丁名和myoglboin为关键词在NCBI中检索获得本次分析中使用到的8条肌红蛋白序列，放到同一个文件中，小实手动去掉终止密码子（有几条序列中不包括终止密码子，计算不知道为什么）（需要原始序列的小实同学可以给我留言）
第二步：codeML输入文件的准备
包括：序列比对；进化树构建，控制文件编辑等
- 使用MEGA进行多序列比对
读入数据，计算选择Alignment——Align by Muscle (Codons)
将比对结果导出为fasta格式 - 构建NJ树
读入上一步比对好的小实数据，构建NJ树：bootstrap1000,计算kumra两参数模型image.png
生成的结果树tree.PNG
将生成的树文件导出为newick格式 - 使用EasyCodeml将比对结果转化为PAML要求的格式
Tools——Seqformat convertor
将上一步比对好的fasta文件拖入，第二个框输出格式选择PAML即可获得结果
参考EasyCodeML帮助文档对进化树分支进行标记 - 按照视频中的小实内容编辑控制文件
image.png
tree.PNG第三步：计算dnds
win+R快捷键输入cmd调出命令行窗口；cd命令切换到保存输出结果的文件夹下，分布输入codeml和控制文件所在的计算路径，回车（win下也可以拖拽）
- 零假设 所有支系共有一个omage值
控制文件内容
这里的参数很多都不明白，前几天大体看了一下paml的文档，里面有很详细的介绍，需要花时间多看
1.PNG关键输出结果
这里的输出结果和视频教程中的略有差别，视频中的ntime和np分别是14和16，暂时还没有搞懂问题出在哪里
- 备择假设：特定支系拥有不同的omega值
只需要将树文件替换为标记分支的树文件，model值改为2即可
lnl_null - lnl_alternative 符合卡方分布
使用R语言的phisq函数计算显著性
自己得到的结果是1，没搞清楚问题出在哪里
步骤应该就是这么个步骤，关于结果的解读和其他模型的设置多看paml的帮助手册
第一种方法
下载地址 http://abacus.gene.ucl.ac.uk/software/paml4.9h.tgz
安装教程 http://abacus.gene.ucl.ac.uk/software/paml.html
第二种方法
直接用conda，conda的安装方法可以参考公众号 生信技能树 的文章价值999的全外显子教学视频--免费送
conda安装好以后直接使用命令
将CDS序列进行基于密码子的比对，这一步可借助mega完成，然后导出比对结果为fasta格式，然后按照视频教程去掉终止密码子并整理成所需要的格式（所需要的格式第一行是tab+数字+tab+数字，第一个数字是几条序列，第二个数字是每条序列有几个位点，自己获得以上信息的做法是将fasta格式的比对文件转换成phylip格式,可以实现比对文件格式转换的小工具http://sing.ei.uvigo.es/ALTER/）；树文件可以用IQ-tree获得，IQ-tree下载及安装http://www.iqtree.org/#download
A、比对：打开mega, File——Open a File 选择文件，提示
1.PNG选择Align,弹出
2.PNG然后选择Alignment——Alignment by ClustalW(codons),然后点OK在点OK,之后Data——Export Alignment——Fasta format,比对完成，之后将比对结果转换成要求的格式
3.PNG这里需要注意的地方是序列名后有两个空格
B、建树：IQ-tree解压出来后bin目录下有可执行的程序iqtree,运行
获得帮助信息
4.PNG建树用到的命令是
C、按照视频教程的内容修改控制文件的内容，然后执行程序即可。获得的结果
5.PNG以上分析用到的CDS序列文件https://pan.baidu.com/s/1xyNj_cvR68An-Lkf2HGBgQ
发现一个小问题：如果使用conda install paml 安装paml 好像找不到需要的控制文件在哪
以上内容基本在linux系统操作完成，如果没有linux系统，目前想到三种解决办法：
1、如果自己电脑是win10系统可以直接安装linux系统，百度搜索相关教程即可；
2、购买腾讯云，学生套餐10块每月还挺划算的，1核2G运行内存，50G存储空间；
3、安装虚拟机，自己之前也尝试过，个人感觉稍微有些麻烦。
关于linux系统的入门学习个人强推 网易云课堂 细说linux 视频教程
配置文件的参数修改
前三个参数用来指定文件位置
runmode = 0 user tree
seqtype = 1 codons
ndata = 自己暂时还不知道这个参数是用来做什么的，视频教程里的配置文件没有这个参数
model = 0
https://www.plob.org/article/6943.html codeml的配置文件
https://www.bilibili.com/video/av10469605?from=search&seid=8859495264060497812 b站教学视频
微信文章： Ka/Ks比值：核酸分子进化的指标
微信文章： 从人见人爱的向日葵说起——Ks与全基因组多倍化事件
微信文章： 价值999的全外显子教学视频--免费送
微信文章： 基于ML法重建系统发育树-IQ-TREE篇
http://www.biotrainee.com/thread-1685-1-1.html 生信技能树
https://wenku.baidu.com/view/39c6d6f30029bd64783e2c60.html 百度文库