封禁一些IPs的SSH登录，保护服务器安全

发表于由
回复

登录服务器的root用户时，可能会收到如下提示：

There were 3225 failed login attempts since the last successful login.

这表示有很黑客在尝试用一些用户和密码登录服务器。为了保障服务器安全，则可以在日志文件/var/log/secure中查看攻击的IPs来源和其失败登录的次数。例如：

Invalid user neisius from 211.144.12.75 port 34041
Failed password for invalid user server from 219.94.99.133 port 15139 ssh2
pam_unix(sshd:auth): authentication failure; logname= uid=0 euid=0 tty=ssh ruser= rhost=106.12.217.39

若某个IP次数很高，则将其IP信息添加到/etc/hosts.deny文件中进行封禁。例如，修改/etc/hosts.deny文件内容：

sshd:101.109.83.140
sshd:103.27.238.202
sshd:104.236.94.202
sshd:106.12.200.13
sshd:106.12.217.39
sshd:106.54.19.67
sshd:111.230.73.133
sshd:117.50.45.254
sshd:118.24.38.53
sshd:122.155.223.119

为了能自动地封禁IPs，编写名为IP_denied_for_ssh_security.pl的perl程序：

#! /usr/bin/perl
use strict;

my $usage = <<USAGE;
Usage:
    $0 <num>

    运行程序后，程序读取/var/log/secure文件中信息。若存在某个IP登录失败(含有Invalid user|Failed password|authentication failure等信息)次数超过指定的次数，则是需要封禁的IP。此外，程序再读取/etc/hosts.deny文件中的IP信息，得到所有需要封禁的IP。然后程序将这些IP信息全部写入到/etc/hosts.deny文件中进行封禁。
    程序需要以root权限运行，因为只有root权限才能读取/var/log/secure文件，对/etc/hosts.deny进行写入。
    程序运行完毕后，会清空/var/log/secure文件内容。

USAGE
if(@ARGV==0){die $usage}

open IN, "/var/log/secure" or die "Can not open file /var/log/secure, $!";
my %num;
while (<IN>) {
    if (m/Invalid user/ or m/Failed password/ or m/authentication failure/) {
        $num{$1} ++ if m/(\d+\.\d+\.\d+\.\d+)/;
    }
}
close IN;

my %deny;
my $deny_num = 0;
foreach (keys %num) {
    if ($num{$_} >= $ARGV[0]) {
        $deny{$_} = 1;
        $deny_num ++;
    }
}
print STDERR "$deny_num IPs were detected for blocking from file /var/log/secure.\n";

open IN, "/etc/hosts.deny" or die "Can not open file /etc/hosts.deny, $!";
while (<IN>) {
    $deny{$1} = 1 if m/(\d+\.\d+\.\d+\.\d+)/;
};
close IN;

open OUT, ">", "/etc/hosts.deny" or die "Can not create file /etc/hosts.deny, $!";
my $total_num = 0;
foreach (sort keys %deny) {
    print OUT "sshd:$_\n";
    $total_num ++;
}
close OUT;
print STDERR "Total $total_num IPs were banned by file /etc/hosts.deny for SSH security.\n";

open OUT, ">", "/var/log/secure" or die "Can not  create file /var/log/secure, $!";
close OUT;

使用root用户执行该perl程序，若登录失败次数>=3次，则封禁对应的IP。

IP_denied_for_ssh_security.pl 3

基因正选择分析原理

发表于由
2

一、 正选择分析的目的：

  1. 两两基因的密码子序列进行比较，从而计算dN/dS，即omega（ω）值。若该值<1，则表示纯化选择；omega = 1，则中性进化；omega > 1，则正选择。若分析基因在两个物种中的序列，可以计算dN/dS的值，若omega > 1，即表明该基因在物种进化过程中，即由其祖先物种分化成这两个物种时，基因受到了正选择。对于两个物种/序列的正选择分析，比较简单。而实际情况中，要分析的物种数量很多，包含多个类群。这个时候的正选择分析相对复杂些。
  2. 对多个物种的基因序列进行正选择分析，若仍然按照两个物种时的要求，即分析该基因在物种进化中是否受到了正选择？这种结果可能不好说清楚。因为该基因可能在某一类群中序列很相似，其两两比较时，omega <= 1；而在另外一类群中两两比较时，很多时候omega > 1。最后软件可以从总体上给一个omega值，该值不可以拿来简单地评价该基因是否受到了正选择。所以，对多个物种进行正选择分析时，没法直接评价该基因是否受到了正选择。正选择只有在进行两两序列比较的时候，才能计算omega值，从而得到结果。
  3. 对基因在多个物种上的正选择分析，分析的目的则是：比较某个分枝上祖先节点和后裔节点（可以理解成，对无根树上某分枝两侧的两组物种进行比较，依然属于两两比较），从而计算该分枝的omega值。而在实际数据中，基因在不同的进化分枝上具有不同的omega值，同时在序列不同的位点也具有不同的omega值。目标分枝两侧的物种数量较多时，可以对序列上的每个位点进行omega值分析，从而鉴定正选择位点。所以，对基因在多个物种上的正选择分析，需要同时分析分析目标分枝的omega值和序列位点的omega值，从而判断基因是否受到正选择压。

二、使用PAM对基因进行正选择分析，有三种方法：

  1. PAML site model: 主要用于检测基因中的正选择位点。该方法分析时，认为进化树中各分枝的omega值是一致的，并比较两种模型：(1)模型m1是null model，认为所有位点的omega值<1或=1; (2)模型m2是正选择模型，存在omega <1、=1或> 1的位点。比较两个模型的似然值（lnL）差异，利用卡方检验（自由度为2）算出p值。若p值 < 0.05，则否定null model，认为存在正选择位点。此外，推荐采用比较模型m7和m8，它们将omega值分成了10类，其p值结果比上一种比较方法更宽松，能检测到更多的正选择基因。使用PAML site model方法能在整体水平上检测基因的正选择位点，而不能表明基因在某个进化分枝上是否受到正选择压。
  2. PAML branch-site model: 主要用于检测基因在某个进化枝上是否存在的正选择位点。该分析方法认为目标分化枝具有一个omega值，其它所有分枝具有一个相同的omega值，然后再检测正选择位点。同样对两种模型进行比较：（1）第一种模型为模型2，将omega值分成<1、=1、>1的三类，这和site model中的一样；（2）第二种模型和前者一致，只是将omega固定成1，作为null model。比较两种模型的似然差异，利用卡方检验（自由度为2）算p值（chi2命令算出的值除以2）。若p值< 0.05，则能通过Bayes Empirical Bayes (BEB)方法计算正选择位点的后验概率，若存在概率值 > 0.95正选择位点，则表示基因在目标分枝上受到正选择压。PAML软件在branch-site模式下，并不给出分枝上的omega值。这表示branch-site模式虽然考虑了目标分枝上具有不同的omega值，但仍然以分析位点上的omega为主。值得注意的是，在branch-site模式下可能检测到正选择位点，但在目标分枝上的omega值仍然可能低于1。可能软件作者基于这点考虑，就没有给出目标分枝上的omega值，以免影响一些人对正选择结果的判断。
  3. PAML branch model: 主要用于检测在某个分枝上，其omega值是否显著高于背景分枝，即基因在目标分枝上进化速度加快。该方法认为基因序列上所有位点的omega值是一致的，对两种模型进行比较：（1）第一种模型为null model，所有分枝具有相同的omega值；（2）第二种模型认为目标分枝具有一个omega值，其它所有分枝具有一个相同的omega值。比较两种模型的似然差异，利用卡方检验（自由度为1）算p值。若p值 <= 0.05，且目标分枝上的omega值高于背景值，则认为该基因为快速进化基因。一般情况下，该方法计算得到的p值会低于第二种方法的结果。

三、其它注意事项

Branch-site model相比于site model的优点是考虑了不同的分枝具有不同的选择压，即具有不同的omega值。该方法让目标分枝具有一个不同的omega值，并没有让所有分枝的omega值独立进行计算（理论上这样是最好的）。这样算法很复杂，程序运行非常非常消耗时间。但其实也没必要这样做，因为正选择分析其实是两条序列比较后，分析dN/dS，再找正选择位点，其分析结果就应该是某个分枝上基因是否受到正选择，在序列那个位点上受到正选择。

若在目标分枝上，其omega值小于1，但是却能找到正选择位点。即该基因在该分枝上的dN/dS < 1，但是在某些位点上，dN/dS > 1。那么该基因是否属于正选择基因？我认为：属于。之所以为正选择基因，主要是因为基因的个别位点或多个位点存在正选择。当只有个别位点受到正选择压时，而其它多个位点存在纯化选择时，可能导致整体上的omega值小于1。此时，该基因也应该是属于正选择基因。

四、参考文献中的正选择分析方法描述

Science文章（https://science.sciencemag.org/content/364/6446/eaav6202）中的正选择基因分析方法：To estimate the lineage-specific evolutionary rate for each branch, the Codeml program in the PAML package (version 4.8) (134) with the free-ratio model (model = 1) was run for each ortholog. Positive selection signals on genes along specific lineages were detected using the optimized branch-site model following the author’s recommendation. A likelihood ratio test (LRT) was conducted to compare a model that allowed sites to be under positive selection on the foreground branch with the null model in which sites could evolve either neutrally and under purifying selection.[ The p values were computed based on Chi-square statistics, and genes with p value less than 0.05 were treated as candidates that underwent positive selection. We identified PSGs at the ancestral branch of Ruminantia (table S22), the ancestral branch of Pecora (table S23), each ancestral family branch of Ruminantia (tableS24), and each ancestral subfamily branch of Bovidae (table S24). We also compared the dN/dS values of Ruminantia families with outgroup mammals (fig. S52).

Science文章（https://science.sciencemag.org/content/364/6446/eaav6202）中快速进化基因分析方法：The branch model in PAML was used, with the null model (model=0) assuming that all branches have been evolving at the same rate and the alternative model (model=2), allowing the foreground branch to evolve under a different rate. An LRT with df =1 was used to discriminate between alternative models for each ortholog in the gene set. Genes with a p value less than 0.05 and a higher ω value for the foreground than the background branches were considered as evolving with a significantly faster rate in the foreground branch.

https://journals.plos.org/plosntds/article?id=10.1371%2Fjournal.pntd.0007463文献中对正选择基因的分析方法：

A calculation of mutational rate ratio ω between two gene sequences was the basis for the positive selection analysis. The ω was calculated as a ratio of nonsynonymous to synonymous mutational rates. The ratio indicates negative purifying selection (0 < ω < 1), neutral evolution (ω = 1), and positive selection (ω > 1) [54]. A set of selected genes from complete genomes was tested relative to positive selection using the maximum likelihood method using the CODEML of the PAML software package [55]. PAML version 4 [56] and its user interface PAMLX [57] were used in our study. For each analyzed gene, its maximum likelihood phylogenetic tree was used as an input tree. The CODEML offers several different codon evolutionary models, and the statistical likelihood ratio test (LRT) was used to compare the codon evolutionary model to the null model. The Bayes empirical Bayes method (BEB) [58] was then used to evaluate the posterior probability of sites considered to have been positively selected.

The CODEML models could produce different results (i.e., a list of sites under positive selection) since they calculate different parameter estimates. Site models allow ω to vary in each site (codon) within the gene. Statistical testing was required for sites with ω > 1. Two pairs of models were predominantly used since their LRTs have low false-positive rates. M1a (nearly neutral evolution) was compared to M2a (positive selection) [58,59] and M7 (beta) was compared to M8 (beta & ω) [60]. Our preliminary testing found that the two model pairs gave the same or very similar results. Therefore we chose to use the M7-M8 model pair. The M7 model is a null model that allows 10 classes of sites with a ω beta-distribution within the interval 0 ≤ ω ≤ 1. Sites with ω > 1 are not allowed. The alternative M8 model adds an eleventh class of sites with ω > 1. Each site was tested to determine the class to which it belongs. The LRT compares twice the log-likelihood difference 2Δl = 2(l1-l0) between the M7 model (log likelihood value l0) and the M8 model (log likelihood value l1) to the χ2 distribution [61]. If the twice log-likelihood difference is above a critical χ2 value, then the null model is rejected, and the positive selection is statistically significant.

A considerable disadvantage of the site models is that ω was calculated as an average over all codons of the site. Therefore, the site models are not suitable for the data where ω also varies between lineages. In contrast, the branch-site models search for positive selection in sites and pre-specifies lineages where different rates of ω may occur [62]. Sequences of lineages are a priori divided into a group of foreground lineages where positive selection may occur and group of background lineages where only purifying selection or neutral evolution occurs. We used branch-site model A, which allows four classes of sites and different setups of foreground lineages to be tested depending on the gene phylogeny. In branch-site model A, all lineages under purifying selection with a low value of ω0 belong to site class 0. Weak purifying selection and neutral evolution with ω1 near to value 1 are allowed in site class 1. In site class 2a, a proportion of class 0 sites in foreground lineages is under positive selection with ω2 > 1. Similarly, site class 2b is a proportion of class 1 sites under positive selection with ω2 > 1. The null model for LRT has ω2 = 1. Critical values of LRT (2Δl) are 2.71 at 5% and 5.41 at 1% [63]. The posterior probabilities of suggested sites under positive selection were calculated using the BEB method.

使用squid搭建代理服务器

发表于由
回复

使用squid可以搭建代理服务器，可以用于翻墙访问国外网站。这里讲解通过在校内服务器上搭建squid代理服务，用于校外通过校内的代理服务器下载文献。

1. 在校内服务器上安装squid软件，并启动squid服务。

使用root权限在学校内网中安装有CentOS7系统的服务器上安装squid软件：

yum install squid

再修改squid配置文件，以让所有其它IP地址的用户可以访问校内服务器：

perl -p -i -e 's/http_access deny all/http_access allow all/' /etc/squid/squid.conf

启动squid服务：

systemctl start squid.service

2. 将代理服务器的3128端口映射到外网Aliyun服务器上的31280端口

squid代理服务器使用的默认端口是3128。若此时将内网服务器的防火墙开放3128端口，则在个人电脑上代理设置上填写内网服务器的固定IP和3128端口（见本文第4部分），在个人电脑上的浏览器访问网络时，就是通过代理服务器来访问了。

若是国外服务器上做的代理，则防火墙开放3128端口，可以在国内访问国外代理服务器进行翻墙。

firewall-cmd --add-port=3128/tcp --permanent 
systemctl restart firewalld.service

由于学校内网一般对外封闭了所有的端口，在校外依然是无法使用代理服务器的。因此，需要通过一个具有外网固定IP的服务器来进行中转。我这儿就使用了Aliyun上购买的一台具有固定IP的云服务器进行中转。中转时不需要校内的服务器开放防火墙3128端口。

中转的方法是使用ssh进行反向隧道连接，在内网服务器上输入反向隧道命令，将内网服务器的3128端口映射到Aliyun服务器的31280端口。这样当访问Aliyun服务器的31280端口时，就等同于在内网服务器上直接访问3128端口了。

ssh -f -N -R 31280:localhost:3128 chenlianfu@115.29.105.xxx

3. 在校外有固定IP的Aliyun服务器防火墙上开放31280端口的访问

使用root用户在Aliyun服务器上操作，开放防火墙31280端口，以利于个人电脑能访问Aliyun服务器的31280端口。

firewall-cmd --add-port=31280/tcp --permanent 
systemctl restart firewalld.service

4. 个人Windows10电脑上设置代理访问

按windows+x键，再按字母n，打开系统设置界面；点击网络和internet，再点击代理；在下半个页面的手动设置代理部分，打开使用代理服务器开关，然后在地址栏填写Aliyun服务器的IP地址115.29.105.xxx，在端口栏填写31280，最后点击保存。

然后打开浏览器，输入网址，即可通过校内的代理服务器访问网络。这时等同于在学校连接校园网后的网络访问。当用IE浏览器访问网址，若代理服务器没生效，IE浏览器会报告正在使用某个IP某个端口进行代理访问，代理不成功等信息。

成功后，示例访问文献https://doi.org/10.1016/j.molp.2020.05.005，若能查看文献全文和下载PDF全文，则表示成功。

使用PAML软件利用密码子序列进行正选择分析

发表于由
1

PAML软件中的codeml命令可以使用Maxmum Likelihood方法对密码子或氨基酸的多序列比对结果进行分析和进化参数计算 。这两种数据共用一个命令，程序运行时的参数绝大部分是公用的，但仍需要分清两种数据各自的参数意义。PAML软件中进行正选择分析的示例数据位于examples/lysozyme/目录下。

1.输入文件（1/3）：含有拓扑结构的树文件

输入 一个树文件input.trees。该文件可以是有根树，也可以是无根树，也可以带有枝长信息。codeml命令用于分析各分枝上的omega值，若是有根树的话，对于分析root节点所在的分枝，其结果可能是不对的。因此，作者推荐使用无根树进行正选择分析。无根树是3分叉的树，最少需要3个物种；而3个物种的无根数有3个分枝（2n – 3），对其进行分析，仅能得到某个物种所对应分枝的omega值，这种情况下没太大意义；而使用4个物种的无根树进行分析，可以得到两个物种共同祖先分枝的omega值。故作者推荐使用至少4~5个物种进行分析。一个树文件内容示例： 

7  1
((Hsa_Human, Hla_gibbon) #1,((Cgu/Can_colobus, Pne_langur) #1, Mmu_rhesus), (Ssc_squirrelM, Cja_marmoset));

文件内容分两行：第一行表述树中有5个物种，共计1个树，两个数值之间用任意长度的空分割。此外，Newick格式的树尾部一定要有分号，没有的话程序可能不能正常运行。

2. 输入文件（2/3）：密码子序列比对结果文件 input.txt

7  390
Hsa_Human        AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Hla_gibbon       AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Cgu/Can_colobus  AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Pne_langur       AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Mmu_rhesus       AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT...
Ssc_squirrelM    AAGGTCTTCGAAAGGTGTGAGTTGGCCAGAACT...
Cja_marmoset     AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...

PAML要求输入的Phylip格式，其物种名和后面的序列之间至少间隔两个空格（是为了允许物种名的属名和种名之间有一个空格）。

3. 输入文件（3/3）：codeml命令的配置文件 codeml.ctl

软件难点在于配置文件的理解，常用示例如下：

*输入输出参数：
      seqfile = input.txt    *设置输入的多序列比对文件路径。
     treefile = input.trees  *设置输入的树文件路径。
      outfile = mlc          *设置输出文件路径。
        noisy = 9            *设置输出到屏幕上的信息量等级：0，1，2，3，9。
      verbose = 0            *设置输出到结果文件中的信息量等级：0，精简模式结果（推荐）；1，输出详细信息，包含碱基序列；2，输出更多信息。
        getSE = 0            *设置是否计算并获得各参数的标准误：0，不需要；1，需要。
 RateAncestor = 0            *设置是否计算序列中每个位点的替换率：0，不需要；1，需要。当设置成1时，会在结果文件rates中给出各个位点的碱基替换率；同时也会进行祖先序列的重构，在结果文件rst中有体现。

*数据使用说明参数：
      runmode = 0      *设置程序获取进化树拓扑结构的方法：0，直接从输入文件中获取进化树拓扑结构（推荐）；1，程序以输入文件中的多分叉树作为起始树，并使用星状分解法搜索最佳树；2，程序直接以星状树作为起始树，并使用星状分解法搜索最佳树；3，使用逐步添加法搜索最佳树；4，使用简约法获取起始树，再进行邻近分枝交换寻找最优树；5，以输入文件中的树作为起始树，再进行邻近分枝交换寻找最优树；-2，对密码子序列进行两两比较并使用ML方法计算DnDs，或对蛋白序列进行两两比较计算ML距离，而不进行其它参数（枝长和omega等）的计算。
* fix_blength = 1      *设置程序处理输入树文件中枝长数据的方法：0，忽略输入树文件中的枝长信息；-1，使用一个随机起始进行进行计算；1，以输入的枝长信息作为初始值进行ML迭代分析，此时需要注意输入的枝长信息是碱基替换率，而不是碱基替换数；2，不需要使用ML迭代计算枝长，直接使用输入的枝长信息，需要注意，若枝长信息和序列信息不吻合可能导致程序崩溃；3，让ML计算出的枝长和输入的枝长呈正比。
      seqtype = 1      *设置输入的多序列比对数据的类型：1，密码子数据；2，氨基酸数据；3，输入数据虽然为密码子序列，但先转换为氨基酸序列后再进行分析。
    CodonFreq = 2      *设置密码子频率的计算算法：0，表示除三种终止密码子频率为0，其余密码值频率全部为1/61；1，程序分别计算三个密码子位点的四种碱基频率，再算四种碱基频率在三个位点上的算术平均值，计算61种密码子频率时，使用该均值进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；2，程序分别计算三个密码子位点的四种碱基频率，计算61种密码子频率时，使用相应位点的碱基频率进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；3，直接使用观测到的各密码子的总的频数/所有密码值的总数，得到所有密码子的频率。一般选择第三种方法，设置CodonFre的值为2。第一种方法让所有密码子频率相等，不符合实际轻卡；第二种方法没有考虑到三种位点各碱基频率的差异，得到的密码子频率不准确；第三种方法能比较准确计算出各密码子的频率；第四种方法由输入数据得到真实的各密码子频率，但有时候有些非终止密码子的频率为0，可能不利于后续的计算。
    cleandata = 1      *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进行数据分析：0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要。
*       ndata = 1      *设置输入的多序列比对的数据个数。
        clock = 0      *设置进化树中各分支上的变异速率是否一致，服从分子钟理论：0，变异速率不一致，不服从分子钟理论（当输入数据中物种相差较远时，各分支变异速率不一致）；1，所有分枝具有相同的变异速率，全局上服从分子钟理论；2，进化树局部符和分子钟理论，程序认为除了指定的分支具有不同的进化速率，其它分枝都具有相同的进化速率，这要求输入的进化树信息中使用#来指定分支；3，对多基因数据进行联合分析。
        Mgene = 0      *设置是否有多个基因的多序列比对信息输入，以及多各基因之间的参数是否一致：0，输入的多序列比对文件中仅包含一个基因时，或多个基因具有相同的Kappa和Pi参数；1，输入文件包含多个基因，这些基因之间是相互独立的（这些基因之间具有不同的Kappa和Pi值，且其进化树的枝长也不相关）；2，输入文件包含多个基因，这些基因具有相同的Kappa值，不同的Pi值；3，输入文件包含多个基因，这些基因具有相同的Pi值，不同的Kappa值；4，输入文件包含多个基因，这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时，多个基因的进化树枝长虽然长度不一样，但是呈正比关系。这点和参数值等于1是不一样的。
        icode = 0      *设置遗传密码。其值1-10和NCBI的1-11遗传密码规则对应：0，表示通用的遗传密码。
   Small_Diff = .5e-6  *设置一个很小的值，一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果，需要该参数值对结果没什么影响。

*位点替换模型参数：
        model = 1            *若输入数据是密码子序列，该参数用于设置branch models，即进化树各分枝的omega值的分布：0，进化树上所有分枝的omega值一致；1，对每个分枝单独进行omega计算；2，设置多类omega值，根据树文件中对分枝的编号信息来确定类别，具有相同编号的分枝具有相同的omega值，没有编号的分枝具有相同的omega值，程序分别计算各编号和没有编号的omega值。
                             *若输入数据是蛋白序列，或数据是密码子序列且seqtype值是3时，该参数用于设置氨基酸替换模型：0，Poisson；1，氨基酸替换率和氨基酸的观测频率成正比；2，从aaRatefile参数指定的文件路径中读取氨基酸替换率信息，这些信息是根据经验获得，并记录到后缀为.dat的配置文件中。这些经验模型(Empirical Models)文件位于PAML软件安装目录中的dat目录下，例如，Dayhoff(dayhoff.dat)、WAG(wag.dat)、LG(lg.dat)、mtMAN(mtman.dat)和mtREV24(mtREV24.dat)等；
   aaRatefile = dat/wag.dat  *当对蛋白数据进行分析，且model = 2时，该参数生效，用于设置氨基酸替换模型。
       aaDist = 0            *设置氨基酸之间的距离。
      NSsites = 0            *输入数据时密码子序列时生效，用于设置site model，即序列各位点的omega值的分布：0，所有位点具有相同的omega值；1，各位点上的omega值小于1或等于1（服从中性进化neutral）；2，各位点上的omega值小于1、等于1或大于1（选择性进化selection）；3，discrete；4，freq；5:gamma；6，2gamma；7，beta；8，beta&w；9，beta&gamma；10，beta&gamma+1；11，beta&normal>1；12，0&2normal>1；13，3normal>0。
                             *可以一次输入多个模型进行计算并比较，其结果输出的rst文件中。
    fix_alpha = 1            *序列中不同的位点具有不同的碱基替换率，服从discrete-gamma分布，该模型通过alpha（shape参数）和ncatG参数控制。该参数设置是否给定一个alpha值：0，使用ML方法对alpha值进行计算；1，使用下一个参数设置一个固定的alpha值。
	                     *对于密码子序列，当NSsites参数值不为0或model不为0时，推荐设置fix_alpha = 1且alpha = 0，即不设置alpha值，认为位点间的变异速率一致，否则程序报错。若设置了alpha值，则程序认为不同密码子位点的变异速率不均匀，且同时所有位点的omega值一致，当然各分枝的omega值也会一致，这时要求NSsites和model参数值都设置为0（这一般不是我们需要的分析，它不能进行正选择分析了）。
        alpha = 0            *设置一个固定的alpha值或初始alpha值（推荐设置为0.5）。该值小于1，表示只有少数热点位置的替换率较高；该值越小，表示位点替换率在各位点上越不均匀；若设置fix_alpha = 1且alpha = 0则表示所有位点的替换率是恒定一致的。
       Malpha = 0            *当输入的多序列比对结果中有多基因时，设置这些基因间的alpha值是否相等：0，分别对每个基因单独计算alpha值；1，所有基因的alpha值保持一致。
        ncatG = 5            *序列中不同位点的变异速率服从GAMMA分布的，ncatG是其一个参数，一般设置为5，4，8或10，且序列条数越多，该值设置越大。
                             *对于密码子序列，当NSites设置为3时，ncatG设置为3；当NSites设置为4时，ncatG设置为5；当NSites值设置>=5时，ncatG值设置为10。
    fix_kappa = 0            *设置是否给定一个Kappa值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的Kappa值。
        kappa = 2            *设置一个固定的Kappa值，或一个初始的Kappa值。
    fix_omega = 0            *设置是否给定一个omega值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的omega值。 
        omega = .4           *设置一个固定的omega值，或一个初始的omega值。
       method = 0            *设置评估枝长的ML迭代算法：0，使用PAML的老算法同时计算所有枝长，在clock = 0下有效；1，PAML新加入的算法，一次对一个枝长进行计算，该算法仅在clock参数值为1，2或3下工作。

4. 运行codeml进行分析

运行程序的命令：

codeml codeml.ctl

5. 正选择基因分析

对基因进行正选择分析，可以按如下步骤进行：

(1) 收集该基因在多个物种中的CDS和protein序列（可以使用orthoMCL分析结果）。
(2) 对该基因的蛋白序列进行多序列比对，再根据CDS序列转换成Codon序列比对结果。这个步骤需要自己编写一些程序来实现。
(3) 根据Codon序列的比对结果构建系统发育树。推荐使用RAxML软件使用ML算法对所有的单拷贝同源基因进行物种树构建。然后，输入物种树信息、Codon多序列比对信息，使用codeml进行omega计算和选择压模型检验。
(4) 快速过滤非正选择基因：设置runmode = -2运行codeml命令对两两序列比较，使用ML方法计算dNdS，若存在omega值大于1，则认为该基因可能属于正选择基因，进行后续分析。
(5) 正选择基因的初步鉴定方法：设置model = 0、NSsites = 1 2运行codeml命令（site models）分别对正选择模型(M2a)和近中性进化模型(M1a)进行检测，若两种模型的似然值相差较大，通过自由度为2的卡方检验，可以确定该基因是否为正选择基因。从程序结果中（例如：lnL(ntime: 11  np: 16):   -870.867266      +0.000000）可以找到M2a模型的似然值l1和M1a模型的似然值l0，再计算2Δl = 2(l1-l0)，通过命令“chi2 2 2Δl”根据卡方检验算出p值。此外，也可以设置NSsites = 7 8进行计算，则表示分别对M8(beta & omega)和M7(beta)模型进行检测，和前者类似，若两种模型的似然值相差较大，通过自由度为2的卡方检验，可以确定该基因是否为正选择基因。根据PAML作者所说，M2a和M1a的比较，比M7-M8的比较更严格。即若想得到更多的正选择基因，可以使用M7_VS_M8分析。此外，根据M2和M8模型的BEB（Bayes empirical Bayes）分析结果，还可以得到在多序列比对结果中第一条序列上的正选择位点。根据PAML作者所说，这种site models之间的比较，能用于检测是否在序列上不同的位点具有不同的omega值，而不是用于正选择检测（We suggest that The M0-M3 comparison should be used as a test of variable ω among sites rather than a test of positive selection）。
(6) 指定分化枝上的正选择基因鉴定方法：经过上一步初步鉴定后，设置model = 2、NSsites = 2、fix_omega = 0、omega = 2.0运行codeml命令（branch-site model A）；再设置model = 2、NSsites = 2、fix_omega = 1、omega = 1运行codeml命令（modified branch-site model A / null model）。进行这两种模型分析时，要求输入的树文件中对目标分化枝进行标注。对这两种模型进行LRT分析，计算2Δl = 2(l1-l0)，注意是前者的似然值减后者（null model）的似然值；再使用自由度为1的卡方检验，通过命令“chi2 1 2Δl”计算出的值再除以2，即得到p值。
(7) 目标分化枝上的快速进化基因（Rapidly evolving gene）鉴定方法：设置model = 0、NSsites = 0运行codeml命令，再设置model = 2、NSsites = 0运行codeml命令（这时要求输入的树文件中对目标分化枝进行标注）。然后比较两种运行模式下的似然值，通过自由度为1的卡方检验，可以鉴定该基因在目标分化枝的omega值和背景差异显著。

使用PAML软件利用核酸序列进行系统发育参数分析

发表于由
回复

PAML软件中的baseml命令可以利用Maxmum Likelihood方法计算系统发育中的相关参数。basml命令能以分支树拓扑结构信息和多序列比对信息作为输入，在参数配置文件中设置一些控制程序运行的参数，来对一些进化过程中的参数进行分析。这些参数主要有：

1. 序列中各碱基频率。
2. 两两物种之间的距离。
3. 系统发育树的枝长信息。
4. 碱基替换矩阵。

PAML软件的baseml示例数据位于软件安装目录下。

1. 输入文件（1/3）：含有拓扑结构的树文件 input.trees

输入一个树文件input.trees。该文件可以是有根树，也可以是无根树，也可以带有枝长信息。例如：

5  1
(((Human:0.1, Chimpanzee:0.2):0.8, Gorilla:0.3):0.7, Orangutan:0.4, Gibbon:0.5);

文件内容分两行：第一行表述树中有5个物种，共计1个树，两个数值之间用任意长度的空分割。此外，Newick格式的树尾部一定要有分号，没有的话程序可能不能正常运行。

2. 输入文件（2/3）：核酸序列比对结果 input.txt

5   895

Human       AAGCTTCACCGGCGCAGTCATTCTCATAATCGCCCACGGACTT...
Chimpanzee  AAGCTTCACCGGCGCAATTATCCTCATAATCGCCCACGGACTT...
Gorilla     AAGCTTCACCGGCGCAGTTGTTCTTATAATTGCCCACGGACTT...
Orangutan   AAGCTTCACCGGCGCAACCACCCTCATGATTGCCCATGGACTC...
Gibbon      AAGCTTTACAGGTGCAACCGTCCTCATAATCGCCCACGGACTA...

PAML要求输入的Phylip格式，其物种名和后面的序列之间至少间隔两个空格（是为了允许物种名的属名和种名之间有一个空格）。

3. 输入文件（3/3）：baseml命令的配置文件baseml.ctl

软件难点在于配置文件的理解，常用示例如下：

*输入输出参数：
      seqfile = input.txt   *输入核酸多序列比对结果文件路径
     treefile = input.trees *输入系统发育树文件路径
      outfile = mlb         *设置输出文件路径
        noisy = 2           *设置屏幕上输出信息模式，值越大则输出越多信息：0，1，2，3。
      verbose = 0           *结果文件中的输出信息方式：1，输出详细信息，包含碱基序列；0，精简模式结果（推荐）。
        getSE = 0           *设置是否计算并获得各参数的标准误：0，不需要；1，需要。
 RateAncestor = 0           *设置是否计算序列中每个位点的碱基替换率：0，不需要；1，需要。当设置成1时，会在结果文件rates中给出各个位点的碱基替换率；同时也会进行祖先序列的重构，在结果文件rst中有体现。

*数据使用说明参数：
      runmode = 0     *设置程序获取进化树拓扑结构的方法：0，直接从输入文件中获取进化树拓扑结构（推荐）；1，程序以输入文件中的多分叉树作为起始树，并使用星状分解法搜索最佳树；2，程序直接以星状树作为起始树，并使用星状分解法搜索最佳树；3，使用逐步添加法搜索最佳树；4，使用简约法获取起始树，再进行邻近分枝交换寻找最优树；5，以输入文件中的树作为起始树，再进行邻近分枝交换寻找最优树。 
* fix_blength = 1     *设置程序处理输入树文件中枝长数据的方法：0，忽略输入树文件中的枝长信息；-1，使用一个随机起始进行进行计算；1，以输入的枝长信息作为初始值进行ML迭代分析，此时需要注意输入的枝长信息是碱基替换率，而不是碱基替换数；2，不需要使用ML迭代计算枝长，直接使用输入的枝长信息，需要注意，若枝长信息和序列信息不吻合可能导致程序崩溃；3，让ML计算出的枝长和输入的枝长呈正比。
    cleandata = 1     *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进行数据分析：0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要。
        clock = 0     *设置进化树中各分支上的变异速率是否一致，服从分子钟理论：0，变异速率不一致，不服从分子钟理论（当输入数据中物种相差较远时，各分支变异速率不一致）；1，所有分枝具有相同的变异速率，全局上服从分子钟理论；2，进化树局部符和分子钟理论，程序认为除了指定的分支具有不同的进化速率，其它分枝都具有相同的进化速率，这要求输入的进化树信息中使用#来指定分支；3，对多基因数据进行联合分析。
        Mgene = 0     *设置是否有多个基因的多序列比对信息输入，以及多各基因之间的参数是否一致：0，输入的多序列比对文件中仅包含一个基因时，或多个基因具有相同的Kappa和Pi参数；1，输入文件包含多个基因，这些基因之间是相互独立的（这些基因之间具有不同的Kappa和Pi值，且其进化树的枝长也不相关）；2，输入文件包含多个基因，这些基因具有相同的Kappa值，不同的Pi值；3，输入文件包含多个基因，这些基因具有相同的Pi值，不同的Kappa值；4，输入文件包含多个基因，这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时，多个基因的进化树枝长虽然长度不一样，但是呈正比关系。这点和参数值等于1是不一样的。
*       ndata = 1     *设置输入的多序列比对文件中的数据个数。
*       icode = 0     *设置遗传密码。当设置参数RateAncestor = 1时生效，有利于根据密码子重构编码蛋白的祖先DNA序列。
   Small_Diff = 7e-6  *设置一个很小的值，一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果，需要该参数值对结果没什么影响。
	  
*位点替换模型参数：
        nhomo = 1     *设置碱基替换模型中的频率参数计算方法：0，使用输入数据中的观测到的所有序列的平均碱基频率；1，使用ML迭代方法计算碱基频率Pi(T)、Pi(C)、Pi(A)和Pi(G)，适合model参数值为F81、F84、HKY85、T92、T93或REV模型；2，进化树所有的分枝具有相同的Kappa值，适合K80、F84和HKY85这三个含有Kappa参数的模型；3，外部节点所在的分枝使用一组碱基频率，所有内部节点所在分枝使用一组碱基频率，root节点所在分枝使用一组碱基频率，也称为N1模型，需要输入有根树，适合model参数值为F84、HKY85和T95模型；4，root节点所在分枝使用一组碱基频率，所有其它分枝使用一组碱基频率，也称为N2模型，需要输入有根树，适合model参数值为F84、HKY85、T95和GTR/REV模型；5，让用户自己定义各组分枝使用的碱基频率，需要在进化树使用#对各分枝进行编号。
        model = 7     *选取碱基替换模型方法：0，JC69；1，K80；2，F81；3，F84；4，HKY85；5，T92；6，TN93；7，REV；8，UNREST；9，REVu；10，UNRESTu。
                      *一般为了得到较为准确的枝长信息，至少选择4，使用HKY85碱基替换模型，程序会使用ML方法计算Kappa值（转换率/颠换率比值）。
                      *推荐选择7，也称为GTR，程序会使用ML计算碱基替换模型。
    fix_alpha = 0     *序列中不同的位点具有不同的碱基替换率，服从discrete-gamma分布，该模型通过alpha（shape参数）和ncatG参数控制。该参数设置是否给定一个alpha值：0，使用ML方法对alpha值进行计算；1，使用下一个参数设置一个固定的alpha值。
        alpha = 0.5   *设置一个固定的alpha值或初始alpha值（推荐设置为0.5）。该值小于1，表示只有少数热点位置的碱基替换率较高；该值越小，表示碱基替换率在各位点上越不均匀；若设置fix_alpha = 1且alpha = 0则表示所有位点的碱基替换率时恒定一致的。
       Malpha = 0     *当输入数据时多基因时，设置这些基因之间的alpha是否不一致：1，不同基因之间具有不同的alpha值；0，所有基因都具有相同的alpha值。
        ncatG = 5     *序列中不同位点的变异速率服从GAMMA分布的，ncatG是其一个参数，一般设置为5，4，8或10，且序列条数越多，该值设置越大。
        nparK = 0     *设置rate-class models：1，不同位点的变异速率是独立的，限定ncatG各类的概率都相等；2，不同位点的变异速率是独立的，不限定ncatG各类的概率都相等；3，邻近位点的变异速率是Markov相关的，限定ncatG各类的概率都相等；4，邻近位点的变异速率是Markov相关的，不限定ncatG各类的概率都相等。此外，设置该参数值不为0，比较消耗计算并且要求nhomo = 0，否则程序报错。
    fix_kappa = 0     *设置是否给定一个Kappa值（该参数在model设置为K80、F84或HKY85时有效）：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的Kappa值。对于JC69和F81，该参数没有效果；对于TN93和REV，则分别对应2个和5个rate参数。
        kappa = 5     *设置一个固定的Kappa值，或一个初始的Kappa值。
       method = 0     *设置评估枝长的ML迭代算法：0，使用PAML的老算法同时计算所有枝长，在clock = 0下有效；1，PAML新加入的算法，一次对一个枝长进行计算，该算法仅在clock参数值为1，2或3下工作。

4. 运行baseml命令对核酸序列进行进化参数计算

程序运行的命令：

baseml baseml.ctl

5. 结果文件解释

程序的主要结果文件是mlb，其主要内容有：

1. 四种碱基观测到的频率及其平均值。
2. 各序列之间的距离矩阵。
3. 各分枝的长度、总的枝长、含有枝长信息的系统发育树。
4. ML计算出来的碱基频率。
5. 碱基替换矩阵。
6. Gamma分布的alpha值。
7. 若使用K81、F84和KEY85模型，则给出Kappa参数值；若选择TN93或REV模型，则分别给出2个和5个rate参数值。

6. baseml命令的意义

PAML软件不擅长搜索进化树的最优拓扑结果，它仅适合于对小数据（序列条数少于10）进行最优树的搜索。一般使用其它ML软件进行系统发育树分析，就可以得到包含枝长信息的最优ML树。再将该树信息输入给basml命令进行分析，又得到一个包含枝长信息的进化树，此外还得到一些没有什么意义的参数值。很多人会疑惑basml命令又有什么用呢？

在我看来，baseml命令的用处：(1) 通过这个软件的使用和练习，可以更加深入了解进化中各参数的意义和系统发育树枝长计算的原理。(2) 再进行物种分歧时间计算时，是需要先通过baseml命令计算相关参数和枝长的。使用mcmctree命令进行物种分歧时间计算时，第一步会调用baseml进行计算。(3) 若使用数百上千个基因构建物种树，可以使用baseml命令对各个基因分别计算进化树枝长，然后通过比较枝长来去除掉和物种树枝长差异较大的基因，再进一步得到更准确的物种树。

GAMMA分布

发表于由
回复

1. 序列不同位点的进化速率服从GAMMA分布

在做进化分析时，很多参数都服从GAMMA distribution（伽马分布）。最典型的参数是序列中各个位点的进化进化速率服从GAMMA分布。序列中的位点数量很多，一般情况下不同位点的进化速率差别较大，比如：密码子第3位碱基进化速率快，第1位和2位进化较慢。

GAMMA分布由alpha(形状)和beta(尺度)两个参数决定。其均值 = alpha / beta，方差 = alpha / (beta的平方)。当alpha > 1时，分布呈钟形，表示大部分位点的进化速率位于均值附近，趋于一致；当alpha <= 1时，分布为一个高度倾斜的L形，表示大部分位点的进化速率非常低，只有少部分位点属于进化的热点。一般情况下，对于真实的多序列比对数据，不同位点其碱基替换率是不一致的，推荐设置alpha = 0.5。

2. 其它服从GAMMA分布的进化参数

(1) rate: 所有位点的平均[碱基/密码子/氨基酸]替换率。

(2) sigma2: 所有位点进化速率取对数后的方差。

(3) kappa: 转换/颠换的比率。当氨基酸替换模型选择HKY时，则需要计算kappa参数。

(4) alpha: 不同位点变异速率不一样从而使用GAMMA分布时，对其形状参数alpha使用GAMMA分布进行计算。

使用PAML进行分歧时间计算

发表于由
11

PAML软件中的mcmctree命令可以使用Bayesian方法估算物种分歧时间。对程序输入带有校准点信息的有根树、多序列比对结果，即可得到进化树各内部节点95%置信区间分歧时间信息。PAML软件也能利用带有枝长信息的树，快速进行分歧时间计算。PAML软件中的练习数据位于文件夹 install_path/paml4.9i/examples/DatingSoftBound目录下。详细说明文档： http://abacus.gene.ucl.ac.uk/software/MCMCtree.Tutorials.pdf

1. 输入文件（1/3）：带有校准点的有根树 input.trees

输入一个带有校准点的有根树文件，示例如下：

7 1
((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

文件内容分两行：第一行表述树中有7个物种，共计1个树，两个数值之间用空分割；第二行则是Newick格式树信息，其中包含有校准点信息。校准点信息一般指95%HPD（Highest Posterior Density）对应的置信区间；校准点单位是100MYA（软件说明文档中使用该单位，也推荐使用该单位，若使用其它单位，后续配置文件中的相关参数也需要对应修改）。此外，Newick格式的树尾部一定要有分号，没有的话程序可能不能正常运行。

2. 输入文件（2/3）：密码子在3个位点的多序列比对结果（Phylip格式） input.txt

7  3331
human               AGACTGTTAGCAGCGCCGGGCAACTCCCTACTCAA...
chimpanzee          AGACTGTTGGCAACGTCGGGCAACTCCCCGCCCAA...
bonobo              AGACTGTTGGCAACGCCGGGCAACTCCCCGCCCAA...
gorilla             AGATTGTTAGCAACGTCGGGTAACCCCCCACTCAA...
orangutan           AGGCTACTAACAGCGCCGGACGATTCCTCGCCTAA...
sumatran            AGACTACTAACAGCGCCGGGCGATTCCTCACCCAA...
gibbon              AGACTATTGACAACGTCGGGCAACTCTCTACTCAA...
7  3331
human               AAATTCCTTCCCTTGTCCCTTTTTTCCTTTCATTA...
chimpanzee          AAATTCCTCCCCTTGTCCCTTTTTTCCTTTCATTA...
bonobo              AAATTCCTCCCCTTGTCCCTTTTTTCCTTTCATTA...
gorilla             AAATTCCTTCCCTTGTCCCTTTTTTCCTTTCATTA...
orangutan           AAATTCCTCCCCTTGTCCCTTTTTTCCTTTCATTA...
sumatran            AAGTTCCTTCCCTTGTCCCTTTTTTCCTTTCATTA...
gibbon              AAATTCCTCCCCTTGTCCCTCTTTTCCTTTCATTA...
7  3331
human               CATGCTACTCCACACACCAAGCTATCTAGCCTCCC...
chimpanzee          CATACTACTCCACACACCAAACTACCTAGCCTCCC...
bonobo              CATGCTACTCCACACACCAAGCTACCTAGTCTCCC...
gorilla             CATACTACTCCACACACCAAATCATCTAGCCTCCC...
orangutan           CATACCACTCCACACCCTATACCATCCAACTTCCC...
sumatran            CATATCACTCCAAACCCCAAACCATCCAGCCTCCC...
gibbon              CATACTACTCCATACACCAAATTATCCAACTCCCC...

PAML要求输入的Phylip格式，其物种名和后面的序列之间至少间隔两个空格（是为了允许物种名的属名和种名之间有一个空格）。

当然，也可以输入一个多序列比对结果进行分析。在软件说明文档中，是密码子三个位点分别的多序列比对结果，程序可以分别计算密码子各个位点的变异速率。

此外，也可以输入密码子或氨基酸序列比对信息，则需要再配置文件中修改相应的参数来表明数据类型。若使用氨基酸序列来进行分析，由于mcmctree命令不能选择较好的氨基酸替换模型进行分析，需要自己手动运行codeml进行分析后，在生成中间文件用于运行mcmctree，比较麻烦。因此，推荐按上述方法使用密码子三位点的碱基序列作为输入信息进行PAML分析。

3. 输入文件（3/3）：mcmctree命令的配置文件 mcmctree.ctl

软件难点在于配置文件的理解，示例如下：

*输入输出参数：
        seed = -1           *设置一个随机数来运行程序。若设置为-1，表示使用系统当前时间为随机数。
     seqfile = input.txt    *设置输入的多序列比对文件路径
    treefile = input.trees  *设置输入的带校准点信息有根树的文件路径 
    mcmcfile = mcmc.txt     *设置输出的mcmc信息文本文件，可以用Tracer软件查看
     outfile = out.txt      *设置输出文件路径，该文件中记录了超度量树和进化速率等参数信息。

*数据使用说明参数：
       ndata = 3         *设置输入的多序列比对的数据个数。这里指密码子3个位置的数据。
     seqtype = 0         *设置多序列比对的数据类型：0，表示核酸数据；1，表示密码子比对数据；2，表示氨基酸数据。根据不同的数据类型，程序调用相应的baseml或codeml进行相应的参数计算。
     usedata = 1         *设置是否利用多序列比对的数据：0，表示不使用多序列比对数据，则不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但是其分歧时间结果是有问题的；1，表示使用多序列比对数据进行likelihood计算，正常进行MCMC，是一般使用的参数; 2，进行正常的approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。此外，由于程序BUG，当设置usedata = 2时，一定要在改行参数后加 *，否则程序报错 Error: file name empty.. 3，程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。由于mcmctree调用baseml/codeml进行计算的参数设置可能不太好（特别时对蛋白序列进行计算时），推荐自己修该软件自动生成的baseml/codeml配置文件，然后再手动运行baseml/codeml命令，再整合其结果文件为out.BV文件。
   cleandata = 0         *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进行数据分析：0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要。 
       clock = 2         *设置分子种方法：1，global clock方法，表示所有分支进化速率一致；2，independent rates方法，各分支的进化速率独立且进化速率的对数log(r)符合正态分布; 3，correlated rates方法，和方法2类似，但是log(r)的方差和时间t相关。
*    TipDate = 1 100     *当外部节点由取样时间时使用该参数进行设置，同时该参数也设置了时间单位。具体数据示例请见examples/TipData文件夹。
     RootAge = '<1.0'    *设置root节点的分歧时间，一般设置一个最大值。
 
*位点替换模型参数：
       model = 4         *设置碱基替换模型：0，JC69；1，K80；2，F81；3，F84；4，HKY85。当设置usedata = 1时，不能使用其它碱基替换模型。若需要选择其它碱基替换模型，则考虑使用usedata = 3参数，就可以设置model参数值为其它的碱基替换模型。此时，程序会将数据进行分割，例如，ndata = 3，则程序分别对3个数据进行分析，生成baseml的输入文件，并调用baseml进行分析（若数据量较大，分析较慢，推荐按ctrl + c 强行终止，程序则终止baseml的运行，进行下一个数据的输入文件生成并运行下个数据的baseml命令，再按 ctrl + c 强行终止，这样可以让程序生成3各数据的输入文件和baseml配置文件，然后手动并行化运行，加快时间）；每个分析结果中会得到rst2结果文件；可以在此处手动修改tmp0001.ctl配置文件，例如，修改其model参数值，推荐修改outfile参数值，从而可以手动并行化运行baseml命令；然后将所有数据的rst2结果使用cat命令合并成文件in.BV，用于下一步设置参数usedata = 2，进行MCMC分析。
       alpha = 0.5       *核酸序列中不同位点，其进化速率不一致，其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。若该参数值设置为0，则表示所有位点的进化速率一致。此外，当userdata = 2时，alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时，不会利用到多序列比对的数据。
       ncatG = 5         *设置离散型GAMMA分布的categories值。
     BDparas = 1 1 0.1   *设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变，则需要相应修改出生率和死亡率的值。例如，时间单位由100Myr变换为1Myr，则要设置成".01 .01 0.1"。
 kappa_gamma = 6 2       *设置kappa（转换/颠换比率）的GAMMA分布参数。
 alpha_gamma = 1 1       *设置GAMMA形状参数alpha的GAMMA分布参数。 

*进化速率参数：
 rgene_gamma = 2 20 1    *设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数：alpha=2、beta=20、初始平均替换率为每100million年（取决于输入有根树文件中的时间单位）1个替换。若时间单位由100Myr变换为1Myr，则要设置成"2 2000 1"。总体上的平均进化速率为：2 / 20 = 0.1 个替换 / 每100Myr，即每个位点每年的替换数为 1e-9。
sigma2_gamma = 1 10 1    *设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。当clock参数值为1时，表示使用全局的进化速率，各分枝的进化速率没有差异，即方差为0，该参数无效；当clock参数值为2时，若修改了时间单位，该参数值不需要改变；当clock参数值为3时，若修改了时间单位，该参数值需要改变。
    finetune = 1: .1 .1 .1 .1 .1 .1    *冒号前的值设置是否自动进行finetune，一般设置成1，然程序自动进行优化分析；冒号后面设置各个参数的步进值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置，该参数不像以前版本那么重要了，可能以后会取消该参数。

*MCMC参数：
       print = 1         *设置打印mcmc的取样信息：0，不打印mcmc结果；1，打印除了分支进化速率的其它信息（各内部节点分歧时间、平均进化速率、sigma2值）；2，打印所有信息。 
      burnin = 2000      *将前2000次迭代burnin后，再进行取样（即打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等）。
    sampfreq = 10        *每10次迭代则取样一次。
     nsample = 20000     *当取样次数达到该次数时，则取样结束（程序也将运行结束）。

a4. 运行mcmctree命令使用进行分歧时间计算

程序运行命令：

mcmctree mcmctree.ctl

5. 结果文件解释

程序在运行过程中，会在屏幕生生成一些信息。比较耗时间的步骤主要在于取样的百分比进度：

第一列：取样的百分比进度。
第2~6列：参数的接受比例。一般，其值应该在30%左右。20~40%是很好的结果，15~70%是可以接受的范围。若这些值在开始时变动较大，则表示burnin数设置太小。
第7~x列：各内部节点的平均分歧时间，第7列则是root节点的平均分歧时间。若有y个物种，则总共有y-1个内部节点，从第7列开始后的y-1列对应这些内部节点。
倒数第3~x列：r_left值。若输入3各多序列比对结果，则有3列。x列的前一列是中划线 - 。
倒数第1~2列：likelihood值和时间消耗。

屏幕信息最后，给出各个内部节点的分歧时间(t)、平均变异速率(mu)、变异速率方差(sigma2)和r_left的Posterior信息：均值(mean)、95%双侧置信区间(95% Equal-tail CI)和95% HPD置信区间(95% HPD CI)等信息。此外，倒数第二列给出了各个内部节点的Posterior mean time信息，可以用于收敛分析。

在当前目录下，生成几个主要结果：

FigTree.tre    生成含有分歧时间的超度量树文件
mcmc.txt       MCMC取样信息，包含各内部节点分歧时间、平均进化速率、sigma2值等信息，可以在Tracer软件中打开。通过查看各参数的ESS值，若ESS值大于200，则从一定程度上表示MCMC过程能收敛，结果可靠。
out.txt        包含由较多信息的结果文件。例如，各碱基频率、节点命名信息、各节点分歧时间、进化速率和进化树等。

6. 使用approximate likelihood方法更快速地计算分歧时间

按照以上方法进行MCMC计算非常消耗计算。为了加快计算速度，可以将分析分成两个步骤：（1）首先，使用Maximum Likelihood方法计算枝长和Hessian信息；（2）再使用MCMC方法计算分歧时间。

运行步骤也相应分成两步：（1）设置配置文件中参数 usedata = 3，然后和上面同样运行mcmctree命令，软件会对各个多序列比对结果调用baseml/codeml命令进行分歧时间计算，生成结果文件rst2，将各个多序列比对结果文件内容直接cat（linux系统的常用命令）合并成out.BV文件；（2）然后将out.BV文件重命名为in.BV，并将配置文件中参数 usedata = 2，再次运行同样的mcmctree命令，程序会运行MCMC分析并得到结果。使用该方法，其MCMC速度快很多。

若输入数据是protein序列，则mcmctree调用codeml命令进行分析时的model参数选择是不正确的，需要自己手动修改codeml命令的参数配置文件后再运行codeml命令，并将其结果文件rst2的内容合并到out.BV文件中。若想使用更复杂的模型，如GTR进行mcmctree分析，则按同样方法手动修改配置文件并允许baseml命令，最后将rst2的内容整合到out.BV文件中。此外，可以将分别来自核酸(非编码RNA)和蛋白的rst2内容同时整合到out.BV中进行计算。

7. 使用infinitesites输入含有枝长信息的系统进化树快速进行分歧时间计算

infinitesites程序进行分歧时间计算时，其假设前提为多序列比对的序列长度为无限长。相比于正常的mcmctree命令，要求额外多输入一个名为 FixedDsClock23.txt 的文件。该文件中存放了相应的带有枝长信息的系统发育树。该文件中系统发育树的拓扑结构要和input.trees一致。推荐使用baseml命令输入input.trees和多序列比对结果计算枝长。

FixedDsClock23.txt 文件内容示例：

7
((((human: 0.029043, (chimpanzee: 0.014557, bonobo: 0.010908): 0.016729): 0.015344, gorilla: 0.033888): 0.033816, (orangutan: 0.026872, sumatran: 0.022437): 0.069648): 0.073309, gibbon: 0.024637);
((((human: 0.012463, (chimpanzee: 0.002782, bonobo: 0.003835): 0.003331): 0.004490, gorilla: 0.014278): 0.006308, (orangutan: 0.010818, sumatran: 0.008845): 0.030551): 0.004363, gibbon: 0.029246);
((((human: 0.270862, (chimpanzee: 0.066698, bonobo: 0.056883): 0.124104): 0.139082, gorilla: 0.310797): 0.391342, (orangutan: 0.152555, sumatran: 0.114176): 0.696518): 0.017607, gibbon: 1.394718);

运行infinitesites命令进行分歧时间计算：

infinitesites mcmctree.ctl

使用infinitesites进行分歧时间计算时，程序要求输入多序列比对文件。虽然程序读取了序列信息，但在计算时会忽略其序列信息。

8. 一些其它注意事项

(1) 如何检测MCMC的结果是否达到收敛状态？

收敛的意思，即经过很多次迭代后，得到的MCMC树的各枝长趋于一个定值，变动非常小。于是，最直接的检测方法是：分别使用不同的Seed值进行mcmctree或infinitesites进行两次或多次分析，然后比较两个结果树是否一致，实际就是比较树文件中各内部节点的Height值（分歧时间 / Posterior time）。比较方法可以有多种: (a) 结果文件mcmc.txt中每行记录了一个取样的结果，包含各内部节点的Posterior time值，即进化树的Height值，对该mcmc.txt文件进行分析，得到每个内部节点的Posterior mean time值。然后作图 （官方说明文档第6页图示），横坐标表示第一次运行的各内部节点的分歧时间均值、纵坐标表示第二次运行的各内部节点的分歧时间均值。该散点图要表现出非常明显的对角线，才认为达到收敛。这时可以考虑计算其相关系数来判断是否符合线性。(b) 我觉得第一种官方文档给的方法比较麻烦且是否符合对角线没有明确的定义，于是就直接比较两次结果树文件中的各枝长。计算各枝长总的偏差百分比，当偏差百分比低于0.1%，则认为两次结果非常吻合，差异低于0.1%，认为达到收敛。(c) 此外，还可以使用Tracer分析mcmc.txt文件，检测其ESS值，一般认为该值高于200，则可能达到收敛。该方法可用于辅助检测。最后，若不收敛，则需要提高burnin、nsample值，重新运行程序。

(2) 如何设置burnin、sampfreq和nsample值？

(a) 一般推荐设置nsample值不小于20k(官方示例数据中设置的值)，只有当该值较大时，求得的均值才有意义。(b) sampfreq表示取样间隔，一般设置为10，100或1000。nsample 和sampfreq 的乘积表示有效迭代的次数，这些迭代越准确，次数越大，则结果越好，越趋于收敛。当然，次数越大，越消耗时间。若数据量很小，可以考虑设置sampfreq为1000；若数据量大，每次迭代耗时多，则考虑设置sampfre为10。(c) 一般最开始的迭代结果很差，需要去除(burnin)其结果，则设置burnin值。要设置足够大的burnin值，保证在程序运行时当迭代比例达到0% (即burnin结束) 时，其参数的接受比例值较好，在30%左右且随迭代次数增多时变动幅度不大。推荐设置burnin的迭代次数为有效次数的10~40%。PAML软件时先burnin后再记录有效迭代的参数值；其它软件如BEAST则记录所有的参数值后，最后汇总时burnin掉一定比例的数据。(d) 总体上，其实就是两个参数：burnin掉的迭代次数和用于计算结果的有效迭代次数。由于需要根据有效迭代数据来进行均值计算，若记录每次迭代的参数，则生成的mcmc.txt文件行数过多，文件太大，汇总时也消耗计算。于是设置没隔一定迭代轮数取样一次，这样生成的mcmc.txt文件不会太大，对最后的均值计算影响不大。(e) 我个人推荐有效迭代次数 1~10M 次，burnin 0.2~4M次。按时间消耗从少到多的参数值：

数据简单时，进行0.5M迭代次数，burnin比例40%：{ burnin 200k; sampfreq 10; nsample 50k } 
数据中等时，进行1M迭代次数，burnin比例40%：{ burnin 400k; sampfreq 10; nsample 100k } 
数据复杂时，进行5M迭代次数，burnin比例20%：{ burnin 1000k; sampfreq 10; nsample 500k }
数据巨大时，进行10M迭代次数，burnin比例20%：{ burnin 2000k; sampfreq 100; nsample 100k }

使用相应参数进行分析后，若不满意其收敛程度，则选用更多迭代次数的参数。

集群中新添加服务器CentOS7系统的准备工作

发表于由
2

在集群中新增加一台服务器，或对集群中的某台服务器重新安装CentOS7系统后。需要对该机器进行相关配置。以下按照我使用的集群实例进行设置

1. 准备网络

进入到新机器CentOS7系统的桌面，鼠标移动到右上角打开控制面板，在NetWork栏目中设置对应网口的IP信息。若打开不了NetWork设置界面，则运行命令 NetworkManager 开启网络设置界面。

IP地址：192.168.30.6
子网掩码：255.255.255.0
网关：192.168.30.4   （网关设置成提供了网络共享的Control服务器，从而可以让本服务器连接外网）
DNS: 211.69.143.1

在IP设置完毕后，修改 /etc/sysconfig/network-scripts/ 目录中对应网口的配置文件，例如，ifcfg-enp4s0f0，使参数值设置为： ONBOOT=yes ，以利于开机启动对应网口。

2. 准备主机名

修改配置文件 /etc/sysconfig/network 内容来设置新机器的主机名：

NETWORKING=yes
HOSTNAME=node6

再修改 /etc/hosts ，设置集群中各机器的主机名，以利于机器间相互识别。

127.0.0.1   localhost localhost.localdomain localhost4 localhost4.localdomain4
::1         localhost localhost.localdomain localhost6 localhost6.localdomain6
192.168.30.1    node1
192.168.30.2    node2
192.168.30.3    node3
192.168.30.4    control
192.168.30.5    node5
192.168.30.6    node6

3. 准备用户

在control服务器的root用户下配置无密码登录到新服务器node6的root用户中：

ssh-keygen -t dsa -P '' -f ~/.ssh/id_dsa
ssh-copy-id -i ~/.ssh/id_dsa.pub root@node6

在control服务器中，使用root权限将用户信息配置文件上传到新的node6服务器中：

scp /etc/passwd /etc/shadow /etc/group node6:/

在新服务器node6中使用root权限修改用户信息：

cat /passwd > /etc/passwd
cat /shadow > /etc/shadow
cat /group > /etc/group
rm /passwd /shadow /group

在新服务器node6中创建用户家目录：

perl -e 'while (<>) { print "mkdir /home/$1; chown -R $1:$1 /home/$1\n" if (m/(\w+):\w+:(\d+)/ && $2 >= 1000 && $2 <= 2000); }' /etc/passwd | sh

在control服务器中，将各用户密钥信息和bash配置文件上传到新服务器node6中：

perl -e 'while (<>) { print "rm /home/$1/.ssh/known_hosts; scp -r /home/$1/.ssh /home/$1/.bash* node6:/home/$1/\n" if (m/(\w+):\w+:(\d+)/ && $2 >= 1000 && $2 <= 2000); }' /etc/passwd | sh

在新服务器node6中修改各用户目录权限： 

perl -e 'while (<>) { print "chown -R $1:$1 /home/$1; chmod 700 /home/$1\n" if (m/(\w+):\w+:(\d+)/ && $2 >= 1000 && $2 <= 2000); }' /etc/passwd | sh

4. 挂载共享文件夹

修改 /etc/fstab 文件，永久挂载其它服务器的共享文件夹。

192.168.30.4:/disks/control_8T    /disks/control_8T   nfs defaults    0   0
192.168.30.4:/disks/control_4T    /disks/control_4T   nfs defaults    0   0
192.168.30.4:/opt/    /opt    nfs defaults    0   0
192.168.30.4:/media/CentOS    /media/CentOS   nfs defaults    0   0

输入命令挂载所有的共享文件夹。

mount -a

5. 准备系统配置文件

修改/etc/profile.d/perl-homedir.sh配置文件，以免每次登录用户，自动在
家目录下生成perl5文件夹。

perl -p -i -e 's/PERL_HOMEDIR=1/PERL_HOMEDIR=0/' /etc/profile.d/perl-homedir.sh

修改 /etc/sudoers 配置文件，将 chenlianfu 用户变成超级管理员用户。

perl -p -i -e 's/^(root.*)/$1\nchenlianfu   ALL=(ALL)       NOPASSWD:ALL /' /etc/sudoers

修改 /etc/selinux/config 配置文件，永久关闭linux的一个安全机制，开启该安全机制会对很多操作造成阻碍。 

perl -p -i -e 's/SELINUX=enforcing/SELINUX=disabled/' /etc/selinux/config
setenforce 0

修改/etc/ssh/sshd_config配置文件，使openssh远程登录更安全，更快速。

perl -p -i -e 's/#RSAAuthentication/RSAAuthentication/' /etc/ssh/sshd_config
perl -p -i -e 's/#PubkeyAuthentication/PubkeyAuthentication/' /etc/ssh/sshd_config
perl -p -i -e 's/#AuthorizedKeysFile/AuthorizedKeysFile/' /etc/ssh/sshd_config
perl -p -i -e 's/.*PermitRootLogin.*/PermitRootLogin no/' /etc/ssh/sshd_config
perl -p -i -e 's/.*Protocol\s+2.*/Protocol 2/' /etc/ssh/sshd_config
perl -p -i -e 's/.*ClientAliveInterval.*/ClientAliveInterval 60/' /etc/ssh/sshd_config
perl -p -i -e 's/.*ClientAliveCountMax.*/ClientAliveCountMax 10/' /etc/ssh/sshd_config
perl -p -i -e 's/.*UseDNS.*/UseDNS no/' /etc/ssh/sshd_config
perl -p -i -e 's/GSSAPIAuthentication yes/GSSAPIAuthentication no/' /etc/ssh/sshd_config
systemctl restart sshd.service

修改用户对系统资源的最大权限：

echo "*	soft	nproc	10240
*	hard	nproc	102400
*	soft	nofile	10240
*	hard	nofile	102400
*	soft	stack	10240
*	hard	stack	102400" >> /etc/security/limits.conf

开放防火墙端口并使之永久生效

firewall-cmd --add-port=80/tcp --permanent
firewall-cmd --add-port=8080/tcp --permanent
firewall-cmd --add-port=3306/tcp --permanent

systemctl restart firewalld.service

systemctl start httpd.service
systemctl start mariadb.service

systemctl enable httpd.service
systemctl enable mariadb.service

设置并启动网页服务器

perl -i -e 'while (<>) { $mo = 1 if m#<Directory />#; $mo = 0 if m#<Files \".ht\*\">#; s/Require all denied/#Require all denied/ if $mo == 1; print; }' /etc/httpd/conf/httpd.conf
perl -p -i -e 's/(AddType text\/html .shtml)/$1\nAddType text\/plain .pl\nAddType text\/plain .py\nAddType text\/plain .sh/' /etc/httpd/conf/httpd.conf
perl -p -i -e 's/#AddHandler cgi-script .cgi/AddHandler cgi-script .cgi/' /etc/httpd/conf/httpd.conf

systemctl restart httpd.service

设置并启动Mysql服务器

/bin/rm -rf /etc/my.cnf /home/mysql/
cp /usr/share/mysql/my-large.cnf /etc/my.cnf 
perl -p -i -e 's/^datadir.*\n$//; s/(\[mysqld\])/$1\ndatadir         = \/home\/mysql/' /etc/my.cnf
mkdir /home/mysql
chown -R mysql:mysql /home/mysql
systemctl restart mariadb.service
/usr/bin/mysqladmin -u root password '123456'
echo "GRANT ALL ON *.* TO 'train'@'localhost' IDENTIFIED BY '123456'; GRANT SELECT ON *.* TO 'pasa'@'localhost' IDENTIFIED BY '123456'; GRANT ALL ON b2gdb.* TO 'blast2go'@'localhost' IDENTIFIED BY 'blast4it'; FLUSH PRIVILEGES;" | mysql -uroot -p123456

6. 安装系统软件

echo "alias yumlocal='yum --disablerepo=* --enablerepo=c7-media'" >> ~/.bashrc
source ~/.bashrc

yumlocal install -y gd* gcc* ftp cmake gsl* lftp screen lm_sensors gnuplot hdparm

wget https://tuxera.com/opensource/ntfs-3g_ntfsprogs-2017.3.23.tgz
tar zxf ntfs-3g_ntfsprogs-2017.3.23.tgz
cd ntfs-3g_ntfsprogs-2017.3.23
./configure && make -j 4 && make install
cd .. && /bin/rm ntfs-3g_ntfsprogs-2017.3.23 -rf
ln -s /usr/sbin/mount.ntfs-3g /usr/sbin/mount.ntfs

系统发育分析软件PAML的安装和简介

发表于由
3

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，则切换回仅显示进程的程序名。

按小写字母i，则将空闲的进程全部不显示。再按一次字母i，则切换显示所有的进程。

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

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

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

按小写字母z，能让切换到彩色显示；再按一次z切换到黑白显示。按大写字母Z可以设置不同的颜色：按S，输入数字0~7设置top头部的统计数据的颜色；按H， 输入数字0~7设置top表头部分的背景颜色； 按T， 输入数字0~7设置top表内进程数据的文字颜色； 按X， 输入数字0~7设置top表内处于运行状态进程数据的文字颜色； 按b，则对处于运行状态进程的背景进行高亮显示；按Enter键确认设置并退出设置界面。

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

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

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

设置好top命令的显示方式后，按W键保存其配置，下次再次运行top则按保存的方式运行。