使用SAMR对蛋白组数据表达量进行差异分析

1. SAMR简介

SAM(Significance Analysis of Microarrays)在基因芯片数据时代中被开发出来进行基因表达量差异分析。该算法也能用于进行RNA-Seq数据的基因表达量差异分析,但貌似较少人会用它进行RNA-Seq数据分析。

最近在一篇对蛋白组数据差异分析软件进行比较的文章中,SAM结果表现最优。本文对基于R软件的SAM算法软件SAMR的使用进行简单讲述。

2. SAMR软件的安装和启动

打开R软件,输入如下命令安装SAMR相关的包:

$ R
> install.packages(c("samr", "matrixStats", "GSA", "shiny", "openxlsx"))
> source("http://bioconductor.org/biocLite.R")
> biocLite("impute")

要注意的是,对openxlsx包的安装可能会失败(我使用的是R-3.2.0版本),则在上述命令中去除对openxlsx的安装,选择手动下载并安装老版本的openxlsx包:

$ wget https://cran.r-project.org/src/contrib/Archive/openxlsx/openxlsx_2.4.0.tar.gz
$ R CMD INSTALL openxlsx_2.4.0.tar.gz

启动SAMR软件:

$ R
> library(shiny)
> runGitHub("SAM", "MikeJSeo")

启动SAMR软件,则会自动打开CentOS 6.8系统自带的火狐浏览器,进入软件的网页界面。

3. SARM软件使用

3.1 输入文件准备

软件的输入文件必须是xlsx格式的EXCEL文件。进行蛋白组表达量数据进行分析,其文件内容要求如下:

1. 第一列是基因Name,第二列是基因ID,每个基因ID应该独一无二。
2. 第一行表示样品名,第一行第一列和第一行第二列是空着的,从第三列开始表示样品名,且样品名仅能是数字1或2,代表两个不同的样品。都使用数字1表示样品1的多个重复,都使用数字2表示样品2的多个重复。虽然SARM能分析多个样品的数据,但是其结果格式不好,一般是进行两两比较。
3. 若有数据缺失,则某行某列不填入数据,或填入非数字(推荐 NA)代替。
4. 若有多组样品进行比较,则需要准备多个excel文件,而不是一个excel文件,多个sheet。

3.2 在网页中进行SAMR操作

按如下步骤进行操作:

1. 点击 Browse,选中目标excel文件。
2. 在 Minimum fold change 一栏中填写一个最小的差异倍数值,比如 1.5 或 2 等。
3. 在 Data type 一栏选择默认的 Array 。
4. 在 Response Type 一栏选择 Two class unpaired。这个选择项有很多,与输入文件的数据格式相关,我们是对两个样品进行比较,故选择此选项。输入文件第一列不是1和2,是不同的数值(比如血压),此处则选择Quantitative。
5. 在 Analysis Type 一栏选择默认的 Standard 即可。
6. 在 Test Statistic 一栏选择默认的 T-statistic 即可。
7. 在 Median Center the arrarys 一栏选择 Yes。选择Yes表示软件会根据中位数来对数据进行标准化,即让每一列的中位数都变成零。这样每个样品的数据都具有相同的中位数,相当于进行了样品间的标准化。软件推荐在输入数据之前进行标准化(不如TMM算法进行标准化),若输入的数据是标准化后的数据,则此栏选择默认的 No。
8. 其它参数选择一般情况下选择默认的即可。
9. 点击左上角的 Run 按钮进行分析,分开得到结果,在右边各个标签栏进行查看结果。
10. 在 Delta Table 页面下查看 Median FDR 值,从上往下找到该值变为0.05时的delta值,再在左侧 Delta 一栏选择 manually Enter Delta, 然后再在 Delta value 一栏填入该值,再进行计算。然后点击Current Settings标签页,看其False Discovery Rate的值,手动调整 Delta 值,直到 FDR 值最接近并低于 0.05 为止。Delta值越大,FDR值越小。
11. 在 Paste the filepath to save the output 一栏填入需要输出的文件夹路径,该路径一定要存在,在其下一栏填入输出文件的前缀,比如 out 。 在点击 Save,然后会输出结果文件 out.xlsx,该excel文件有多个sheet分别多个标签页中的结果。

4. SARM 软件算法和原理

这个比较复杂,我也不怎么搞懂,包括结果中一些Delta, Score(d), q-value, FDR, localfdr等,有点糊涂。等以后有时间搞更明白了,再添加解释。

Delta value:
Score(d):T检验的的d值,由 Numerator(r) / Denominator(s+s0) 得到。

FALCON

1. 安装FALCON

$ cd /opt/biosoft
$ export GIT_SYM_CACHE_DIR=~/.git-sym-cache # to speed things up
$ git clone git://github.com/PacificBiosciences/FALCON-integrate.git
$ cd FALCON-integrate
$ git checkout master  # or whatever version you want
$ git submodule update --init # Note: You must do this yourself! No longer via `make init`.
$ make init
$ source env.sh
$ make config-edit-user
$ make -j all
$ make test  # to run a simple one

2. FALCON运行示例

先设置好环境变量
$ cd /opt/biosoft/FALCON-integrate/
$ source env.sh

下载基因组大小只有200kb的greg200k数据,以及基因组大小有4.6Mb的E. coli数据
$ cd FALCON-examples/.git-sym/
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/greg200k-sv2.tar | tar xvf -
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar | tar xvf -

对200kb基因组进行组装
$ cd /opt/biosoft/FALCON-integrate/FALCON-examples/run/greg200k-sv2
$ fc_run.py fc_run.cfg

对4.6Mb基因组进行组装
$ cd /opt/biosoft/FALCON-integrate/FALCON-examples/run/ecoli
$ perl -p -i -e 's/^#job_type/job_type/' fc_run.cfg
$ fc_run.py fc_run.cfg

生成的最终结果文件为 2-asm-falcon/p_ctg.fa 。

3. FALCON工作原理

1. 对原始subreads进行Overlapping分析
2. 进行预组装,对Subreads进行校正
3. 对校正后的Subreads进行Overlapping分析
4. 过滤不可靠的Overlaps结果
5. 根据可靠的Overlap构建Grahp
6. 由Grahp得到contigs序列

4. FALCON的配置文件

FALCON软件的输入文件是Pacbio测序得到的Fasta文件。将这些Fasta文件的路径写入到文件 input.fofn 中用于软件输入。

FALCON软件运行的难点在于其配置文件fc_run.cfg的编写。一个fc_run.cfg示例如下:

[General]
## 设置任务提交方式
# jobtype的默认值应该是SGE,表示使用SGE集群进行计算;
# local表示使用本地主机运行FALCON,兼容性最好,毕竟集群的部署非常麻烦;
job_type = local

## 数据输入
# 设置输入文件为input.fofn,该文件中包含有PacBio数据的fasta文件。此外也可以输入dexta格式的文件。dexta格式是Pacbio测序结果h5的一种压缩结果,能极大压缩测序数据文件的大小。可以使用FALCON软件附带的命令dexta将fasta文件转换为dexta格式。
input_fofn = input.fofn
# 设置输入数据为原始测序数据,或者为修正后的数据。
input_type = raw
#input_type = preads

## 对PacBio数据进行校正:对Pacbio raw subreads进行overlapping、consensus和pre-assembly分析,对reads进行校正。
# 选择长度高于指定阈值的reads进行分析
length_cutoff = 12000
# 选择长度高于指定阈值的reads进行预组装
length_cutoff_pr = 12000
# 调用fasta2DB命令将reads数据分割成多份Blocks。
# -s50 参数表示每份数据含有50Mbp的数据量,该参数默认值为200,默认参数情况下,使用daligner进行比对需要消耗约16G内存。
# -x500 参数表示read长度低于500bp的会被忽略掉。
# -a 参数表示来同一个自零级波导孔的第二个read不会被忽略掉。
pa_DBsplit_option = -x500 -s50
# 调用daligner对所有Blocks进行Overlapping分析。
# -v 参数表示输出daligner程序运行信息
# -B4 参数表示每个daligner命令对4个Blocks进行Overlapping分析。该参数的值越大,则每个daligner命令的计算量越大,但是总的任务数越少。该参数等同于以前的-dal
# -k -w -h 参数设置匹配的kmer相关参数,其默认值分别为 14,6,35 。
# -T4 表示每个daligner使用4个线程进行计算,该参数默认值是4,该参数可以设置成2,4,8,16,32...。
# -M32 表示每个daligner命令使用32G内存,加入该参数起到限制内存使用的作用,对大基因组比较有用。
# -t16 参数表示过滤掉覆盖度高于16的kmer,这些kmer会导致内存使用过多。默认设置下,daligner可以根据-M参数的值自动计算本参数的值。
# -l1000 参数表示忽略长度低于1000bp的reads。
# -s1000 参数表示记录比对结果时以每1000bp为一个记录点,相比于默认值100,能少很多记录点。
# 使用daligner的默认参数能很好地处理raw pacbio数据。
# 而对corrected pacbio数据,推荐使用-k20 -h60 -e.85参数。
pa_HPCdaligner_option = -v -B4 -k14 -T8 -t16 -e.70 -l1000 -s1000
# FALCON使用fc_consensus.py调用C语言写的程序来根据daligner进行Overlapping分析的结果进行consensus分析,从而对subreads进行校正。
# --min_cov 参数设置当read序列上某位点覆盖度低于指定阈值时,对read进行打断或截短,默认值是6。
# --min_cov_aln 参数设置当read序列的平均覆盖度低于指定阈值时,直接略过该read,默认值是10。
# --min_n_read 和 --max_n_read 参数设定比对结果中包含的reads数在此范围内才能得到consensus结果,其默认值分别是10和500。对于基因组重复程度较高的情况,要设置较低的--max_n_read值来减少对重复区域进行consensus分析的计算消耗。
# --min_idt 参数设置最小identity的比对结果能用于reads校正。
# --n_core 参数设置允许的线程数,默认值是24。
falcon_sense_option = --output_multi --min_idt 0.70 --min_cov 4 --max_n_read 200 --n_core 8
# 设置运行daligner任务的并发数。注意的是daligner和fc_consensus.py任务本身可以多线程,因此总的计算需求是线程数*并发数。若是总的计算需要远远超过服务器的计算资源,容易导致宕机。
da_concurrent_jobs = 10
la_concurrent_jobs = 2
# 设置运行fc_consensus.py任务的并发数。注意的是daligner和fc_consensus.py任务本身可以多线程,因此总的计算需求是线程数*并发数。若是总的计算需要远远超过服务器的计算资源,容易导致宕机。
cns_concurrent_jobs = 20
# 设置在SGE集群运行的并发数
sge_option_da = -pe smp 8 -q jobqueue
sge_option_la = -pe smp 2 -q jobqueue
sge_option_fc = -pe smp 80 -q jobqueue
sge_option_cns = -pe smp 16 -q jobqueue

## 对校正后的reads进行overlapping分析,其参数和上一个步骤的参数一致。
ovlp_DBsplit_option = -x500 -s50
ovlp_HPCdaligner_option = -v -B4 -k20 -h60 -T8 -t32 -e.96 -l500 -s1000
# 设置对校正后reads运行daligner任务的并发数。
pda_concurrent_jobs = 10
pla_concurrent_jobs = 2
sge_option_pda = -pe smp 8 -q jobqueue
sge_option_pla = -pe smp 2 -q jobqueue

## 过滤overlaps
# 若reads首尾两端的覆盖度比平均覆盖度大很多,则表明reads首尾是重复序列;若reads首尾两端的覆盖度比平均覆盖度相差较小很多,则表明reads首尾出现错误的可能性很大。需要过滤掉这种reads的overlaps结果。该步骤的参数设置不对,容易导致overlaps全部被过滤掉,得不到基因组组装的结果。
# --bestn设置报告reads此数目的最优overlaps。
# --min_cov和--max_cov表示所允许的reads首尾的覆盖度范围。对于通过length_cutoff得到 >20x 校正后的数据进行的基因组组装,可以设置--min_cov值为5,设置--max_cov为平均覆盖度的3倍。若用于基因组组装的数据量较少,可以设置该值为1或2。
# --max_diff设置所允许的首尾覆盖度值的最大差异。设置该参数的值为平均覆盖度的2倍。
# 对于特殊情况,可以设置更高的--max_cov和--max_diff值。
# 可以使用在1-preads_ovl目录下运行 fc_ovlp_stats.py --fofn merge-gather/las.fofn 导出overlap结果首尾的覆盖度结果,从而帮助设置以上参数
overlap_filtering_setting = --max_diff 50 --max_cov 75 --min_cov 5 --bestn 10

5. FALCON结果文件

0-rawreads/
    该目录存放对raw subreads进行overlpping分析与校正的结果;
    以job_前缀的文件夹的表示daligner任务的结果;
    以m_前缀的文件夹表示整合任务的结果,对raw subreads进行overlapping分析的比对结果是las文件,位于这些文件夹中;
    merge-gather中存放着las.fofn文件,该文件中存放着las文件的路径信息;
    preads目录下保存着校正后的reads信息,位于cns_前缀目录下的fasta文件。

1-preads_ovl/
    该目录存放对校正后reads进行overlapping的结果;
    以job_前缀的文件夹的表示daligner任务的结果;
    以m_前缀的文件夹表示整合任务的结果,对校正后reads进行overlapping分析的比对结果是las文件,位于这些文件夹中;
    merge-gather中存放着las.fofn文件,该文件中存放着las文件的路径信息;

2-asm-falcon/
    该目录是最终结果目录,包含draft contigs结果文件p_ctg.fa。

SGE部署

使用SGE,能使用户在单一的控制节点上投放多个任务,而不必考虑这些任务运行在哪个节点,能方便用户的使用。

修改主机名

对控制节点(master)的hostname进行修改:

修改配置文件 /etc/sysconfig/network 内容:

NETWORKING=yes
HOSTNAME=control

修改配置文件 /proc/sys/kernel/hostname 内容:

control

修改配置文件 /etc/hosts 内容:

127.0.0.1 localhost
192.168.30.1 control
192.168.30.2 node1

搭建NFS服务

NFS(Network File System)可以方便地将计算机上的指定文件夹共享给网络上的其它计算机。例如准备将 /share 目录共享

# mkdir /share
# chmod 1777 /share

通过修改配置文件 /etc/exports 来共享指定的文件夹。在该配置文件中添加以下一行内容:

/share	192.168.30.0/24(rw,sync,no_root_squash,no_subtree_check)

当前主机的IP是192.168.30.1,以上配置信息表示将当前主机的 /share 目录共享给同一个局域网内指定IP段的计算机,且其权限如下:

rw
  可读可写。
ro
  只读。
sync
  将数据同步写入到内存和磁盘中。
async
  将数据会先暂存于内存中,必要时才写入磁盘。
no_root_squash
  若客户端使用root用户操作共享文件夹的时候,具有最大权限。
root_squash(默认)
  若客户端使用root用户操作共享文件夹的时候,将其身份设定成匿名用户nfsnobody,降低权限。
no_all_squash(默认)
  访问用户先与本机用户匹配,匹配失败后再映射为匿名用户或用户组。
all_squash
  客户端的使用者用户统一被转换成匿名用户nfsnobody。
subtree_check (默认)
  若输出目录是一个子目录,则nfs服务器将检查其父目录的权限。
no_subtree_check
  即使输出目录是一个子目录,nfs服务器也不检查其父目录的权限,这样可以提高效率。

由于NFS服务的目的是能对多台服务端计算机共享nfs服务器的指定目录,它需要随机取用多个小于1024的端口来传输数据。而这些端口的开放是不固定的,为了让客户端能连接上正确的端口,则需要开启nfs服务的同时开启RPC(Remote Procedure Call)服务。
NFS服务的开启需要设置开放一些端口,可以在 /etc/sysconfig/nfs 中看到,修改该配置文件,将所有端口设置行前的注释取消掉

RQUOTAD_PORT=875
LOCKD_TCPPORT=32803
LOCKD_UDPPORT=32769
MOUNTD_PORT=892
STATD_PORT=662
STATD_OUTGOING_PORT=2020
RDMA_PORT=20049

此外,NFS服务本身的端口是2049,RPC服务的端口是111,因此需要在防火墙中开放以上9个端口。修改防火墙配置文件/etc/sysconfig/iptables,在正确位置添加:

-A INPUT -p tcp -s 192.168.30.0/24 -m multiport --dport 111,2049,875,32803,32769,892,662,2020,20049 -j ACCEPT
-A INPUT -p udp -s 192.168.30.0/24 -m multiport --dport 111,2049,875,32803,32769,892,662,2020,20049 -j ACCEPT

最后,启动RPC和NFS服务:

# /etc/init.d/rpcbind restart
# /etc/init.d/nfs restart

# chkconfig rpcbind on
# chkconfig nfs on

在客户端192.168.30.2计算机上使用nfs服务器的共享文件夹

# mkdir /share
# mount -t nfs 192.168.30.1:/share /share/

以上通过NFS共享了 /share 目录,该目录在nfs服务器上的权限和服务端的权限是一致的,其权限是根据UID来识别的,因此客户端和服务端有共同的用户名和UID,且用户名和UID是完全匹配的,才有利于文件的共享。当机器数目很多的时候,为了能保证多台计算机上具有相同的用户名和密码等设置,则可以使用NIS服务来解决。

若在Master机器上修改了/etc/exports配置文件信息,则使用命令exportfs命令使修改生效

# exportfs -rv

NIS服务

在Master和Slaves机器上都安装NIS(Network Information Service)软件

# yum install ypserv

在Master和Slaves机器上修改/etc/sysconfig/network,尾部添加如下两行,使它们具有相同的NIS域名。

NISDOMAIN=chenlianfuNIS
YPSERV_ARGS="-p 1011"

以上第一行的值可以随意设定,但要求Master和Slaves机器上的该值一致
以上第二行表示NIS启动在1011端口上

在NIS配置文件 /etc/ypserv.conf 尾部添加一行来设定Slaves机器的权限

*                         : * : * : none

以上设置表示所有机器都都有最大权限。虽然限制不严格,但是可以通过iptables防火墙来进行安全控制。

修改 /etc/hosts,将Master和Slaves机器上的配置文件都设置正确

127.0.0.1 localhost
192.168.30.1 control
192.168.30.2 node1

修改 /etc/sysconfig/yppasswdd 配置文件,设置开放端口

YPPASSWDD_ARGS="--port 1012"

修改 /etc/sysconfig/iptables 防火墙配置,开放1011和1012端口

-A INPUT -p tcp -s 192.168.30.0/24 -m multiport --dport 111,2049,875,32803,32769,892,662,2020,20049,1011,1012 -j ACCEPT
-A INPUT -p udp -s 192.168.30.0/24 -m multiport --dport 111,2049,875,32803,32769,892,662,2020,20049,1011,1012 -j ACCEPT

启动服务

# /etc/init.d/iptables restart
# /etc/init.d/ypserv start
# /etc/init.d/yppasswdd start

ypserv用于启动NIS服务
yppasswdd用于启动NIS客户端密码修改服务

# chkconfig ypserv on
# chkconfig yppasswdd on

将Master机器上账号转换成数据库

# /usr/lib64/yp/ypinit -m

进入交互式界面,直接按control + d结束选择,再按y同意,程序则根据Master机器内的用户来创建数据库。

在Slaves机器上操作来启动NIS服务

# setup

进入交互式界面,选择Authentication configuration, 使用TAB切换并用空格选中Use NIS,使用TAB切换选择Next, 填写正确的NISDomain和NIS服务器的IP,使用TAB切换选择OK,选择退出后Slave机器会开启NIS服务。

在Slave机器上成功开启NIS服务后,则可以使用Master机器上的用户名远程登录Slave机器了。若该用户名在Slave机器中存在,则直接用该Slave机器中的用户直接登录;若该用户名在Slave机器中不存在,则会使用NIS服务,使用Master机器中的用户登录,但是该用户没家目录及其配置文件。

若需要更多用户用于NIS服务,则在Master机器中新建用户用,重新运行ypinit命令进行初始化即可。

使用yptest命令检验数据库信息

# yptest

在Test 9中会给出NIS Master机器上的所有账户信息。如果给出信息正常,表示验证成功。

使用ypwhich检验数据库文件

# ypwhich -x

可以看到相关文件名,这些文件名存放在 /var/yp/chenlianfuNIS/ 目录下。

使用ypcat读取数据库文件内容

# ypcat passwd.byname

使用yppasswd在Slave机器上修改Master机器上的用户名,其用法和passwd用法一致。

# yppasswd

SGE

先开启SGE需要的端口6444,修改配置文件/etc/sysconfig/iptables

-A INPUT -p tcp -s 192.168.30.0/24 -m multiport --dport 111,2049,875,32803,32769,892,662,2020,20049,1011,1012,6444 -j ACCEPT
-A INPUT -p udp -s 192.168.30.0/24 -m multiport --dport 111,2049,875,32803,32769,892,662,2020,20049,1011,1012,6444 -j ACCEPT

重启iptables服务

# /etc/init.d/iptables restart

在Master机器和Slaves机器上都使用yum安装gridengin软件

需要使用epel源安装gridengin
# rpm -Uvh http://dl.fedoraproject.org/pub/epel/epel-release-latest-6.noarch.rpm
# yum install gridengine-*

在Master机器上部署SGE

# cd /usr/share/gridengine/
# ./install_qmaster
进入交互式界面,基本全部Enter即可,需要输入密码,我就选择123456。
在Master上,这一个命令同时部署了控制进程和执行进程。

在Slaves机器上部署SGE,需要先将Master机器上的/usr/share/gridengine/文件夹拷贝到Slaves机器上的相同路径上。

在Master机器上进行打包操作,注意软链接的拷贝。
# cd /usr/share/gridengine
# mkdir /share/gridengine
# cp -aL * /share/gridengine
# cd /share/
# rm gridengine/ -rf
# tar zcf gridengine.tar.gz gridengine/

在Slave机器上下载压缩包,然后解压缩覆盖其/usr/share/gridengine/文件夹
# cd /usr/share
# tar zxf /share/gridengine.tar.gz

在Slave机器上部署SGE,仅需要部署执行进程。
# cd gridengine
# ./install_execd
进入交互式界面,全部enter即可。

启动SGE

在部署SGE的时候已经启动了SGE服务
启动控制进程
# /etc/init.d/sgemaster restart
启动执行进程
# /etc/init.d/sge_execd restart

# chkconfig sgemaster on
# chkconfig sge_execd on

在控制节点上同时启动以上2个进程,而在其它计算节点上仅启动执行进程。

SGE的使用原理:

集群中的主机分2种:控制节点(mater)和计算节点(slave)。其中控制节点只在一台机器上部署,该控制节点也同时作为计算节点;其它主机全部是计算节点。

计算资源是由host的slots构成。可以选取集群中部分的hosts,定义为host用户组。

队列则表示集群中计算资源的容器。例如,名称叫all.q的队列对应着集群中全部的计算资源。若不想让某些用户使用集群全部的计算资源,则定义一个新的队列名,且该队列仅能使用集群部分的计算资源。

使用SGE集群进行计算的时候,为了进行并行化计算,需要设置并行化参数。

SGE的使用

qconf -ae hostname
    添加执行主机
qconf -de hostname
    删除执行主机
qconf -sel
    显示执行主机列表

qconf -ah hostname
    添加管理主机
qconf -dh hostname
    删除管理主机
qconf -sh
    显示管理主机列表

qconf -as hostname
    添加提交主机
qconf -ds hostname
    删除提交主机
qconf -ss
    显示提交主机列表

qconf -ahgrp groupname
    添加主机用户组
qconf -mhgrp groupname
    修改主机用户组
qconf -shgrp groupname
    显示主机用户组成员
qconf -shgrpl
    显示主机用户组列表

qconf -aq queuename
    添加集群队列
qconf -dq queuename
    删除集群队列
qconf -mq queuename
    修改集群队列配置
qconf -sq queuename
    显示集群队列配置
qconf -sql
    显示集群队列列表

qconf -ap PE_name
    添加并行化环境
qconf -mp PE_name
    修改并行化环境
qconf -dp PE_name
    删除并行化环境
qconf -sp PE_name
    显示并行化环境
qconf -spl
    显示并行化环境名称列表

qstat -f
    显示执行主机状态
qstat -u user
    查看用户的作业
qhost
    显示执行主机资源信息

通过使用命令qconf -mq queuename来对队列进行配置。修改hostlist来配置该队列可以使用执行主机;修改slots来配置各台执行主机可使用的线程数。从而对队列的计算资源进行设置。
部署完毕SGE后,会生成一个默认主机用户组@allhosts,它包含所有的执行节点;生成一个默认的all.q队列名,它包含所有节点所有计算资源。默认的队列包含的计算资源是最大的。

使用qsub提交作业

qsub简单示例:
$ qsub -V -cwd -o stdout.txt -e stderr.txt run.sh

其中run.sh中包含需要运行的程序,其内容示例为如下三行:
#!/bin/bash
#$ -S /bin/bash
perl -e 'print "abc\n";print STDERR "123\n";'

qsub的常用参数:
-V
    将当前shell中的环境变量输出到本次提交的任务中。
-cwd
    在当前工作目录下运行程序。默认设置下,程序的运行目录是当前用户在其计算节点的家目录。
-o
    将标准输出添加到指定文件尾部。默认输出文件名是$job_name.o$job_id。
-e
    将标准错误输出添加到指定文件尾部。默认输出文件名是$job_name.e$job_id。
-q
    指定投递的队列,若不指定,则会尝试寻找最小负荷且有权限的队列开始任务。
-S
    指定运行run.sh中命令行的软件,默认是tcsh。推荐使用bash,设置该参数的值为 /bin/bash 即可,或者在run.sh文件首部添加一行#$ -S /bin/bash。若不设置为bash,则会在标准输出中给出警告信息:Warning: no access to tty (Bad file descriptor)。
-hold_jid
    后接多个使用逗号分隔的job_id,表示只有在这些job运行完毕后,才开始运行此任务。
-N
    设置任务名称。默认的job name为qsub的输入文件名。
-p
    设置任务优先级。其参数值范围为 -1023 ~ 1024 ,该值越高,越优先运行。但是该参数设置为正数需要较高的权限,系统普通用户不能设置为正数。
-j y|n
    设置是否将标准输出和标准错误输出流合并到 -o 参数结果中。
-pe
    设置并行化环境。

任务提交后的管理

$ qstat -f
    查看当前用户在当前节点提交的所有任务,任务的状态有4中情况:qw,等待状态,刚提交任务的时候是该状态,一旦有计算资源了会马上运行;hqw,该任务依赖于其它正在运行的job,待前面的job执行完毕后再开始运行,qsub提交任务的时候使用-hold_jid参数则会是该状态;Eqw,投递任务出错;r,任务正在运行;s,被暂时挂起,往往是由于优先级更高的任务抢占了资源;dr,节点挂掉后,删除任务就会出现这个状态,只有节点重启后,任务才会消失。

$ qstat -j jobID
    按照任务id查看

$ qstat -u user
    按照用户查看

$ qdel -j jobID
    删除任务

使用openmpi来支持SGE的并行化环境,首先使用–with-sge参数来安装openmpi。

$ tar zxf openmpi-1.8.6.tar.gz
$ cd openmpi-1.8.6
$ ./configure --prefix=/opt/sysoft/openmpi-1.8.6 --with-sge
$ make -j 4
$ make install
$ cd .. && rm -rf openmpi-1.8.6/
$ echo 'export PKG_CONFIG_PATH=/opt/sysoft/openmpi-1.8.6/lib/pkgconfig/:$PKG_CONFIG_PATH
export LD_LIBRARY_PATH=/opt/sysoft/openmpi-1.8.6/lib/:$LD_LIBRARY_PATH
export C_INCLUDE_PATH=/opt/sysoft/openmpi-1.8.6/include:$C_INCLUDE_PATH
export PATH=/opt/sysoft/openmpi-1.8.6/bin/:$PATH' >> ~/.bashrc.openmpi
$ source ~/.bashrc.openmpi

然后,将openmpi的安装结果完全复制到所有节点的相同路径下。推荐使用NFS来搞定。

再添加并行化环境

# qconf -ap mpi

进入vi编辑界面,修改其中两项:
slots              200
allocation_rule    $fill_up
然后:wq保存退出

运行qsub则可以使用并行化环境(即可以使用-pe参数)了。

使用HTSeq进行有参转录组的表达量计算

1. HTSeq简介

HTSeq是使用Python编写的一支用于进行基因Count表达量分析的软件,能根据SAM/BAM比对结果文件和基因结构注释GTF文件得到基因水平的Counts表达量。HTSeq进行Counts计算的原理非常简单易懂,容易上手。

2. HTSeq安装

PYPI下载HTSeq的Python包
$ wget https://pypi.python.org/packages/46/f7/6105848893b1d280692eac4f4f3c08ed7f424cec636aeda66b50bbcf217e/HTSeq-0.7.2.tar.gz
$ tar zxf HTSeq-0.7.2.tar.gz
$ cd HTSeq-0.7.2
$ /opt/sysoft/Python-2.7.11/bin/python setup.py build
$ /opt/sysoft/Python-2.7.11/bin/python setup.py install
$ cd ../ && rm HTSeq-0.7.2 -rf

3. HTSeq使用

3.1 HTSeq的Count模式

HTSeq计算counts数有3种模式,如下图所示(ambiguous表示该read比对到多个gene上;no_feature表示read没有比对到gene上):
HTSeq Count模式

3.2 HTSeq的使用命令

HTseq安装完毕后,在Python软件的bin目录下生成htseq-count命令。
htseq-count运行简单示例:

对于非链特异性真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于链特异性测序真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s reverse -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于非链特异性原核生物转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m intersection-strict bowtie2.sam genome.gtf > counts_out.txt

htseq-count命令的常用参数:

-f | --format     default: sam
  设置输入文件的格式,该参数的值可以是sam或bam。
-r | --order     default: name
  设置sam或bam文件的排序方式,该参数的值可以是name或pos。前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。
-s | --stranded     default: yes
  设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序;若是单端测序,yes表示read比对到了基因的正义链上;若是双末端测序,yes表示read1比对到了基因正义链上,read2比对到基因负义链上;reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解,一般情况下双末端链特异性测序,该参数的值应该选择reverse(本人暂时没有测试该参数)。
-a | --a     default: 10
  忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
-t | --type     default: exon
  程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。
-i | --idattr     default: gene_id
  设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。
-m | --mode     default: union
  设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。
-o | --samout 
  输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。
-q | --quiet
  不输出程序运行的状态信息和警告信息。
-h | --help
  输出帮助信息。

3.3 HTSeq使用注意事项

HTSeq的使用有如下注意事项,否则得到的结果是错误的:

1. HTSeq是对有参考基因组的转录组测序数据进行表达量分析的,其输入文件必须有SAM和GTF文件。
2. 一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析,而不是一个样品内部基因的表达量比较。因此,HTSeq设置了-a参数的默认值10,来忽略掉比对到多个位置的reads信息,其结果有利于后续的差异分析。
3. 输入的GTF文件中不能包含可变剪接信息,否则HTSeq会认为每个可变剪接都是单独的基因,导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous,从而不能计算到基因的count中。即使设置-i参数的值为transcript_id,其结果一样是不准确的,只是得到transcripts的表达量。

3.4 HTSeq的结果

HTSeq将Count结果输出到标准输出,其结果示例如下:

gene00001	0
gene00002	9224
gene00003	880
...
gene12300	1043
gene12301	200
__no_feature	127060
__ambiguous	0
__too_low_aQual	4951
__not_aligned	206135
__alignment_not_unique	0

使用DBG2OLC对二、三代混合数据进行基因组组装

1. DBG2OLC软件简介

DBG2OLC能利用二代和三代混合数据组装大基因组。其文章于2016年发表在Scientific Reports上。

2. DBG2OLC软件下载与安装

使用git下载DBG2OLC软件

$ cd /opt/biosoft/
$ git clone https://github.com/yechengxi/DBG2OLC.git
$ cd /opt/biosoft/DBG2OLC
按照说明中对软件进行编译,编译出的3个可执行程序全部都是DBG2OLC命令
$ g++ -O3 -o SparseAssebmler DBG2OLC.cpp
$ g++ -O3 -o DBG2OLC *.cpp
$ g++ -O3 -o Sparc *.cpp
直接拷贝作者编译好的程序即可
$ chmod 755 compiled/*
$ cp compiled/* .
$ echo 'PATH=$PATH:/opt/biosoft/DBG2OLC' >> ~/.bashrc
$ source ~/.bashrc

DBG2OLC程序第三步需要blasr, sparc/pbdagcon软件支持。其中sparc在DBG2OLC安装文件夹下。
安装blasr

下载BLASR
$ git clone https://github.com/PacificBiosciences/blasr.git /opt/biosoft/blasr
$ cd /opt/biosoft/blasr/
下载libcpp和pbbam两个submodules
$ make update-submodule

blasr编译需要hdf5支持,从hdf5官网下载适合centos6的二进制包并安装 
$ wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0-patch1/bin/linux-centos6-x86_64-gcc447/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz
$ tar zxf hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz -C /opt/sysoft/

可以使用cmake, pitchfork和make三种方式对BLASR进行编译,以下使用常规的make方法进行编译,需要高版本gcc支持
对BLASR进行configure
$ ./configure.py --shared --sub --no-pbbam HDF5_INCLUDE=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/include/ HDF5_LIB=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/
对submodules进行configure
$ make configure-submodule
对submodules进行make
$ make build-submodule -j 4
对BLASR进行make
$ make blasr -j 4
对其它工具,例如pls2fasta, loadPulses, sawriter等进行编译,其结果文件在utils文件夹中
$ make -j 4

可选步骤:手动将有用的命令和库文件放置到指定的地方
blasr的正常运行需要依赖libcpp里面三个库文件和hdf5软件中的库文件
$ mkdir bin lib
$ cp blasr bin/
$ find utils -maxdepth 1 -perm 775 -exec cp {} bin/ \;
$ cp ./libcpp/pbdata/libpbdata.so ./libcpp/hdf/libpbihdf.so ./libcpp/alignment/libblasr.so lib/
$ echo 'export LD_LIBRARY_PATH=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/:/opt/biosoft/blasr/lib/:$LD_LIBRARY_PATH
PATH=/opt/biosoft/blasr/bin/:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

若DBG2OLC流程第三步选择使用pbdagcon进行运算,则需要安装pbdagcon软件

pbdagcon软件的编译需要高版本gcc支持
$ git clone https://github.com/PacificBiosciences/pbdagcon.git /opt/biosoft/pbdagcon
$ cd /opt/biosoft/pbdagcon
$ ./configure.py --boost --gtest --sub --no-pbbam
$ make init-submodule
$ make -j 4
$ make check
$ mkdir bin
$ cp src/cpp/dazcon src/cpp/pbdagcon bin/
$ echo 'PATH=/opt/biosoft/pbdagcon/bin:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

3. 程序运行

使用DBG2OLC软件利用二代和三代数据混合的基因组组装,其运行流程分3步。

3.1 使用SparseAssembler利用二代数据进行DBG组装

首先,利用Illumina小片段文库数据使用SparseAssembler命令组装出contigs序列。此外,也可以使用其他基因组组装软件组装出contigs序列后,直接跳到DBG2OLC的第二个步骤。值得注意的是:输入到第二步骤的contigs必须是没有经过repeat resolving的原始序列;绝大部分基因组组装软件为了获得更完整连续的contigs序列,都牺牲了contigs的准确性,其结果不能用于DBG2OLC软件的第二步,否则最终结果会很差;作者推荐可以直接用于第二步的其它contig组装软件有Platanus和Meraculous。
一般情况下,输入到SparseAssembler命令中~50x的Illumina数据,能获得较好的contigs结果。
SparseAssembler命令参数:

常用参数:
LD 
    是否载入k-mer graph。第一次运行SparseAssembler命令的时候,该参数的值必须是0;若为了使用SparseAssembler得到更好的contigs结果,则需要调整NodeCovTh和EdgeCovTh参数;调整这些参数的时候,不需要再次计算k-mer graph,设置该参数为1来跳过这个步骤,从而节约很多时间。
k 
  设置使用DBG方法计算时的Kmer大小,支持的Kmer大小为15-127。
g 
  number of skipped intermediate k-mers, support 1-25.该参数软件示例中使用的值是15。
NodeCovTh 
  设置k-mers覆盖度阈值,去除覆盖度较低的k-mers。该值设定范围为0-16,默认值为1。
EdgeCovTh 
  设置link覆盖度阈值,去除覆盖度较低的links。该值设定范围为0-16,默认值为0。
GS 
  设置一个基因组大小的值。该参数用于决定预先占用的内存量。推荐设置得比基因组大,例如设置为2倍基因组大小。
f 
  输入单端测序数据的路径。输入文件可以是fasta或fastq文件。若有多个输入文件,则使用多个f参数。
i1  i2 
  输入inward paired-end数据。若有多组paired-end数据,则多次使用i1/i2参数。

其它参数:
o1  o2 
  输入outward paired-end数据。
i1_mp  i2_mp 
  输入插入片段长度>10kb的inward paired-end数据。
o1_mp  o2_mp 
  输入插入片段长度>10kb的outward paired-end数据。
PathCovTh 
  设置path覆盖度阈值,去除覆盖度较低的paths。该值设定范围为0-100。根据经验,不推荐添加该参数。
TrimLen 
  将所有过长的序列截短到此指定的长度。
TrimN 
  若read中的碱基N数目超过此设定的值,则去除该read数据。
TrimQual 
  从尾部截短质量低于此值的碱基。
QualBase 
  设置Fastq文件中最低碱基质量对应的ASCII码符号。默认值是'!',等同于Pred33。
Denoise 
  设置是否对reads进行修正。默认值是0,表示不对reads进行修正。
H 
  混合模式。默认值是0,表示对reads的尾部进行截短,来保证高质量的数据进行reads修正。
CovTh 
  覆盖度 < 此设定值的k-mer会被检测,从而被校正。若该参数值设置为0,则软件会自动计算该值。
CorrTh 
  覆盖度 >= 此设定值的k-mer可以用来对reads做校正。若该参数值设置为0,则软件会自动计算该值。

SparseAssembler运行示例:

对某物种Illumina小片段文库测序的PE150bp数据使用trimmomatic质控,再使用FindErrors进行修正,再运行SparseAssembler:
$ SparseAssembler LD 0 k 95 g 15 NodeCovTh 1 EdgeCovTh 0 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq
$ cp Contigs.txt Contigs.txt.00
增大NodeCovTh和EdgeCovTh参数后,再次运行SparseAssembler,并比较两次结果。第二次运行较第一次运行,耗时少了很多很多。
$ SparseAssembler LD 1 k 95 g 15 NodeCovTh 2 EdgeCovTh 1 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq

SparseAssembler在当前目录下生成了18个文件结果,其中Contigs.txt文件是Fasta格式的Contigs序列。
运行SparseAssembler的注意事项:

1. SparseAssembler只可以简单地对Fastq文件进行质控和错误修正。推荐使用其它软件进行reads质控和修正,以获得更好的结果。
2. 参数k设置了k-mer的大小,该参数的值对结果的影响较大。若基因组较小,推荐设置多个k-mer值进行多次计算,从而选择最优k-mer值。个人经验,PE150bp数据的最优的k-mer值约为91~99。
3. 选定了k-mer大小后,使用默认的NodeCovTh和EdgeCovTh参数(默认参数一般能得到很好的结果)运行一遍SparseAssembler。然后尝试增大NodeCovTh和EdgeCovTh参数,设置LD 1参数再次运行SparseAssembler,以获得最优的Contigs结果。
4. 可能是先使用了最小的NodeCovTh和EdgeCovTh参数做运算后,才能再次使用更大的参数进行运算。
5. SparseAssembler虽然也有输入大片段文库数据的参数和Scaffolding参数,但是不推荐输入大片段文库数据进行Scaffolding操作,没太大意义。
6. 虽然SparseAssembler命令的文件输入方式有多种,若是仅进行contigs组装,没有利用到paired信息,因此使用i1 i2参数输入文件和使用f参数输入文件的结果是一模一样的。

3.2 使用DBG2OLC找Contigs序列和Pacbio reads的Overlap并进行Layout

DBG2OLC通过比较contigs和Pacbio reads之间的overlap,将contigs序列定位到Pacbio reads上,将DBG的contigs结果运用到OLC算法中。
DBG2OLC命令参数:

主要参数:
LD 
  是否载入compressed reads information。第一次运行DBG2OLC命令的时候,该参数的值必须是0;若为了得到更好的结果,则需要调整其它参数;调整这些参数的时候,设置该参数为1来跳过这个步骤,从而节约很多时间。
k 
  设置k-mer大小。k-mer用来比较contig和pacbio read之间的重叠,而不是用于基因组组装,推荐设置为 17 即可。
AdaptiveTh 
  若contig和pacbio read之间匹配的k-mers数目 < AdaptiveTh * contig长度,则认为contig和pacbio read没有重叠。推荐设置为0.001-0.02。
KmerCovTh 
  若contig和pacbio read之间匹配k-mers的覆盖度 < KmerCovTh,则认为contig和pacbio read没有重叠。推荐设置为2-10。
MinOverlap 
  两条Pacbio read之间匹配的k-mers数目 < MinOverlap,则认为它们之间没有重叠。推荐设置为10-150。
RemoveChimera 
  去除嵌合体Pacbio reads。若Pacbio数据覆盖度大于10x,推荐设置该参数为 1 。
Contigs 
  输入contigs序列文件路径
f 
  输入Pacbio测序Fasta/Fastq文件路径。

其它参数:
MinLen 
  设置能用于计算的最小Pacbio reads长度。
ChimeraTh 
  该参数默认值是 1 ;若Pacbio数据覆盖度大于100x,则推荐加入该参数并设置为 2 。
ContigTh 
  该参数默认值是 1 ;若Pacbio数据覆盖度大于100x,则推荐加入该参数并设置为 2 。

DBG2OLC运行示例:

$ DBG2OLC LD 0 k 17 AdaptiveTh 0.001 KmerCovTh 2 MinOverlap 20 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq
$ DBG2OLC LD 1 k 17 AdaptiveTh 0.005 KmerCovTh 3 MinOverlap 30 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq

DBG2OLC的结果文件很多,其主要结果文件是backbone_raw.fasta和DBG2OLC_Consensus_info.txt,是第三步的输入文件。
运行DBG2OLC的注意事项:

1. AdaptiveTh, KmerCovTh和minOverlap这3个计算Overlap的参数对结果的影响最大。对于10x/20x PacBio数据:KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01;对于50x-100x PacBio数据:KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02。
2. 不推荐对Pacbio数据就行修正后再运行DBG2OLC。可以比较使用修正前的数据用于DBG2OLC的结果,一般情况下使用未修正的Pacbio数据能获得更好的结果。此外,DBG2OLC运行过程中有一步多序列比对模块来进行错误修正。
3. 可能是先使用了最小的AdaptiveTh, KmerCovTh和minOverlap参数做运算后,才能再次使用更大的参数进行运算。

3.3 Call consensus

本步骤需要使用/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本,来调用blasr和consensus模块Sparc(可以考虑使用pbdagcon作为consensus模块,但DBG2OLC没有提供相应的脚本)进行运算。

先将/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本拷贝到当前目录
$ cp /opt/biosoft/DBG2OLC/utility/*.sh /opt/biosoft/DBG2OLC/utility/*.py ./
若使用了最新版本的blasr软件,其参数书写方法有一个中划线变成了两个中划线,因此需要修改.sh文件中blasr命令的参数书写方法。
此外,也需要修改.sh文件中Sparc命令路径,或者将Sparc命令也拷贝到当前目录。

将Contigs序列和Pacbio reads数据合并成一个文件ctg_pb.fasta
$ cp Contigs.txt ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell01.fastq >> ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell02.fastq >> ctg_pb.fasta

运行脚本程序split_and_run_sparc.sh
$ ./split_and_run_sparc.sh backbone_raw.fasta DBG2OLC_Consensus_info.txt ctg_pb.fasta ./ 2 > cns_log.txt
结果会输出到 ./ 目录下。最后的结果文件是final_assembly.fasta。

使用AWStats对网站流量进行统计

1. 安装AWStats

# yum install awstats

2. 使用AWStats

常用示例:
# /var/www/awstats/awstats.pl --config=localhost.localdomain -update -output > /var/www/awstats/index.html

--config=virtualhostname
    该参数用于导入配置文件。配置文件位于 /etc/awstats 或 /usr/local/etc/awstats 目录。程序会导入awstats.virtualhostname.conf或awstats.conf配置文件。
-update
    对数据统计结果进行更新。
-output
    输出HTML结果文件。

对Perl代码进行编译与加密

我写了一些Perl程序。为了防止程序在传播扩散过程中遭人随意篡改或出售而引起版权纠纷,于是需要对一些程序进行编译和加密处理。

1. 使用perlcc命令对perl代码进行编译

我安装的时CentOS 6系统,该系统中默认能的Perl版本时5.10版本。该版本中取消了perlcc命令以及相应的B::C, B::CC, B::Bytecode等模块。只有不搞于5.9.4版本的perl才会有perlcc命令。详情请见:http://perldoc.perl.org/perl5100delta.html。
perlcc的使用方法(http://search.cpan.org/~nwclark/perl-5.8.9/utils/perlcc.PL):

$ perlcc -o hello hello.pl

2. 使用pp命令对perl代码进行编译

perlcc命令可能对perl代码编译不成功,或成功后不能正常运行。推荐使用pp命令来进行该项工作。

2.1 安装 pp 和 PAR::Filter::Crypto 模块

pp 模块用于perl程序的编译和打包:http://search.cpan.org/~rschupp/PAR-Packer-1.035/lib/pp.pm
PAR::Filter::Crypto 模块用于加密: http://search.cpan.org/~shay/Filter-Crypto-2.07/lib/PAR/Filter/Crypto.pm

$ sudo cpan -i pp
$ sudo cpan -i PAR::Filter::Crypto

2.2 对perl程序进行编译和加密

$ pp -f Crypto -F Crypto -M Filter::Crypto::Decrypt -o hello hello.pl

CentOS系统NAT共享上网

现在有服务器A通过PPPOE联网,而服务器B直接和服务器A使用网线连接。若需要使B能正常上网,则需要将A的网络共享给B。此外,服务器A和B都具有多网口,并都是CentOS系统。将服务器A的网络共享给B,其对两台服务器的设置如下:

1. 服务器A的设置

1.1 服务器A的第一个网口进行pppoe连接

将网线插入到服务器A的第一个网口eth0,然后设置服务器A的PPPOE连接:

安装pppoe软件
# yum install rp-pppoe
配置pppoe设置,填写上网账号和密码,该pppoe配置名称为ppp0,保证对应的网口为eth0,设置。
# pppoe-setup
关闭ppp0的连接
# ifdown ppp0
开启ppp0的连接
# ifup ppp0
若发现ppp0连接不上,输入下面命令后再连接
# pppoe-stop

1.2 服务器A的第二个网口的IP设置

再将服务器A的eth1口和服务器B的eth1口进行连接。对服务器A的eth1口进行设置:

# setup
通过setup命令配置eth1的IP,设置其:
IP地址:192.168.1.1
子网掩码:255.255.255.0
网关:192.168.1.1
DNS1:211.69.143.1
DNS2:8.8.8.8

其中DNS1是我们学校提供的DNS服务器网址,DNS2是google提供的DNS网址。可能在不同的地方其DNS网址不一样。
然后,修改 eth0 和 eth1 的配置文件,设置这两个网口开机启动:
# vi /etc/sysconfig/network-scripts/ifcfg-eth0
# vi /etc/sysconfig/network-scripts/ifcfg-eth1
ONBOOT=yes

1.3 将服务器A的ppp0网络进行NAT共享

首先,修改配置文件/etc/sysctl.conf的一个参数来开启使用NAT进行IP转发。

$ perl -p -i -e 's/net.ipv4.ip_forward =.*/net.ipv4.ip_forward = 1/' /etc/sysctl.conf

再使用iptables命令来对ppp0进行网络共享。

生成文件 /usr/local/bin/ishare 并使其可执行,执行该命令,即可共享 ppp0 网络。
# echo '#!/bin/bash
## Internet connection shating script
sysctl -w net.ipv4.ip_forward=1
sysctl -p
iptables -X
iptables -F
iptables -t nat -X
iptables -t nat -F
iptables -I INPUT -m state --state RELATED,ESTABLISHED -j ACCEPT
iptables -I FORWARD -m state --state RELATED,ESTABLISHED -j ACCEPT
iptables -t nat -I POSTROUTING -o ppp0 -j MASQUERADE' > /usr/local/bin/ishare
# chmod 755 /usr/local/bin/ishare
# /usr/local/bin/ishare

若需要开机则启动该命令:
# echo '/usr/local/bin/ishare' >> /etc/rc.d/rc.local

1.4 重启服务器A的网络设置

服务器A的配置修改完毕,然后重启服务器A的网络,使配置生效:

# /etc/init.d/network restart
保证看到pppoe和eth1网络的正常启动。

2. 服务器B的设置

配置服务器B的 eth1 网口:

# setup
通过setup命令配置eth1的IP,设置其:
IP地址:192.168.1.2
子网掩码:255.255.255.0
网关:192.168.1.1
DNS1:211.69.143.1
DNS2:8.8.8.8

修改 eth1 的配置文件,设置这该网口开机启动:
# vi /etc/sysconfig/network-scripts/ifcfg-eth1
ONBOOT=yes

再重启重启服务器B的网络设置
# /etc/init.d/network restart
保证看到eth1网络的正常启动。

3. 会出现的问题

以上设置完毕后,通过 ping 命令来检测服务器B是否能联网。若服务器ping不通192.168.1.1,则表示服务器B和服务器A连接不通;若能ping通192.168.1.1而不能联外网,则表示服务器A没有开启共享或服务器A没有联网。需要按照上述教程逐步排查。

mycology_学科杂志分区

生物大学科分区	MYCOLOGY小学科分区	杂志名称	2014-2015影响因子
1区	1区	STUDIES IN MYCOLOGY	13.25
2区	2区	FUNGAL DIVERSITY	6.221
2区	2区	PERSOONIA	5.3
3区	2区	MYCORRHIZA	3.459
3区	3区	Fungal Ecology	2.929
3区	3区	FUNGAL GENETICS AND BIOLOGY	2.587
4区	3区	MYCOLOGIA	2.471
4区	4区	Fungal biology	2.342
3区	3区	MEDICAL MYCOLOGY	2.335
4区	4区	MYCOSES	2.239
3区	3区	World Mycotoxin Journal	2.157
4区	4区	MYCOLOGICAL PROGRESS	1.913
4区	4区	MYCOPATHOLOGIA	1.528
4区	4区	CRYPTOGAMIE MYCOLOGIE	1.524
4区	4区	LICHENOLOGIST	1.454
4区	4区	Mycoscience	1.418
4区	4区	Revista iberoamericana de micología : órgano de la Asociación Espa ola de Especialistas en Micología1.056
4区	4区	SYDOWIA	1.021
4区	4区	INTERNATIONAL JOURNAL OF MEDICINAL MUSHROOMS	0.927
4区	4区	MYCOTAXON	0.705
4区	4区	JOURNAL DE MYCOLOGIE MEDICALE	0.573

blast进行重复序列屏蔽

1. 构建数据库的时候屏蔽参考序列的重复

segmasker 屏蔽氨基酸的低复杂序列
dustmasker 屏蔽核算序列的低复杂序列
windowmasker 按照序列重复的次数来屏蔽
convert2blastmask 根据小写字母来屏蔽

这几个都可以先得到一个含有屏蔽信息的文件。然后进行 makeblastdb 的时候输入这个文件,就可以相应的 masked 数据库了。

参考:http://www.ncbi.nlm.nih.gov/books/NBK279681/

2. 比对的时候对query序列的重复进行屏蔽

blast 比对的时候,可以对 query 序列进行屏蔽。 这几个参数估计这样理解:
-seg blastp的参数,是否对query 序列使用 segmasker来屏蔽低复杂重复,默认 no
-dust blastn的参数,是否对query 序列使用 dustmasker来屏蔽低复杂重复,默认 no
-lcase_masking 对query序列的小写部分进行屏蔽
-soft_masking 是否进行软屏蔽。软屏蔽则是不会使用屏蔽的序列进行种子比对,但是可以延长时候比对。硬屏蔽,则是直接不对屏蔽序列部分进行比对。blastn的默认值是yes,blastp的默认值是no