使用PAML进行正选择基因分析

PAML (Phylogenetic Analysis by Maximum Likelihood),是杨子恒开发的一款利用DNA或Protein数据使用最大似然法进行系统发育分析的软件。

该软件不擅长构树,但能用于评估系统进化过程中的参数和假设检验,有利于PAUP、PHYLIP、MOLPHY、PHyML和RAxML等其它软件的构树。 PAML软件的baseml和codeml命令中用于搜索系统发育树的算法非常原始且效果较差,当数据量较大,例如物种数多于10个时,作者推荐使用其他更好的软件来构建系统发育树。 软件的详细用法请参考PDF说明文档常见问题

下载并安装PAML软件

$ wget http://abacus.gene.ucl.ac.uk/software/paml4.9i.tgz
$ tar zxf paml4.9i.tgz -C /opt/biosoft/
$ cd /opt/biosoft/paml4.9i/
$ rm bin/*
$ cd src
$ make -f Makefile
$ cp baseml basemlg chi2 codeml evolver infinitesites mcmctree pamp yn00 ../bin
$ echo 'PATH=$PATH:/opt/biosoft/paml4.9i/bin' >> ~/.bashrc
$ source ~/.bashrc

PAML软件包含多支程序,来实现不同的系统发育分析功能:

baseml 用于对核酸序列进行最大似然法分析。
codeml 用于对密码子和蛋白序列进行最大似然法分析。该程序整合了以前的两支程序codonml和aaml,前者用于对密码子序列进行分析,后者对蛋白序列进行分析。在新的codeml程序的配置文件中设置seqtype=1,表示使用codonml命令,设置seqtype=2,表示使用aaml命令。
evolver 用于模拟核酸、密码子或蛋白序列。
basemlg 相比于baseml,该程序应用了连续的gamma模型。当数据中的物种数多于6到7个时,该程序运行非常慢以至得不到结果。这时,需要使用设置离散gamma模型的baseml命令。
mcmctree 该程序应用Bayesian MCMC算法计算物种分歧时间。
pamp 用于简约法分析。
yn00 用于对两条编码蛋白的DNA序列进行比较并计算dnds。
chi2 用于在likelihood ratio test中计算卡方临界值和p值。

PAML使用前的注意事项:

1. PAML要求的输入文件是多序列比对结果。PAML软件不能进行多序列比对,需要使用其他软件进行分析。
2. PAML软件进行基于密码子的分析要求输入的序列必须是去除了Introns等非编码区的DNA序列,且其长度是3的整数倍,且第一个核酸位点是密码子的第一位。PAML软件不能进行基因预测,密码子的比对结果需要自行编写程序分析得到。
3. 对于较大数据(物种数>=10)量的分析,输入的树文件最好使用其他软件计算。

1. 使用codeml和进行正选择分析

CodonFreq

程序根据输入的序列计算出密码子3个位点上各4种碱基的频率,例如:

Codon position x base (3x4) table, overall
position 1: T:0.20110 C:0.13736 A:0.32308 G:0.33846
position 2: T:0.20000 C:0.17143 A:0.33077 G:0.29780
position 3: T:0.32747 C:0.22857 A:0.25165 G:0.19231
Average T:0.24286 C:0.17912 A:0.30183 G:0.27619

然后再计算64种密码子的频率。计算方法有四种:
(1)设置CodonFreq参数值为0,表示除三种终止密码子频率为0,其余密码值频率全部为1/61;
(2)设置CodonFreq参数值为1(F1X4),表示除三种终止密码子频率为0,计算61种密码子时,使用上面Average的值进行计算,例如TTC密码子的频率为0.24286*0.24286*0.17912;这61种密码子频率之和不等于1,然后对其数据标准化,使其和为1,得到所有密码子的频率。
(3)设置CodonFreq参数值为2(F3X4),表示除三种终止密码子频率为0,计算61种密码子时,使用上面三种position的的值进行计算,例如TTC密码子的频率为 0.20110*0.20110*0.22857;这61种密码子频率之和不等于1,然后对其数据标准化,使其和为1,得到所有密码子的频率。
(4)设置CodonFreq参数值为3(codon table),表示直接使用观测到的各密码子的总的频数/所有密码值得总数,得到所有密码子的频率。

一般选择第三种方法,设置CodonFre的值为2。第一种方法让所有密码子频率相等,不符合实际轻卡;第二种方法没有考虑到三种位点各碱基频率的差异,得到的密码子频率不准确;第三种方法能比较准确计算出各密码子的频率;第四种方法由输入数据得到真实的各密码子频率,但有时候有些非终止密码子的频率为0,可能不利于后续的计算。

Linux Top命令常用技巧

Linux系统中top命令类似Windows系统的资源管理器,能查看当前系统运行状况。本文介绍top命令使用的一些小技巧。最全面的top使用方法参考其说明文档:

$ man top

1. top页面信息

第1行:显示当前时间、开机时间、用户登录数量、系统负载(1分钟、5分钟、15分钟)。

第2行:显示任务运行数量。

第3行:显示CPU消耗状态。us和sy之和等于总的CPU消耗百分比。

第4行:物理内存消耗情况,最后一个数值表示剩余的物理内存数量。

第5行:虚拟内存消耗情况。

第6行:用来输入命令对top显示方式进行设置。

第7行:top主体部分的表头。

>=8行:top主体部分显示的用户进程运行监视状态。

2. top命令小技巧

按小写字母t,能切换CPU消耗的显示方式,将数字显示换成条形显示。

按小写字母m,能切换内存消耗的显示方式,将数字显示换成条形显示。

按小写字母c,则将最后一列COMMAND信息显示出完整的命令行。再按一次字母c,则切换回仅显示进程的程序名。

按小写字母d,再输入一个浮点数,例如1.0,则能修改刷新时间,表示每隔1.0秒刷新一次进程运行信息。默认每3秒刷新一次。

按小写字母u,再输入某用户名,则仅显示该用户名下的进程信息。

按大写字母P,则根据CPU消耗进行排序; 按大写字母M,则根据内存消耗进行排序;按大写字母N,则根据PID编号进行排序;按大写字母T,则按照消耗的计算时间进行排序。

按小写字母z,能让显示带颜色。

按小写字母f,进入top命令的列管理页面,能实现较丰富的功能:(1)指定用于排序的列。使用方向键的上和下移动标注栏,比如移动到PID,然后按小写字母s,表示按PID进行排序;再按q键推出列管理页面,返回到top的进程监视页面,这是发现是按照PID进行了排序;同样道理,可以按照任何一列进行排序。(2)更改各列的顺序。在列管理界面,按方向键右,则可以将标注栏扩大,这时再按方向键上和下,可以移动标注对象,用于改变进程监视中各列的位置。(3)选择显示更多列或减少不需要的列。按方向键上或下移动标注栏,按空格键选中或反选指定列,从而增加或减少列。

按大写字母X,再输入数字3,则能增大每列的宽度,可以让USER栏显示全用户名。

按数字键2,可以显示每颗物理CPU的使用率;按数字键1,则显示每个CPU线程的使用率;按数字键3,再输入数字选择物理CPU,以显示对应CPU的各线程使用率。

解决Virtualbox安装CentOS7鼠标左键失效情况

使用Virtualbox安装CentOS7系统,在安装增强包工具后,鼠标左键在桌面窗口中失灵。此时,需要在root权限下安装kernel软件,重启后恢复正常。

# yum install -y kernel*
# shutdown -r now

也有人说CentOS7安装完毕默认安装了open-vm-tools软件。将其默认安装的该软件删除后,再安装增强包工具后重启,能解决该问题。

R软件常规操作

在Linux系统上安装R软件

wget https://mirrors.tuna.tsinghua.edu.cn/CRAN/src/base/R-3/R-3.5.3.tar.gz -P ~/software/
tar zxf ~/software/R-3.5.3.tar.gz
cd R-3.5.3/
./configure --prefix=/opt/sysoft/R-3.5.3
make -j 4
make install
cd .. && rm R-3.5.3/ -rf
echo 'PATH=/opt/sysoft/R-3.5.3/bin/:$PATH' >> ~/.bashrc
source ~/.bashrc

自己设置R包的安装源并安装R包

设置CRAN源,用于安装R包。

R
> options(repos="http://mirrors.tuna.tsinghua.edu.cn/CRAN/")
> install.packages(c('ggplot2', 'gplots', 'ape'))
> install.packages("devtools", repos="http://mirrors.tuna.tsinghua.edu.cn/CRAN")

设置 Bioconductor源,用于安装生信R包。

> source("http://bioconductor.org/biocLite.R")
> options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
> biocLite(c('edgeR', 'fastcluster', 'limma', 'DESeq2', 'Biobase', 'qvalue', 'seqLogo', 'pathview'))

下载并安装软件包时提示错误status was ‘SSL connect error’

原因是nss版本较低,更新nss即可。

# yum -y update nss