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进行分歧时间计算

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系统的准备工作

在集群中新增加一台服务器,或对集群中的某台服务器重新安装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的安装和简介

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则按保存的方式运行。

解决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

使用Supernova对10x genomics测序数据进行基因组De novo组装

1. Supernovo软件下载和安装

在10x genomics官网的Supernovo下载页面填写资料,会给出下载命令。

curl -o supernova-2.1.1.tar.gz "http://cf.10xgenomics.com/releases/assembly/supernova-2.1.1.tar.gz?Expires=1558654333&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cDovL2NmLjEweGdlbm9taWNzLmNvbS9yZWxlYXNlcy9hc3NlbWJseS9zdXBlcm5vdmEtMi4xLjEudGFyLmd6IiwiQ29uZGl0aW9uIjp7IkRhdGVMZXNzVGhhbiI6eyJBV1M6RXBvY2hUaW1lIjoxNTU4NjU0MzMzfX19XX0_&Signature=Pnfu7Dty-bxyagrAeSN5Amqz~25hhExNHni4MIbtccZWrP590nE0MsmfxLWJsE9508Ae8rio4JmvlHw5LCk9blp~BGdLkkMqg5r1kgWOshEfKnsksvgNCzk6moOfdHEIRyUV1DwYFiy-M3vEDdRBFcNMVvGTbhGzRb~lmu5VVHV8JCRl6FCT1BvjtIe0xy4Jq6gpIgI3NRAREJZ3-Bbrko7JE0sJItsTwx6MpND-nmxvNkZ~iY69lx6zOWEqsRnxm4BlvLAkjuJ~6am0ROBaxAdDup8iaxYvLVGzZzW44GN7-f5OyRCCcbd6zo86QtRPzBNtZJbeqFg9R5W35VQI0g__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"

Supernovo软件当前最新版本为v2.1.1,解压缩安装后,直接是二进制预编译版本,可以直接在CentOS系统上使用。

tar zxf supernova-2.1.1.tar.gz -C /opt/biosoft/
echo 'PATH=$PATH:/opt/biosoft/supernova-2.1.1/' >> ~/.bashrc
source ~/.bashrc

2. 示例数据下载和Supernovo使用

Rfam数据库的Small nuclear RNA分类统计

常见的Non-coding RNA基因主要分四类:rRNA基因、tRNA基因、snRNA基因和miRNA基因。tRNA基因一般使用tRNAScan-SE软件进行分析;rRNA基因一般使用RNAmmer软件或BLAST软件进行分析;对于snRNA基因和miRNA基因,则使用Infernal软件将基因组序列比对到Rfam数据库进行分析得到。但Rfam数据库分析能得到很多Non-conding RNA信息,提取其中的snRNA信息则有点麻烦。具体分析方法如下:

1. snRNA简介

snRNA一类较小的存在细胞核中的RNA分子,通常~150个核苷酸长度。可以分成四类:

1. Sm-class snRNA: 由RNA聚合酶II转录,包含U1, U2, U4, U4atac, U5, U7, U11,和U12几类。该类snRNA能从核孔转移到细胞质中,最后形成splicesome,用于对mRNA前体进行加工,去除introns。

2. Lsm-class snRNA:和第一类snRNA类似,但由RNA聚合酶III转录,不会离开细胞核。

3. C/D box snoRNA: 存在核仁中的RNA分子,包含两个保守的序列motifs,C (RUGAUGA) and D (CUGA),通常能对rRNA进行甲基化加工。

4. H/ACA box snoRNA:存在核仁中的RNA分子,包含两个保守的序列motifs,H box (consensus ANANNA) and the ACA box (ACA),通常能对tRNA进行假尿苷化加工。

2. snoRNA分类和RF编号的获取

相比于snRNA,snoRNA家族庞大。若需要统计snoRNA信息,首先需要了解那些Rfam中的编号属于那种类型的snoRNA。

在Rfam数据库文件Rfam.cm中以关键词Small nucleolar RNA搜索,从而确定属于snoRNA的Rfam编号。

perl -e '$/ = "\n//"; while (<>) { if (m/Small nucleolar RNA/i) { my $name = $1 if m/NAME\s+(\S+)/; my $acc = $1 if m/ACC\s+(\S+)/; print "$acc\t$name\n"; } }' Rfam.cm > snoRNA_1.list

然后,根据Rfam编号在Rfam网站中鉴定所属的snoRNA分类,C/D box或 H/ACA box。此步骤需要编写如下程序snoRNA_type.pl,来自动下载相应Rfam编号的网页信息,并解析其所属分类。

#!/usr/bin/perl
use strict;

my $info = `curl http://rfam.xfam.org/family/$ARGV[0]`;

if ($info =~ m/Gene; snRNA; snoRNA; ([^\s;]+)/) {
    print "$ARGV[0]\t$1\n";
}
else {
    print "$ARGV[0]\tother\n";
}

批量化运行解析程序,并合并结果得到snoRNA的Rfam编号和分类信息。

for i in `cut -f 1 snoRNA_1.list`
do
    echo "perl snoRNA_type.pl $i > $i.type"
done > command.snoRNA_type.list
ParaFly -c command.snoRNA_type.list -CPU 50

cat *.type > snoRNA_2.list
rm *.type
cut -f 2 snoRNA_1.list > aa
paste snoRNA_2.list aa > snoRNA.list
rm aa snoRNA_1.list snoRNA_2.list snoRNA_type.pl

3. 编写程序统计snRNA信息

程序名称:snRNA_stats_from_Rfam.pl

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 [options] rfam_out.tab snoRNA_type.txt

    --out-prefix <string>    default: snRNA_out
    设置输出文件前缀。默认输出文件:
    (1)snRNA_out.stats 对各类snRNA的统计信息。
    (2)snRNA_out.txt 对Rfam数据库进行Infernal分析中属于snRNA的表格结果。
    (3)snRNA_out.gff3 将第2个结果转换为GFF3格式的结果。

    程序通过读取Rfam的infernal分析结果,统计出small nuclear RNA的信息。snRNA_out.stat结果文件包含多列:
    第1列:snRNA分类
    第2列:snRNA名称
    第2列:RF编号
    第3列:基因家族数量
    第4列:序列平均长度
    第5列:序列总长度

    small nuclear RNA分四大类:
    Sm-class snRNA: (U1, U2, U4, U4atac, U5, U6, U7, U11, U12)
    Lsm-class snRNA: (U6, U6atac)
    C/D box snoRNA
    H/ACA box snoRNA

    snoRNA是最复杂的一类,在Rfam v14.1中包含565个家族,其中388个属于C/D box,177个属于H/ACA box。通过读取snoRNA_type.txt文件中RF编号和snoRNA分类对应的信息,来统计两类snoRNA的数量。

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

my ($out_prefix);
GetOptions(
    "out-prefix:s" => \$out_prefix,
);
$out_prefix ||= "snRNA_out";

# Sm-class或Lsm-class类型的snRNA在Rfam数据库中名称和RF编号获取:
# perl -e '$/ = "\n//"; while (<>) { if (m/small nuclear RNA/i or m/spliceosomal RNA/i) { my $name = $1 if m/NAME\s+(\S+)/; my $acc = $1 if m/ACC\s+(\S+)/; print "$acc\t$name\n"; } }' Rfam.cm | uniq

my %snRNA_RF2Name = (
    "RF00003" => "U1",
    "RF00004" => "U2",
    "RF00007" => "U12",
    "RF00015" => "U4",
    "RF00020" => "U5",
    "RF00026" => "U6",
    "RF00066" => "U7",
    "RF00488" => "U1_yeast",
    "RF00548" => "U11",
    "RF00618" => "U4atac",
    "RF00619" => "U6atac",
    "RF01844" => "SmY",
    "RF02491" => "Gl_U1",
    "RF02492" => "Gl_U2",
    "RF02493" => "Gl_U4",
    "RF02494" => "Gl_U6"
);

my %snRNA_RF2Class = (
    "RF00003" => "Sm-class snRNA",
    "RF00004" => "Sm-class snRNA",
    "RF00007" => "Sm-class snRNA",
    "RF00015" => "Sm-class snRNA",
    "RF00020" => "Sm-class snRNA",
    "RF00026" => "Lsm-class snRNA",
    "RF00066" => "Sm-class snRNA",
    "RF00488" => "Sm-class snRNA",
    "RF00548" => "Sm-class snRNA",
    "RF00618" => "Sm-class snRNA",
    "RF00619" => "Lsm-class snRNA",
    "RF01844" => "Sm-class snRNA",
    "RF02491" => "Sm-class snRNA",
    "RF02492" => "Sm-class snRNA",
    "RF02493" => "Sm-class snRNA",
    "RF02494" => "Sm-class snRNA"
);

open IN, $ARGV[1] or die "Can not open file $ARGV[1], $!\n";
while (<IN>) {
    next if m/^#/;
    next if m/^\s*$/;
    chomp;
    @_ = split /\t/;
    $snRNA_RF2Name{$_[0]} = $_[2];

    if ($_[1] eq "CD-box") {
        $snRNA_RF2Class{$_[0]} = 'C/D box snoRNA';
    }
    elsif ($_[1] eq "HACA-box") {
        $snRNA_RF2Class{$_[0]} = 'H/ACA box snoRNA';
    }
    else {
        $snRNA_RF2Class{$_[0]} = 'other';
    }
}
close IN;

my %class;
foreach (keys %snRNA_RF2Class) {
    $class{$snRNA_RF2Class{$_}}{$_} = 1;
}

open OUT1, ">", "$out_prefix.txt" or die "Can not create file $out_prefix.txt, $!\n";
open OUT2, ">", "$out_prefix.gff3" or die "Can not create file $out_prefix.gff3, $!\n";
open OUT3, ">", "$out_prefix.stats" or die "Can not create file $out_prefix.stats, $!\n";
my (%stats, %gff3, $total_snRNA_number);
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
while (<IN>) {
    if (m/^#/) { print OUT1; next; }
    next if m/^\s*$/;

    @_ = split /\s+/;
    if (exists $snRNA_RF2Name{$_[3]}) {
        $total_snRNA_number ++;
        print OUT1;
        $gff3{$_} = 1;

        $stats{$_[3]}{"num"} ++;
        $stats{$_[3]}{"total_length"} += (abs($_[8] - $_[7]) + 1);
    }
}
close IN;

my @class = ("Sm-class snRNA", "Lsm-class snRNA", 'C/D box snoRNA', 'H/ACA box snoRNA');
push @class, "other" if exists $class{"other"};
foreach my $class (@class) {
    my @rf = keys %{$class{$class}};
    my (%out1 , %out2);
    foreach (keys %{$class{$class}}) {
        my $num = $stats{$_}{"num"};
        next unless $num;
        my $total_length = $stats{$_}{"total_length"};
        my $avg = int(($total_length / $num) * 100 + 0.5) / 100;
        my $out = "$class\t$snRNA_RF2Name{$_}\t$_\t$num\t$avg\t$total_length";
        $out1{$out} = $num;
        $out2{$out} = $total_length;
    }
    foreach (sort {$out1{$b} <=> $out1{$a} or $out2{$b} <=> $out2{$a}} keys %out1) {
        print OUT3 "$_\n";
    }
}

my (%sortGFF1, %sortGFF2, %sortGFF3, %sortGFF4, %gff3_name, %gff3_RF, $number);
foreach (keys %gff3) {
    @_ = split /\s+/, $_;
    my ($seqID, $start, $end, $score, $strand, $name, $rf_code) = ($_[0], $_[7], $_[8], $_[14], $_[9], $_[2], $_[3]);
    if ($start > $end) {
        my $tmp_vaule = $start;
        $start = $end;
        $end = $tmp_vaule;
    }

    my $out = "$seqID\t.\tsnRNA\t$start\t$end\t$score\t$strand\t\.";
    $sortGFF1{$out} = $seqID;
    $sortGFF2{$out} = $start;
    $sortGFF3{$out} = $end;
    $sortGFF4{$out} = $strand;
    $gff3_name{$out} = $name;
    $gff3_RF{$out} = $rf_code;
}

foreach (sort {$sortGFF1{$a} cmp $sortGFF1{$b} or $sortGFF2{$a} <=> $sortGFF2{$b} or $sortGFF3{$a} <=> $sortGFF3{$b} or $sortGFF4{$a} cmp $sortGFF4{$b}} keys %sortGFF1) {
    $number ++;
    my $id = '0' x (length($total_snRNA_number) - length($number)) . $number;
    print OUT2 "$_\tID=snRNA_$id;Name=$gff3_name{$_};RF=$gff3_RF{$_}\n";
}
close OUT1; close OUT2; close OUT3;

4. 编写程序统计miRNA信息

程序名称:miRNA_stats_from_Rfam.pl

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 [options] rfam_out.tab Rfam.cm

    --out-prefix <string>    default: miRNA_out
    设置输出文件前缀。默认输出文件:
    (1)miRNA_out.stats 对各类microRNA的统计信息。
    (2)miRNA_out.txt 对Rfam数据库进行Infernal分析中属于microRNA的表格结果.
    (3)miRNA_out.gff3 将第2个结果转换为GFF3格式的结果。

    程序通过读取Rfam数据库信息,使用关键词microRNA搜索,找到属于miRNA的编号,再分析infernal结果文件,统计miRNA信息。

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

my ($out_prefix);
GetOptions(
    "out-prefix:s" => \$out_prefix,
);
$out_prefix ||= "miRNA_out";

open IN, $ARGV[1] or die "Can not open file $ARGV[1], $!\n";
my %miRNA_RF2Name;
$/ = "\n//";
while (<IN>) {
    if (m/microRNA/i) {
        my $name = $1 if m/NAME\s+(\S+)/;
        my $acc = $1 if m/ACC\s+(\S+)/;
        $miRNA_RF2Name{$acc} = $name;
    }
}
close IN;
$/ = "\n";

open OUT1, ">", "$out_prefix.txt" or die "Can not create file $out_prefix.txt, $!\n";
open OUT2, ">", "$out_prefix.gff3" or die "Can not create file $out_prefix.gff3, $!\n";
open OUT3, ">", "$out_prefix.stats" or die "Can not create file $out_prefix.stats, $!\n";
my (%stats, %gff3, $total_miRNA_number);
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
while (<IN>) {
    if (m/^#/) { print OUT1; next; }
    next if m/^\s*$/;

    @_ = split /\s+/;
    if (exists $miRNA_RF2Name{$_[3]}) {
        $total_miRNA_number ++;
        print OUT1;
        $gff3{$_} = 1;

        $stats{$_[3]}{"num"} ++;
        $stats{$_[3]}{"total_length"} += (abs($_[8] - $_[7]) + 1);
    }
}
close IN;

foreach (sort {$stats{$b}{"num"} <=> $stats{$a}{"num"}} keys %stats) {
    my $num = $stats{$_}{"num"};
    next unless $num;
    my $total_length = $stats{$_}{"total_length"};
    my $avg = int(($total_length / $num) * 100 + 0.5) / 100;
    print OUT3 "$miRNA_RF2Name{$_}\t$_\t$num\t$avg\t$total_length\n";
}

my (%sortGFF1, %sortGFF2, %sortGFF3, %sortGFF4, %gff3_name, %gff3_RF, $number);
foreach (keys %gff3) {
    @_ = split /\s+/, $_;
    my ($seqID, $start, $end, $score, $strand, $name, $rf_code) = ($_[0], $_[7], $_[8], $_[14], $_[9], $_[2], $_[3]);
    if ($start > $end) {
        my $tmp_vaule = $start;
        $start = $end;
        $end = $tmp_vaule;
    }

    my $out = "$seqID\t.\tmiRNA\t$start\t$end\t$score\t$strand\t\.";
    $sortGFF1{$out} = $seqID;
    $sortGFF2{$out} = $start;
    $sortGFF3{$out} = $end;
    $sortGFF4{$out} = $strand;
    $gff3_name{$out} = $name;
    $gff3_RF{$out} = $rf_code;
}

foreach (sort {$sortGFF1{$a} cmp $sortGFF1{$b} or $sortGFF2{$a} <=> $sortGFF2{$b} or $sortGFF3{$a} <=> $sortGFF3{$b} or $sortGFF4{$a} cmp $sortGFF4{$b}} keys %sortGFF1) {
    $number ++;
    my $id = '0' x (length($total_miRNA_number) - length($number)) . $number;
    print OUT2 "$_\tID=miRNA_$id;Name=$gff3_name{$_};RF=$gff3_RF{$_}\n";
}
close OUT1; close OUT2; close OUT3;

服务器集群设置

使用服务器集群能有利于快速地进行大数据分析,但集群服务器管理相对复杂。根据本人经验讲解服务器集群初始安装和配置。

1. 服务器集群硬件信息和CentOS7系统安装

现有服务器集群共12台服务器:包含1台管理节点node1、胖计算节点1台node10,普通计算节点10台,从node11到node20。此外:有一台直连式存储连接到管理节点node1、所有12台服务器连接到infiniband交换机上。

对这12台服务器统一安装CentOS7_1810系统,并最大化安装系统。安装系统的时候可以设置各服务器的主机名。此外,安装系统时,可以在node1服务器上,设置将存储挂载到 /disk 目录下,也可以在安装系统完毕后,通过修改 /etc/fstab 文件实现:

cat << EOF >> /etc/fstab
/dev/mapper/mpatha      /disk                   ext4    defaults        1 2
EOF

2. 服务器CentOS系统安装完毕后的基本配置

使用root用户在管理节点和计算节点上进行以下操作:

  1. 修改/etc/profile.d/perl-homedir.sh配置文件,以免每次登录用户,自动在家目录下生成perl5文件夹
  2. 修改/etc/sudoers配置文件,将自己的用户(例如 train)变成超级管理员用户
  3. 修改/etc/selinux/config配置文件,永久关闭linux的一个安全机制,开启该安全机制会对很多操作造成阻碍。
  4. 修改/etc/ssh/sshd_config配置文件,使openssh远程登录更安全,更快速
  5. 增加系统资源权限
perl -p -i -e 's/PERL_HOMEDIR=1/PERL_HOMEDIR=0/' /etc/profile.d/perl-homedir.sh
echo 'eval "$(perl -Mlocal::lib=$HOME/.perl5)"' >> ~/.bashrc

perl -i.bak -e 'while (<>) { if (/^root/) { print; print "train   ALL=(ALL)       NOPASSWD:ALL\n"; last; } else { print } }' /etc/sudoers

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

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

perl -p -i -e 's/^\*.*\n$//' /etc/security/limits.conf
cat << EOF >> /etc/security/limits.conf
*       soft    nofile  10240
*       hard    nofile  102400
*       soft    stack   10240
*       hard    stack   102400
*       soft    core    unlimited
*       hard    core    unlimited
*       soft    nproc   10240
*       hard    nproc   102400
EOF

3. 配置集群中各服务器的主机名和IP地址

使用root用户在管理节点 和计算节点服务器上对infiniband网口进行配置,修改 /etc/sysconfig/network-scripts/ifcfg-ib0 配置文件内容:

BOOTPROTO=none
ONBOOT=yes
IPADDR=12.12.12.xx # xx要修改成对应节点的编号。例 node11节点这要修改成11
PREFIX=24
GATEWAY=12.12.12.1

修改好ifcfg文件后,重启网络服务,使生效:

systemctl restart network.service

各节点服务器在infiniband网络之间的联通需要在控制节点node1上安装一些相关的系统软件,并启用相应服务:

yum install opensm* opensm-devel* infiniband-diags perftest* gperf* opensm*
systemctl restart opensm.service
systemctl enable rdma.service
systemctl enable opensm.service

然后将所有节点服务器的 /etc/hosts 文件内容修改成同样的内容:

cat << EOF > /etc/hosts
12.12.12.1      node1
12.12.12.10     node10
12.12.12.11     node11
12.12.12.12     node12
12.12.12.13     node13
12.12.12.14     node14
12.12.12.15     node15
12.12.12.16     node16
12.12.12.17     node17
12.12.12.18     node18
12.12.12.19     node19
12.12.12.20     node20
EOF

4. 将控制节点的以太网共享给计算节点

控制节点通过电信100M宽带连接外网,通过网线将node1控制节点连接到电信网关(光猫和路由器合一的电信盒子)上。设置网口自动使用HDCP方法分配IP地址即可。在外网可以正常连接的情况,可以将该网络通过infiniband网卡共享给其它计算节点。

在node1控制节点上使用root用户进行操作:

  • 开启NAT转发
  • 开放DNS使用的53端口并重启防火墙,否则可能导致内网服务器虽然设置正确的DNS,但是依然无法进行域名解析。
  • 控制节点上是在ens5网口连接外网,对其网络进行共享。
firewall-cmd --permanent --zone=public --add-masquerade

firewall-cmd --zone=public --add-port=53/tcp --permanent
systemctl restart firewalld.service

echo 'net.ipv4.ip_forward=1' >> /etc/sysctl.conf
sysctl -p
firewall-cmd --permanent --direct --passthrough ipv4 -t nat -I POSTROUTING -o ens6 -j MASQUERADE -s 12.12.12.0/24
systemctl restart firewalld.service

在计算节点上对infiniband网卡进行IP设置时,将网关设置成提供网络的主机IP即可,即将网关设置成node1管理节点的IP地址12.12.12.1。这在上一步已经设置过了。

5. 将控制节点的大容量存储共享给计算节点

在控制节点node1服务器上,修改NFS配置文件/etc/sysconfig/nfs配置文件,打开所有带有PORT的注释行,表示使用相应的防火墙端口,并修改防火墙配置,开放对应端口:

perl -p -i -e 's/^#(.*PORT)/$1/' /etc/sysconfig/nfs

firewall-cmd --add-port=32803/udp --permanent
firewall-cmd --add-port=32803/tcp --permanent
firewall-cmd --add-port=32769/udp --permanent
firewall-cmd --add-port=32769/tcp --permanent
firewall-cmd --add-port=892/udp --permanent
firewall-cmd --add-port=892/tcp --permanent
firewall-cmd --add-port=662/udp --permanent
firewall-cmd --add-port=662/tcp --permanent
firewall-cmd --add-port=2020/udp --permanent
firewall-cmd --add-port=2020/tcp --permanent
firewall-cmd --add-port=875/udp --permanent
firewall-cmd --add-port=875/tcp --permanent

systemctl restart firewalld.service

然后,在控制节点node1服务器上,启动NFS服务,并设置成开机启动:

systemctl restart rpcbind.service
systemctl restart nfs.service

systemctl enable rpcbind.service
systemctl enable nfs.service

继续,在控制节点node1服务器上, 修改/etc/exports文件内容,添加被共享的文件夹信息,并使配置生效:

cat << EOF >> /etc/exports
/disk   12.12.12.0/24(rw,sync,no_root_squash,no_subtree_check)
/opt    12.12.12.0/24(rw,sync,no_root_squash,no_subtree_check)
EOF

exportfs -rv

在各计算节点服务器上,使用root用户修改配置文件/etc/fstab,对node1的共享文件夹进行挂载:

mkdir /disk

cat << EOF >> /etc/fstab
12.12.12.1:/disk        /disk   nfs     defaults        0       0
12.12.12.1:/opt /opt    nfs     defaults        0       0
EOF

mount -a

6. 在集群计算机上创建新用户

首先,生成文件/disk/users.txt。该文件每行一个待生成的用户名。

然后,在所有节点服务器中进行操作,生成用户并使用create_random_passwd.pl命令赋予随机密码:

cd /disk
for i in `cat users.txt`
do
    useradd $i 2> /dev/null
    ./create_random_passwd.pl $i
done

在控制节点node1服务器中进行操作:在大容量存储对应的共享文件夹中建立新用户的专属文件夹;使用root用户生成新用户的ssh密钥对数据和授权文件信息并放入到各新用户的家目录下。

/bin/rm /disk/ssh_info/ -rf
mkdir -p /disk/ssh_info/
for i in `cat users.txt`
do
    mkdir /disk/ssh_info//$i /disk/$i
    chown -R $i:$i /disk/$i
    chmod 700 /disk/$i
    ssh-keygen -t dsa -P '' -f /disk/ssh_info/$i/id_dsa
    chown -R $i:$i /disk/ssh_info/$i
    mkdir /home/$i/.ssh
    /bin/cp -a /disk/ssh_info/$i/* /home/$i/.ssh
    chown -R $i:$i /home/$i/.ssh
    chmod 700 /home/$i/.ssh
    cat /disk/ssh_info/$i/id_dsa.pub >> /home/$i/.ssh/authorized_keys
    chown -R $i:$i /home/$i/.ssh/authorized_keys
    chmod 600 /home/$i/.ssh/authorized_keys
done

在各个计算节点服务器中使用root用户将上一步生成的ssh密钥对数据和授权文件信息放入到计算节点服务器中各新用户的家目录下:

cd /disk

for i in `cat users.txt`
do
    useradd $i 2> /dev/null
    ./create_random_passwd.pl $i
done

for i in `cat users.txt`
do
    mkdir /home/$i/.ssh 
    /bin/cp -a /disk/ssh_info/$i/* /home/$i/.ssh
    chown -R $i:$i /home/$i/.ssh
    chmod 700 /home/$i/.ssh
    cat /disk/ssh_info/$i/id_dsa.pub >> /home/$i/.ssh/authorized_keys
    chown -R $i:$i /home/$i/.ssh/authorized_keys
    chmod 600 /home/$i/.ssh/authorized_keys
done

create_random_passwd.pl程序代码:

#!/usr/bin/perl
#use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 [options] username

    使用root用户执行该程序,输入用户名,则能调用passwd命令给该用户创建一个随机密码。并将用户名及其密码输出到标准输出。
    --length <int>    default:10
    设置生成密码的字符长度。

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

my $length;
GetOptions(
    "length:i" => \$length,
);
$length ||= 10;
my @cha = ('!', '@', '#', '$', '%', '^', '&', '*', '.', '_');
foreach (0..9) {
    push @cha, $_;
}
foreach (a..z) {
    push @cha, $_;
}
foreach (A..Z) {
    push @cha, $_;
}

my $passwd;
for (1..$length) {
    my $cha_num = rand(@cha);
    $passwd .= $cha[$cha_num];
}
print "$ARGV[0]\t$passwd\n";

my $cmdString = "echo \'$passwd\' | passwd --stdin $ARGV[0] &> /dev/null";
(system $cmdString) == 0 or die "Faield to excute: $cmdString, $!\n";

7. 远程桌面软件vncserver安装和使用

由于控制节点node1是连接到了电信网关上,没有固定IP地址,推荐使用vnc来对内网服务器使用图形化桌面方法进行控制。

首先,使用root用户在node1服务器上进行操作,安装vncserver软件并开放相应的防火墙5901,5902,5903端口:

yum install vcn vnc-server

firewall-cmd --zone=pulic --add-port=5901/tcp --permanent
firewall-cmd --zone=pulic --add-port=5902/tcp --permanent
firewall-cmd --zone=pulic --add-port=5903/tcp --permanent
systemctl restart firewalld.service

然后,使用普通用户(例如,train)开启vncserver服务:

vncserver
# 第一次启动需要输入密码

进行其它vnc操作并修改桌面分辨率,提供更好的vnc体验:

查看当前开启的vncserver桌面列表
vncserver -list

查看第一个vncserver桌面的端口号
cat ~/.vnc/node1\:1.log

关闭第一个vncserver桌面
vncserver -kill :1

修改vncserver桌面的分辨率
cat << EOF >> .vnc/config
geometry=2000x1052
EOF

关闭后再次启动vncserver桌面,则分辨率变得更好了
vncserve

为了让vnc能在外网对node1进行控制。需要将node1控制节点服务器和公网服务器使用ssh进行连接,开启反向隧道,并进行端口转发,在node1服务器上进行操作。以下命令将node1服务器VNC服务对应的5901端口映射到公网服务器115.29.105.12的4497端口上:

ssh -N -f -R 4497:localhost:5901 train@115.29.105.12

注意,以上命令需要在公网服务器115.29.105.12上拥有train用户的密码,才能ssh连接成功;并且,还需要使用该公网服务器的root用户开启4497防火墙端口,同时在ssh配置文件设置允许端口转发,才能使vnc访问生效。

最后,在windows系统下下载nvcviewer软件,然后安装并打开软件,输入115.29.105.12:4497,再输入之前设置的密码,即可访问远程桌面。

8. 在控制节点上控制计算节点的开机和关机

在控制节点上,对计算节点可以使用ssh连接并导入shutdown指令的方法进关机。基于此原理,编写名为 guanji 的Perl程序来对指定的节点进行关机。该程序代码:

#!/usr/bin/perl
use strict;

my $usage = <<USAGE;
Usage:
    $0 node10 node11 node12 ...

    使用此命令关闭目标节点。该命令后可以输入1个或多个主机名,关闭相应的计算节点。
    若命令后输入的主机名中有一个是all,则会关闭所有的计算节点(从node11到node20)。
    此外,支持node11-node15这样中间带有中划线的输入方法,表示多个连续的节点。

For example:
    $0 node11 node13-node16 node20

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

my @node = qw/node10 node11 node12 node13 node14 node15 node16 node17 node18 node19 node20/;
my %node;
foreach (@node) { $node{$_} = 1; }

my %target;
foreach (@ARGV) {
    if ($_ eq "all") {
        foreach (@node) { $target{$_} = 1; }
        last;
    }
    elsif (m/(\d+)-node(\d+)/) {
        foreach ($1 .. $2) {
            $target{"node$_"} = 1;
        }
    }
    else {
        if (exists $node{$_}) {
            $target{$_} = 1;
        }
        else {
            print STDERR "Warning: $_不是能控制的目标节点。\n";
        }
    }
}

foreach (sort keys %target) {
    &guanji($_);
}

sub guanji {
    print STDERR "正在检测到 $_ 的连接\n";
    my $ping = `ping $_ -c 1`;
    if ($ping =~ m/Unreachable/) {
        print STDERR "Warning: $_连接失败,可能已经处于关机状态。\n";
    }
    else {
        my $cmdString = "ssh $_ 'sudo shutdown -h now' &> /dev/null";
        system $cmdString;
        print "对主机 $_ 已经发送关机指令\n";
    }
}

在控制节点node1上,可以使用wol软件基于网络唤醒的方法对计算节点进行开机。基于此原理,编写名为 kaiji 的Perl程序对指定节点进行开机。该程序代码:

#!/usr/bin/perl
use strict;

my $usage = <<USAGE;
Usage:
    $0 node10 node11 node12 ...

    使用此命令开启目标节点。该命令后可以输入1个或多个主机名,开启相应的计算节点。
    若命令后输入的主机名中有一个是all,则会开启所有的计算节点(从node11到node20)。
    此外,支持node11-node15这样中间带有中划线的输入方法,表示多个连续的节点。

For example:
    $0 node11 node13-node16 node20

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

my @node = qw/node10 node11 node12 node13 node14 node15 node16 node17 node18 node19 node20/;
my %node = ("node10" => "00:e0:ec:27:e9:a0",
"node11" => "e8:61:af:11:e9:4b",
"node12" => "e8:61:af:11:e8:3f",
"node13" => "e8:61:af:1b:ec:80",
"node14" => "e8:61:af:1b:ed:84",
"node15" => "e8:61:af:1b:ec:9e",
"node16" => "e8:61:af:1b:ed:0e",
"node17" => "e8:61:af:1b:ed:b4",
"node18" => "e8:61:af:1b:ec:94",
"node19" => "e8:61:af:1b:ec:5a",
"node20" => "e8:61:af:1b:eb:d0");

my %target;
foreach (@ARGV) {
    if ($_ eq "all") {
        foreach (@node) { $target{$_} = 1; }
        last;
    }
    elsif (m/(\d+)-node(\d+)/) {
        foreach ($1 .. $2) {
            $target{"node$_"} = 1;
        }
    }
    else {
        if (exists $node{$_}) {
            $target{$_} = 1;
        }
        else {
            print STDERR "Warning: $_不是能控制的目标节点。\n";
        }
    }
}

foreach (sort keys %target) {
        &kaiji($_);
}

sub kaiji {
        print "对主机 $_ 已经发送开机指令\n";
        my $cmdString = "/opt/sysoft/wol-0.7.1/bin/wol --host=10.10.10.255 $node{$_}";
        system $cmdString;
}

9. 在集群上部署SGE网格系统

开启防火墙端口并重启防火墙:

firewall-cmd --add-port=992/udp --permanent
firewall-cmd --add-port=6444/tcp --permanent
firewall-cmd --add-port=6445/tcp --permanent
systemctl restart firewalld.service

在控制节点node1服务器上安装SGE软件:

yum install csh java-1.8.0-openjdk java-1.8.0-openjdk-devel gcc ant automake hwloc-devel openssl-devel libdb-devel pam-devel libXt-devel motif-devel ncurses-libs ncurses-devel
yum install ant-junit junit javacc

wget https://arc.liv.ac.uk/downloads/SGE/releases/8.1.9/sge-8.1.9.tar.gz -P ~/software/
tar zxf ~/software/sge-8.1.9.tar.gz
cd sge-8.1.9/source
./scripts/bootstrap.sh
./aimk -no-herd -no-java

mkdir /opt/sysoft/sge
export SGE_ROOT=/opt/sysoft/sge
./scripts/distinst -local -allall -noexit # 虽然普通用户在目标文件夹有写权限,但是程序要对一些文件进行权限修改,使用root用户不会报错。
cd ../../ && rm sge-8.1.9/ -rf

在控制节点node1服务器上部署SGE:

cd $SGE_ROOT
./install_qmaster
# 注意在某个交互式界面输入/etc/hosts文件,以便让集群中所有服务器能倍识别。

在控制节点node1启动SGE,载入SGE软件环境变量:

echo 'source /opt/sysoft/sge/default/common/settings.sh' >> ~/.bashrc
source ~/.bashrc

/opt/sysoft/sge/default/common/sgemaster

在其它计算节点服务器上部署SGE:

cd $SGE_ROOT
./install_execd

其它计算节点上可能不会自动开机启动sge服务。手动设置计算节点上的sge开机启动:

cat <<EOF > sge.service
[Unit]
Description=SGE sgeexecd
After=network.target remote-fs.target nss-lookup.target

[Service]
Type=forking
User=chenlianfu
Group=chenlianfu
ExecStart=/opt/sysoft/sge/default/common/sgeexecd start
ExecReload=/opt/sysoft/sge/default/common/sgeexecd restart
ExecStop=/opt/sysoft/sge/default/common/sgeexecd stop
PrivateTmp=true

[Install]
WantedBy=multi-user.target
EOF

chmod 755 sge.service

# 在每台计算节点上执行如下命令
sudo cp sge.service /usr/lib/systemd/system
sudo systemctl restart sge.service
sudo systemctl enable sge.service