1. 对目标服务器BIOS进行设置，允许网络唤醒。一般情况下大部分服务器BIOS设置默认是开启了运行网络唤醒。

使用ethtool命令查看服务器网卡设置：

$ sudo ethtool enp4s0f0
Settings for enp4s0f0:
     Supported ports: [ TP ]
     Supported link modes:   10baseT/Half 10baseT/Full 
                             100baseT/Half 100baseT/Full 
                             1000baseT/Half 1000baseT/Full 
     Supported pause frame use: No
     Supports auto-negotiation: Yes
     Supported FEC modes: Not reported
     Advertised link modes:  100baseT/Half 100baseT/Full 
                             1000baseT/Half 1000baseT/Full 
     Advertised pause frame use: No
     Advertised auto-negotiation: Yes
     Advertised FEC modes: Not reported
     Speed: Unknown!
     Duplex: Unknown! (255)
     Port: Twisted Pair
     PHYAD: 1
     Transceiver: internal
     Auto-negotiation: off
     MDI-X: Unknown
     Supports Wake-on: g
     Wake-on: g
     Link detected: no  

保证Wake-on参数的值为g，则表示服务器通过该网卡支持网络唤醒。

2. 确定目标服务器网卡的MAC地址

使用ifconfig查看目标网卡的MAC地址。

$ ifconfig
enp4s0f0: flags=4163  mtu 1500
         inet 192.168.30.6  netmask 255.255.255.0  broadcast 192.168.30.255
         inet6 fe80::574:abf3:4745:5bd6  prefixlen 64  scopeid 0x20

         ether d8:9d:67:77:c0:8c  txqueuelen 1000  (Ethernet)
         RX packets 1604247143  bytes 837254279822 (779.7 GiB)
         RX errors 0  dropped 0  overruns 0  frame 0
         TX packets 906929717  bytes 13220325308560 (12.0 TiB)
         TX errors 0  dropped 0 overruns 0  carrier 0  collisions 0
         device interrupt 43  

以上命令查知目标服务器的目标网卡MAC地址为d8:9d:67:77:c0:8c，broadcast值为192.168.30.255，这是后续唤醒该目标服务器所需要的参数值。

3. 在控制服务器上安装wol软件

$ wget https://sourceforge.net/projects/wake-on-lan/files/wol/0.7.1/wol-0.7.1.tar.gz
$ tar zxf wol-0.7.1.tar.gz
$ cd wol-0.7.1
$ ./configure --prefix=/opt/sysoft/wol-0.7.1/ && make -j 4 && make install
$ cd .. && rm -rf wol-0.7.1
$ echo 'PATH=$PATH:/opt/sysoft/wol-0.7.1/bin/' >> ~/.bashrc
$ source ~/.bashrc

4. 在控制服务器上运行wol命令唤醒目标服务器

控制服务器一般是全天运行的，该服务器和其它用于计算的目标服务器处于同一个内部局域网内。这些局域网内的服务器具有相似的IP地址，前3个IP字段（192.168.30）相同，最后1个字段不同。此外，控制服务器一般能连接外网，具有一个能联外网的网口，该网口具有一个其他IP段的地址。

控制服务器有两个网口在同时工作。一个用于连接外网，一个用于和内网机器连接。使用控制服务器运行wol命令对目标服务器进行唤醒时，需要指定使用内网网口来对内网的目标机器进行唤醒。

$ wol -h 192.168.30.255 d8:9d:67:77:c0:8c
直接对指定MAC地址的局域网内服务器进行唤醒。

-h <string>
    广播到该IP地址或主机，设置该参数的值来选择使用正确的局域网网口来对内网机器进行唤醒。该参数要设置成目标服务器网卡的broadcast值。特别是当控制服务器有多个网口连接时，不设置该参数的值，或设置错误的值，可能会导致不能唤醒目标服务器。

CentOS7的第三方源EPEL不再提供SGE（Son of Grid Engine）软件了，需要自己手动安装SGE。

1. 设置防火墙，开放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

2. 从SGE官网下载最新版本的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

下载SGE软件并进行编译

$ 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
$ ./aimk -no-herd       # 进行图形化安装，不太推荐。我尝试使用该方式编译，导致后续没能部署成功

将编译好的SGE安装到指定的文件夹（需要使用root用户执行）

# mkdir /opt/sysoft/sge
# export SGE_ROOT=/opt/sysoft/sge
# ./scripts/distinst -local -allall -noexit # 虽然普通用户在目标文件夹有写权限，但是程序要对一些文件进行权限修改，使用root用户不会报错。
# cd ../../ && rm sge-8.1.9/ -rf
# echo 'export SGE_ROOT=/opt/sysoft/sge' >> ~/.bashrc
# echo 'PATH=$PATH:/opt/sysoft/sge/bin/:/opt/sysoft/sge/bin/lx-amd64/' >> ~/.bashrc
# source ~/.bashrc    #普通用户也需要进行变量设置

3. 部署SGE前设置主机名

部署SGE前，需要设置好各个节点的主机名，需要修改3个文件。修改配置文件 /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 node2
192.168.30.3 node3
192.168.30.4 node4 

4. 在所有节点上部署SGE

在管理节点上部署SGE：

cd $SGE_ROOT
./install_qmaster

运行部署命令后，会进入交互式界面。基本上全部都按Enter键使用默认设置即可。需要注意的事项是：

1. 有一步骤是启动Grid Engine qmasster服务，可能会启动不了导致失败。原因是多次运行该命令进行部署，第一次会成功运行qmaster daemon，以后重新运行该程序进行部署则会失败。需要删除相应的sge_qmaster进程再进行部署。 
2. 启动Grid Engine qmasster服务，要提供部署SGE的节点主机名信息，按y和Enter键使用一个文件来提供主机信息，输入文件路径/etc/hosts提供主机信息。

只有先进行一个控制节点部署后，才能对各个计算节点进行部署。计算节点的部署比较简单，交互过程全部按Enter即可。

./install_execd

4. 启动SGE软件

部署完毕后，若需要使用SGE软件，则执行如下命令载入SGE的环境变量信息：

$ source /opt/sysoft/sge/default/common/settings.sh

或将该信息添加到~/.bashrc从而永久生效：

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

启动SGE软件方法：

$ /opt/sysoft/sge/default/common/sgemaster    # 控制节点启动
$ /opt/sysoft/sge/default/common/sgeexecd     # 计算节点启动

查看SGE软件运行日志文件：

Qmaster:      /opt/sysoft/sge/default/spool/qmaster/messages
Exec daemon:  /opt/sysoft/sge/default/spool/<hostname>/messages

5. 使用SGE软件

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

使用qconf命令对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
    显示执行主机资源信息

使用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
    删除任务

1. EggNOG数据库简介

eggNOG ( Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups) 公共数据库V5.0版本搜集了5090个生物（477真核生物、4445个代表性细菌和168古菌）和2502个病毒的全基因组蛋白序列。将这些物种分成了379类（taxonomic levels），每类的编号以NCBI的分类编号表示，例如fungi的编号4751。对各个类别物种的全基因组蛋白序列进行同源基因分类，总共得到4.4M个同源基因类（orthologous groups / OGs）。对每个同源基因类进行了多序列比对、系统发育树构建、HMM文件构建和功能注释。

eggNOG数据库是NCBI的COG数据库的扩展，它收集了更全面的物种和更大量的蛋白序列数据。同样进行了同源基因聚类分析和对每个同源基因类的描述和功能分类。eggNOG更强大的功能在于：

1. 对更全面的物种和更大量蛋白序列进行分类。相比于COG数据库纯人工且较为准确的分类，eggNOG数据库扩大物种和序列数据量，采用了非监督聚类方法进行计算。
2. 对每个同源基因类进行了系统发育树构建、HMM模型构建、GO注释、KEGG Pathway注释、SMART/FPAM结构域注释、CAZyme注释等。
3. 提供了本地化软件和网页工具进行eggNOG注释。

eggNOG数据库的对379类物种，每类物种都生成其数据库文件：

0. e5.proteomes.faa文件：收集所有物种的全基因组蛋白序列文件。一个总的数据库文件。

1. members文件：记录OG编号及对应的蛋白序列登录号。编写程序可以从0号文件中提取目标类的蛋白序列数据。
2. annotations文件：记录OG编号的描述信息及所属大类信息。
3. hmms文件：已经构建好的hmm数据库文件。

4. raw_algs文件：多序列比对文件
5. trimmed_algs文件：截短后的多序列比对文件
6. tress文件：统发育树文件。
7. stats文件：统计同源基因类的个数。

2. 使用网页工具进行eggNOG注释

目前（20190414），eggNOG虽然提供了第5版数据信息的下载，但是网页工具仅提供4.5.1版本的分析。网页工具支持DIAMOND和HMMER两种算法进行eggNOG注释。虽然HMMER方法一般准确性更高，但仅支持单次提交不超过5000条序列的FASTA文件。

推荐使用网页工具进行eggNOG注释，简单方便，结果认可度高。以下情况推荐使用
官方提供的软件eggnog-mapper进行本地化注释：

（1）注释的数据量过大，比如宏基因组数据；
（2）使用更准确的HMMER算法一次性注释超过5000条序列；
（3）想更快地进行eggNOG注释；
（4）网络条件较差。

3. 下载并安装eggnog-mapper软件和数据库文件

安装eggnog-mapper软件：

$ wget https://github.com/eggnogdb/eggnog-mapper/archive/1.0.3.tar.gz -O ~/software/eggnog-mapper-1.0.3.tar.gz
$ tar zxf ~/software/eggnog-mapper-1.0.3.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/eggnog-mapper-1.0.3
$ python setup.py install
$ echo 'PATH=$PATH:/opt/biosoft/eggnog-mapper-1.0.3' >> ~/.bashrc
$ source ~/.bashrc
$ download_eggnog_data.py euk bact arch viruses
    下载eggNOG v4.5版本的数据库文件eggnog.db和eggnog_proteins.dmnd，并额外下载指定的HMM数据库。可以根据需要选择指定的数据库即可，下载多个数据库消耗更大的磁盘空间：比如euk消耗76G磁盘空间、bact消耗26G磁盘空间、arch消耗7G磁盘空间、viruses消耗468M磁盘空间。

在eggnog-mapper软件安装目录下的data文件夹下存放有数据库文件：

eggnog_proteins.dmnd    所有蛋白序列的DIMOND数据库，用于DIMOND快速序列比对
eggnog.db    功能注释数据库，用于根据比对结果进行功能注释
hmmdb_levels/    存放有各个物种类别的HMM数据库
OG_fasta/    存放各个物种类别的蛋白序列数据库。

4. 使用eggnog-mapper软件进行eggNOG注释

eggnog-mapper软件运行的原理：（1）程序调用HMMER算法将query序列比对到HMM数据库，得到匹配上的OGs结果。然后将query序列和最优OG中的所有蛋白序列进行比对，得到最优匹配序列（seed ortholog）结果。（2）也可以选择DIAMOND算法将query序列比对到蛋白序列数据库，得到比对结果，选择最优的比对结果作为seed ortholog。（3）基于eggNOG数据库中的对应OG的系统发育树，根据seed ortholog能得到详细的同源基因进化关系，从而剔除亲缘关系较远的直系同源基因，甚至可以剔除一对多的同源基因而仅保留一对一的直系同源基因，将最准确的同源基因的功能注释结果转移给query序列，从而对query序列进行功能注释。

使用HMMER算法进行eggNOG注释：

使用eggnog-mapper软件，推荐使用软件自带的HMMER和DIAMOND程序。我使用最新版本的HMMER软件，出现错误提示：“truct.error: unpack requires a string argument of length ”并得不到结果。
$ export PATH=/opt/biosoft/eggnog-mapper-1.0.3/bin/:$PATH

由于eggNOG数据库较大，推荐将数据库放到内存中，以加快程序运行速度。需要注意的是euk数据比较大，比较费内存，占用了约100G内存左右。
$ emapper.py -d euk -i proteins.fasta -o eggNOG_hmmer --cpu 80 --usemem --no_file_comments

若有多组数据需要进行eggNOG注释，可以分两步运行。先将数据库放到内存中，再分别对多组数据进行运算，这样只读取一次数据库到内存中节约运算时间。
$ emapper.py -d euk --cpu 80 --servermode
待数据库加载完毕后，可以在其它终端中执行数据分析。
$ emapper.py -d euk:localhost:51400 -i proteins1.fasta -o eggNOG1 --no_file_comments
$ emapper.py -d euk:localhost:51400 -i proteins2.fasta -o eggNOG2 --no_file_comments

使用DIAMOND算法 进行eggNOG注释 ：

$ python emapper.py -m diamond -i proteins.fasta -o eggNOG_diamond --cpu 80 --no_file_comments

分两步运行DIAMOND比对和功能注释。自己手动运行DIAMOND，设置一些性能参数和严格阈值能更有效率。
$ diamond blastp --db /opt/biosoft/eggnog-mapper-1.0.3/data/eggnog_proteins --query proteins.fasta --out eggNOG.tab --outfmt 6 --sensitive --max-target-seqs 20 --evalue 1e-5 --id 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1
$ parsing_blast_result.pl -type tab --max-hit-num 1 --db-query proteins.fasta --db-subject /opt/biosoft/eggnog-mapper-1.0.3/eggNOG_Database_5.0/eggnog_proteins.fasta --HSP-num 1 eggNOG.tab | cut -f 1,2,11,12 > eggNOG.emapper.seed_orthologs
$ perl -e '<>; while (<>) { print; }' eggNOG.emapper.seed_orthologs > aa; mv aa eggNOG.emapper.seed_orthologs
$ emapper.py --annotate_hits_table eggNOG.emapper.seed_orthologs -o eggNOG

emapper.py程序常用参数：

输入参数：
-h | --help
    打印帮助信息。
-d | --database <string>
    输入使用的数据库名称。一般选择euk，bact，arch。或输入本地计算机上的HMM数据库路径。
-i <string>
    输入query序列FASTA文件。
--translate
    添加该参数，认为输入的FASTA文件是核酸序列，程序能将核酸序列转换为氨基酸序列后用于比对和注释。
--annotate_hits_table <string>
    输入seed ortholog比对信息。该文件分四列：queryID、subjectID、evalue和bit score。输入该信息，则不需要使用-i参数输入fasta文件了，程序直接运行最后一步功能注释。
--data_dir <string>
    输入eggNOG数据库文件所在的路径。

输出参数：
--output_dir <string>    default: ./
    设置输出文件夹。
-o | --out <string>
    设置输出文件前缀。
--resume
    添加该参数来继续运行程序。
--override
    添加该参数覆盖已有的结果文件。
--no_search
    不进行HMM比对，使用已有的HMM比对结果。默认情况下，程序使用HMMER算法运行分三步：HMM比对、找seed ortholog和功能注释。
--no_refine
    不利用HMM的比对结果进行序列比对寻找seed ortholog，而结束程序，仅仅报告HMM比对结果。
--no_annot
    不进行功能注释，仅报告HMM和seed ortholog的比对结果。
--temp_dir <string>    default: ./
    设置临时文件夹路径。
--no_file_comments
    添加该参数，在输出结果中不输出注释行信息。推荐添加该参数，程序有BUG容易导致seed_ortholog文件中含有#的注释行出现在结果中部，导致程序运行失败。
--keep_mapping_files
    添加该参数，不删除中间文件中的HMMER或DIAMOND比对结果文件。
--report_orthologs
    将用于功能注释的ortholog信息输出到一个单独的文件中。

比对参数：
-m <string>    default: hmmer
    设置比对算法。可以选择的值为hmmer和diamond。
--hmm_maxhits <int>    default: 1
    设置HMMER算法得到的hit结果数量。
--hmm_evalue <float>    default: 0.001
    设置HMMER算法的evalue阈值。
--hmm_score <int>    default: 20
    设置HMMER算法的score阈值。
--hmm_maxseqlen <int>    default: 5000
    忽略长度超过此阈值的query序列。
--hmm_qcov <float>    default: disabled
    设置HMMER算法中query序列的最小覆盖率阈值。
--dmnd_db <string>
    当使用DIAMOND算法时，设置DIAMOND数据库路径。
--matrix <string>    default: BLOSUM62
    设置DIAMOND算法的计分举证。
--gapopen <int>
    设置DIAMOND算法中打开Gap罚分值。
--gapextend <int>
    设置DIAMOND算法中延长Gap罚分值。
--seed_ortholog_evalue <float>    default: 0.001
    进行seed ortholog鉴定时，设定的evalue阈值。
--seed_ortholog_score <int>    default: 60
    进行seed ortholog鉴定时，设定的score阈值。 

功能注释参数：
--tax_scope <string>
    仅使用指定物种类下的同源基因进行功能注释。默认设置下是分别对每条query序列进行自动调节。
--target_orthologs <string>
    进行功能注释时，设置选择什么类型的同源基因进行功能注释。可以选择的值有：one2one,many2one,one2many,many2many,all。
--go_evidence <string>
    设置使用什么类型的GO terms进行功能注释。可以设置的值有：experimental，仅使用源自实验证据支持的GO terms；non-electronic，仅使用非电子计算归纳的Go terms。

性能参数：
--cpu <int>
    设置程序使用的CPU线程数。
--usemem
    添加该参数，能将HMM数据库读取到内存中，增加计算性能。在程序运行结束后，释放内存。
--servermode
    添加该参数，会自动添加--usemem参数，将HMM数据库读取到内存中，并一直将HMM数据库维持到内存中，并提供数据库服务，直到强行结束程序。

5. eggNOG结果解读和分类图绘制

1. PHI-base数据库简介

PHI-base数据库从文献中收集经过实验验证了的致病基因和效应基因的序列。目前（20190411）数据库版本为4.6版本，从3011篇文献中收集了263种致病菌（细菌、真菌、原生动物和线虫）的6438个基因和194种宿主（植物占~70%、脊椎动物、昆虫、线虫和真菌）的11340种相互关系，其中包含510中疾病。PHI-base将收集到的参考文献信息、 基因信息、病原和宿主信息、疾病信息、表型和相互关系等记录到数据库中，并提供关键词进行搜索。同时，PHI-base还提供BLAST网页工具，用于将序列比对到数据库中，从而根据BLAST结果进行致病基因分析。此外，PHI-base也提供了数据库下载，能更方便地进行本地化批量注释。

PHI-base将数据库中基因的表型分为多类：

1. Loss of pathogenicity: 导致病原菌致病能力丧失的基因。
2. Reduced virulence： 导致病原菌致病能力减弱的基因。
3. Unaffected pathogenicity： 对病原菌致病性没有影响的基因。
4. Increased pathogenicity (Hypervirulence)： 导致病原致病能力增强的基因。
5. Effector (plant avirulence determinant)：致病效应基因。效应基因能直接或间接地让宿主识别病原，并激活植物防御反应机制，导致致病失败。此外，少数效应基因是一些易感宿主致病所必须的基因，大部分易感宿主致病不需要效应基因。
6. Enhanced antagonism: 导致植物内生菌表型出症状，病原菌在和宿主互作中占优势。
7. Lethal：致死因子。缺少致死基因或减少其表达能导致病原菌无法存活。这类基因的产物是生命所必不可少的。注意该类基因不是直接的致病基因，仅是保证病原生存的基因，和增强病原菌的致病能力没有关系。 
8. chemistry_target resistance to chemical: 抗药基因。
9. chemistry_target sensitivity to chemical: 感药基因。

对全基因组基因型PHI注释，大部分基因属于1~3类。属于第4类Increased pathogenicity (Hypervirulence)类型的基因才是关键致病基因；属于5类Effector类型和第6类Enhanced antagonism基因和致病也比较相关。

2. 下载并部署PHI-base数据库

貌似目前（201904）该数据库不能提供下载。我保留了之前4.5版本的文件并进行了BLAST数据创建：

$ mkdir /opt/biosoft/PHI-base4.5
$ cd /opt/biosoft/PHI-base4.5
# 数据库文件phi45.fas和phibase45.csv
$ grep -v -P "^#" phi45.fas > phi45.fasta
    .fas文件首行含有#符号，需要删除，才能构建BLAST数据库。
$ /opt/biosoft/ncbi-blast-2.2.28+/bin/makeblastdb -in phi45.fasta -dbtype prot -title protein -parse_seqids -out phi -logfile phi.protein.log
    .fasta文件头部的序列名超过了50个字符长度，新版本的makeblastdb不支持。
$ cp phi* /opt/biosoft/wwwblast/db/

3. 使用BLASTP进行PHI注释，并编写程序提炼结果

对全基因组蛋白序列和PHI数据库进行比对后，一条序列可能有多个Hits结果，而每个Hits结果可能有多种表型。为了得到更准确的PHI注释结果，推荐按如下方法进一步分析：首先对BLAST结果设定严格的evalue、identity和coverage阈值进行过滤；然后对每个基因BLAST结果中各种PHI表型的次数进行排序；最后选次数最高、BLAST比对得分最大的表型作为基因的注释结果，可以得到一个基因对应一个表型的结果；或者在前者基础上，额外增加次数出现也较多（出现的次数不低于最优表型次数的60%）的表型，可以得到一个基因对应多个表型的结果。

$ blast.pl --CPU 40 /opt/biosoft/wwwblast/db/phi proteins.fasta > phi.xml
$ parsing_blast_result.pl --evalue 1e-5 --identity 0.3 --subject-coverage 0.3 --query-coverage 0.3 --HSP-num 1 phi.xml > phi.tab
$ perl -e '<>; while (<>) { chomp; my @field = split /\t/, $_; $field[1] =~ m/.*#([^\t]+)/; my $annot = $1; my @annot = split /__/, $annot; foreach (@annot) { $annot{$field[0]}{$_} ++; $annot_score{$field[0]}{$_} = $_[11] if $annot_score{$field[0]}{$_} < $_[11];} } foreach my $gene (sort keys %annot) { my @annot = sort { $annot{$gene}{$b} <=> $annot{$gene}{$a} or $annot_score{$gene}{$b} <=> $annot_score{$gene}{$a} } keys %{$annot{$gene}}; print "$gene\t$annot[0]\n"; $stats{$annot[0]}{$gene} = 1; } foreach (sort keys %stats) { my @gene = sort keys %{$stats{$_}}; my $num = @gene; my $gene = join ",", @gene; print STDERR "$_\t$num\t$gene\n"; }' phi.tab > phi.annot 2> phi.stats
$ perl -e 'print "geneID\tPhenotype\tReportNum\n"; <>; while (<>) { chomp; my @field = split /\t/, $_; $field[1] =~ m/.*#([^\t]+)/; my $annot = $1; my @annot = split /__/, $annot; foreach (@annot) { $annot{$field[0]}{$_} ++; $annot_score{$field[0]}{$_} = $_[11] if $annot_score{$field[0]}{$_} < $_[11];} } foreach my $gene (sort keys %annot) { my @annot = sort { $annot{$gene}{$b} <=> $annot{$gene}{$a} or $annot_score{$gene}{$b} <=> $annot_score{$gene}{$a} } keys %{$annot{$gene}}; my $max_num = $annot{$gene}{$annot[0]}; print "$gene\t$annot[0]\t$annot{$gene}{$annot[0]}\n"; $stats{$annot[0]}{$gene} = 1; shift @annot; foreach (@annot) { if (($annot{$gene}{$_} / $max_num) >= 0.6) { print "$gene\t$_\t$annot{$gene}{$_}\n"; $stats{$_}{$gene} = 1;} } } foreach (sort keys %stats) { my @gene = sort keys %{$stats{$_}}; my $num = @gene; my $gene = join ",", @gene; print STDERR "$_\t$num\t$gene\n"; }' phi.tab > phi.annot 2> phi.stats 

NCBI已经不再维护并下架了ePCR软件，转而推荐使用其Primer-BLAST网页工具。这对于单个引物设计任务比较方便，但不利于基因组较大的非模式物种的特异性引物设计或批量化的引物设计。

本教程讲解使用NCBI-ePCR和Primer3进行引物批量化设计。

1.下载并安装NCBI-ePCR和Primer3软件

$ wget http://ftp.debian.org/debian/pool/main/e/epcr/epcr_2.3.12-1.orig.tar.gz -P ~/software
$ tar zxf ~/software/epcr_2.3.12-1.orig.tar.gz -C /opt/biosoft/
$ make -j 4
$ echo 'PATH=$PATH:/opt/biosoft/e-PCR-2.3.12/' >> ~/.bashrc
$ source ~/.bashrc

$ wget https://sourceforge.net/projects/primer3/files/primer3/2.4.0/primer3-2.4.0.tar.gz -P ~/software/
$ tar zxf ~/software/primer3-2.4.0.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/primer3-2.4.0/src/
$ make all
$ echo 'PATH=$PATH:/opt/biosoft/primer3-2.4.0/src/' >> ~/.bashrc
$ source ~/.bashrc

2. 使用ePCR进行引物验证

首先，使用famap命令和fahash命令分两步将基因组序列转换为哈希数据库。

$ famap -t N -b genome.famap genome.fasta
     程序将FASTA格式的序列转换为famap数据库文件。

常用参数：
-t <string>
    设置碱基类型。可以设置4种值：n，允许含有小写碱基atcgn，其它字符转换为n或N；nx，允许含有小写碱基和兼并碱基字符，其它字符转换为n或N；N，仅允许大写碱基，atcgn自动转换为大写，其它字符转换为N；NX，允许大写碱基和兼并碱基字符，其它字符转换为N。 -b <string>
    设置输出的famap数据库文件路径。

$ fahash -b genome.hash -w 12 -f 3 genome.famap
程序进一步将famap文件转换为hash数据库文件。

常用参数：
-b <string>
    输出hash数据库文件。
-w <int>
    设置wordsize长度。
-f <int>
    设置wordcnt长度。

然后，使用re-PCR将引物和数据库进行比对，寻找引物比对结果。

分别对多个引物进行比对，得到各个引物的匹配结果：
$ re-PCR -p genome.hash -n 1 -g 1 ACTATTGATGATGA AGGTAGATGTTTTT …

输入一对引物，并设置产物长度期望值，得到一对引物的匹配结果：
$ re-PCR -s genome.hash -n 1 -g 1 ACTATTGATGATGA AGGTAGATGTTTTT 50-1000
常用参数：
-p 
     输入hash数据库，在命令行中直接输入引物序列，将引物和数据库进行比对。
-s 
     输入hash数据库，在命令行中必须输入一对引物和产物长度期望范围，得到一对引物的匹配结果。
-n 
     设置允许的错配碱基数。
-g 
     设置允许的gap数。

3. Primer3软件的使用

Primer3是命令行形式的引物设计软件，能很好地用于引物的批量设计。 其主程序是primer3_core。 该命令的输入是Boulder-IO格式，适合于软件读入数据；输出文件默认下是适合人类阅读的格式，也可以是Boulder-IO格式，有助于下一步引物结果的批量操作。 Boulder-IO格式以文本形式记录着引物设计信息，且每个引物信息的结尾使用“等于换行符”分隔。每个记录由多个标签和对应的值构成，用于指定输入信息或输出结果。例如：

SEQUENCE_ID=example
SEQUENCE_TEMPLATE=GTAGTCAGTAGACNATGACNACTGACGATGCAGACNACACACACACACACAGCACACAGGTATTAGTGGGCCATTCGATCCCGACCCAAATCGATAGCTACGATGACG
SEQUENCE_TARGET=37,21
PRIMER_TASK=generic
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=1
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_OPT_SIZE=18
PRIMER_MIN_SIZE=15
PRIMER_MAX_SIZE=21
PRIMER_MAX_NS_ACCEPTED=1
PRIMER_PRODUCT_SIZE_RANGE=75-100
P3_FILE_FLAG=1
SEQUENCE_INTERNAL_EXCLUDED_REGION=37,21
PRIMER_EXPLAIN_FLAG=1

Primer3使用示例和参数：

$ primer3_core -p3_settings_file p3_settings_file -strict_tags -format_output input_file result.p3.out
     程序的输入文件要求是Boulder-IO格式，如果没有输入文件，则会从标准输入读取数据。 

常用参数：
-p3_settings_file
    用于输入primer3的配置文件。这是因为primer3可设置的参数太多了，将大部分参数放入到该配置文件，有利于参数的输入。primer3的默认参数不好，特别是默认参数下没有开启热力学计算。推荐使用配置文件来根据自己的需求来设定相关参数。此外，输入文件中的参数设置能取代此文件中的设定值。
    此文件格式：第1行固定为"Primer3 File - http://primer3.sourceforge.net"；第2行固定为"P3_FILE_TYPE=settings"；第3行是个空行；从第4行开始则是标准的Boulder-IO格式内容。
-format_output
    让primer3_core产生人类易读的结果，否则产生机器易读的结果（Boulder-IO格式结果）。
-strict_tags
    要求输入文件中的标签要全部正确。设置该参数后，如果有标签不能被程序识别，则报错并停止执行。不设置此参数，则软件忽略不识别的参数。
-p3_settings_file=file_path
    指定 primer_core 的配置文件，该配置文件的设定会取代默认设置。当然，
-echo_settings_file
    打印出p3_settings_file中的设置信息。如果没有指定设置文件，或含有-format_output，则该参数失效。
-output=file_path
    指定输出文件路径，如果不指定，则输出到标准输出。
-error=file_path
   指定错误信息输出路径，如果不指定，则输出到stderr中。

Primer3使用的难点是根据自身需求准备p3_settings_file文件。我的一个示例如下（使用时需要去掉#注释部分和空行，且必须第三行留空行）：

Primer3 File - http://primer3.sourceforge.net
P3_FILE_TYPE=settings

P3_FILE_ID=P3 Settings from Lianfu Chen
PRIMER_FIRST_BASE_INDEX=1
PRIMER_TASK=generic
PRIMER_NUM_RETURN=5
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=0
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_PICK_ANYWAY=1
PRIMER_THERMODYNAMIC_PARAMETERS_PATH=/opt/biosoft/primer3-2.4.0/src/primer3_config/

### 引物 TM 值设定: ###
PRIMER_TM_FORMULA=1                 # TM 的计算方法， 1 表示使用 the SantaLucia parameters (Proc Natl Acad Sci 95:1460-65)
PRIMER_MIN_TM=55.0  
PRIMER_OPT_TM=60.0
PRIMER_MAX_TM=65.0
PRIMER_PAIR_MAX_DIFF_TM=5.0         # 两个引物之间的 TM 值最多相差 5 摄氏度
PRIMER_WT_TM_LT=0
PRIMER_WT_TM_GT=0
PRIMER_PAIR_WT_DIFF_TM=0.0

### 引物长度设定: ###
PRIMER_MIN_SIZE=18
PRIMER_OPT_SIZE=20
PRIMER_MAX_SIZE=22
PRIMER_WT_SIZE_LT=0
PRIMER_WT_SIZE_GT=0

### 引物 GC 含量设定：###
PRIMER_MIN_GC=30.0
PRIMER_MAX_GC=70.0
PRIMER_WT_GC_PERCENT_LT=0.0
PRIMER_WT_GC_PERCENT_GT=0.0

### 引物的热力学计算: ###
# 开启热力学计算，开启后，根据热力学的值 TH 值来计算罚分。 TH 的罚分方法为：罚分系数 * (1 / (引物TM - 4 - TH值))。
# 该方法好处是，TH 值越大，罚分力度越重。
PRIMER_THERMODYNAMIC_OLIGO_ALIGNMENT=1

# 引物自身进行反向互补，罚分系数计算为 9 / ( 1/(60-4-45) - 1/(60-4-0) ) = 123.2
# 这样计算罚分系数的原则是，若 TM 为最佳 60 摄氏度的时候，所允许的最大罚分值减去最小罚分值为 9，和 3' 端碱基的稳定性的罚分额度一致。
PRIMER_MAX_SELF_ANY=8.00
PRIMER_WT_SELF_ANY=0.0
PRIMER_MAX_SELF_ANY_TH=45.00
PRIMER_WT_SELF_ANY_TH=123.2

# 引物自身进行 3' 端反向互补形成引物二聚体，罚分系数计算为 9 / ( 1/(60-4-35) - 1/(60-4-0) ) = 302.4
PRIMER_MAX_SELF_END=3.00
PRIMER_WT_SELF_END=0.0
PRIMER_MAX_SELF_END_TH=35.00
PRIMER_WT_SELF_END_TH=302.4

# left primer 和 right primer 序列的反向互补
PRIMER_PAIR_MAX_COMPL_ANY=8.00
PRIMER_PAIR_WT_COMPL_ANY=0.0
PRIMER_PAIR_MAX_COMPL_ANY_TH=45.00
PRIMER_PAIR_WT_COMPL_ANY_TH=123.2

# left primer 和 right primer 进行 3' 端反向互补形成引物二聚体
PRIMER_PAIR_MAX_COMPL_END=3.00
PRIMER_PAIR_WT_COMPL_END=0.0
PRIMER_PAIR_MAX_COMPL_END_TH=35.00
PRIMER_PAIR_WT_COMPL_END_TH=302.4

# 发夹结构，罚分系数计算为 9 / ( 1/(60-4-24) - 1/(60-4-0) ) = 672
PRIMER_MAX_HAIRPIN_TH=24.00
PRIMER_WT_HAIRPIN_TH=672

# 3' 端碱基的稳定性
PRIMER_MAX_END_STABILITY=9.0
PRIMER_WT_END_STABILITY=1

### 碱基序列设定： ###
PRIMER_LOWERCASE_MASKING=0              # 模板序列中包含小写字符不影响引物设计
PRIMER_MAX_POLY_X=4                     # 引物序列中不能包含单核苷酸连续长度超过 4 bp
PRIMER_MAX_NS_ACCEPTED=0                # 引物中允许的 N 的数目
PRIMER_WT_NUM_NS=0.0                    # 每个 N 的罚分
PRIMER_MAX_END_GC=5                     # 引物中 3' 端 5bp 碱基中允许的最大 Gs 或 Cs 的数目
PRIMER_GC_CLAMP=0                       # 引物中 3' 端碱基中不能出现连续的 Gs 和 Cs 序列
PRIMER_LIBERAL_BASE=1                   # 是否允许有诸如 N A R Y 等类型的碱基。必须在设定PRIMER_MAX_NS_ACCEPTED 不为 0 后方有效。
PRIMER_LIB_AMBIGUITY_CODES_CONSENSUS=0  # 如果设置为 1，则 C 能与 S 完美匹配，任意碱基能和 N 完美匹配。

### 碱基质量设定： ###
PRIMER_MIN_QUALITY=0                    # primer 序列允许最小的碱基质量
PRIMER_MIN_END_QUALITY=0
PRIMER_QUALITY_RANGE_MIN=0
PRIMER_QUALITY_RANGE_MAX=100
PRIMER_WT_SEQ_QUAL=0.0
PRIMER_WT_END_QUAL=0.0

## 引物位置设置： ###
PRIMER_SEQUENCING_LEAD=50               # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 表明引物的 3' 端距目标区域有 50bp。
PRIMER_SEQUENCING_SPACING=500           # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 该值决定了在同一条链上的两个引物的距离.
PRIMER_SEQUENCING_INTERVAL=250          # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 该值决定了在不同链上的两个引物的距离。
PRIMER_SEQUENCING_ACCURACY=20           # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 该值决定了引物设计的区间.
PRIMER_OUTSIDE_PENALTY=0
PRIMER_INSIDE_PENALTY=-1.0
PRIMER_WT_POS_PENALTY=0.0

### PCR 反应体系设定: ###
PRIMER_SALT_MONOVALENT=50.0                 # 单价盐离子浓度(mM)
PRIMER_SALT_CORRECTIONS=1                   # PRIMER_SALT_CORRECTIONS=1 means use the salt correction in SantaLucia et al 1998
PRIMER_SALT_DIVALENT=1.5                    # 二价镁离子浓度(mM)
PRIMER_DNTP_CONC=0.6                        # 总dNTPs浓度(mM)
PRIMER_DNA_CONC=50.0                        # DNA产物浓度(mM)

### PCR 产物的设定： ###
PRIMER_PRODUCT_MIN_TM=-1000000.0
PRIMER_PRODUCT_OPT_TM=0.0
PRIMER_PRODUCT_MAX_TM=1000000.0
PRIMER_PAIR_WT_PRODUCT_TM_LT=0.0
PRIMER_PAIR_WT_PRODUCT_TM_GT=0.0
PRIMER_PRODUCT_OPT_SIZE=0
PRIMER_PAIR_WT_PRODUCT_SIZE_LT=0.0
PRIMER_PAIR_WT_PRODUCT_SIZE_GT=0.0

## 罚分因子： ###
PRIMER_PAIR_WT_PR_PENALTY=1.0                # left primer和right primer的罚分之和。此和乘以此系数，再加其它罚分作为最终罚分。
PRIMER_PAIR_WT_IO_PENALTY=0.0                # 将 internal oligo 的罚分乘以此系数，加入到引物的最终罚分中。

### internal oligo 的设定：###
# TM 值
PRIMER_INTERNAL_MIN_TM=57.0
PRIMER_INTERNAL_OPT_TM=60.0
PRIMER_INTERNAL_MAX_TM=63.0
PRIMER_INTERNAL_WT_TM_LT=1.0
PRIMER_INTERNAL_WT_TM_GT=1.0

# 长度
PRIMER_INTERNAL_MIN_SIZE=18
PRIMER_INTERNAL_OPT_SIZE=20
PRIMER_INTERNAL_MAX_SIZE=27
PRIMER_INTERNAL_WT_SIZE_LT=1.0
PRIMER_INTERNAL_WT_SIZE_GT=1.0

# GC 含量
PRIMER_INTERNAL_MIN_GC=20.0
PRIMER_INTERNAL_MAX_GC=80.0
PRIMER_INTERNAL_OPT_GC_PERCENT=50.0
PRIMER_INTERNAL_WT_GC_PERCENT_LT=0.0
PRIMER_INTERNAL_WT_GC_PERCENT_GT=0.0

# 热力学
PRIMER_INTERNAL_MAX_SELF_ANY=12.00
PRIMER_INTERNAL_WT_SELF_ANY=0.0
PRIMER_INTERNAL_MAX_SELF_END=12.00
PRIMER_INTERNAL_WT_SELF_END=0.0

# 碱基序列
PRIMER_INTERNAL_MAX_POLY_X=5
PRIMER_INTERNAL_MAX_NS_ACCEPTED=0

# 碱基质量
PRIMER_INTERNAL_MIN_QUALITY=0
PRIMER_INTERNAL_WT_END_QUAL=0.0
PRIMER_INTERNAL_WT_SEQ_QUAL=0.0

# 非引物数据库
#PRIMER_INTERNAL_MAX_LIBRARY_MISHYB=12.00
#PRIMER_INTERNAL_WT_LIBRARY_MISHYB=0.0

# PCR 反应体系
PRIMER_INTERNAL_SALT_MONOVALENT=50.0
PRIMER_INTERNAL_SALT_DIVALENT=1.5
PRIMER_INTERNAL_DNA_CONC=50.0
PRIMER_INTERNAL_DNTP_CONC=0.0
=

4. 编写程序调用primer3进行引物设计再调用e-PCR软件进行特异性验证

编写程序primer3_with_ePCR_validation.pl，输入模板序列的FASTA文件，即可对所有的模板序列批量化进行引物设计。若同时输入全基因组的序列，进一步可以对设计出的引物进行特异性筛选。

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

my $usage = <<USAGE;
Usage:
    $0 [options] seq.fasta > primer3.out

    程序对seq.fasta文件中的序列调用primer3进行引物批量化设计。并将结果输出到标准输出。

    --min_product_length <INT>  default: 80
    设置引物得到的最小产物长度。

    --max_product_length <INT>  default: 150
    设置引物得到的最大产物长度。

    --primer-num <int>    default: 5
    设置使用primer3软件设计的引物对数量。

    --p3_setting_file <STRING>
    程序使用Primer3进行引物批量设计，该参数设置所使用的Primer3配置文件路径。该参数是必须的，以利于primer3设置正确有效的引物。

    --CPU <int>    default: 1
    设置并行化运行primer3的任务数。程序调用ParaFly命令进行并行化运算，需要能直接在terminal中直接运行ParaFly命令。

    --db <string>
    若使用该参数输入FASTA格式的数据库序列，程序能额外对primer3得到的引物进行e-PCR验证。

    --max-mismatch <int>    default: 1
    --max-gap <int>    default: 1
    设置e-PCR过程中允许引物和数据库序列最大的错配和gap数量。设置更高能得到更好的特异性。

    --max-hit <int>    default: 1
    设置允许e-PCR对数据库搜索得到的最大结果数量。若模板序列存在数据库中，推荐设置为1，若模板序列不存在数据库中，推荐设置为0。

    --no-primer-seq <string>
    若设置该参数的值，则生成相应的文件，输出没有引物结果的序列。

    --tmp-dir <string>    default: primer3_with_ePCR_validation.tmp
    设置临时文件夹，程序的中间文件放入该文件夹中。

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

...

1. 动物转录因子数据AnimalTFDB3.0简介

动物转录因子数据库AnimalTFDB3.0对97个动物基因组的转录因子（Transcription Factor）和转录辅助因子（Transcription cofactor）进行了归纳整理。基于DNA结合结构域，将动物转录因子分成了73个基因家族，将转录辅助因子分成了83个基因家族。此外，动物转录因子分为六大类（Basic Domain Group、Zinc-Coordinating Group、Beta-Scaffold Factors、Helix-turn-helix、Other Alpha-Helix Group和Unclassified Structure），动物转录辅助因子也分为六大类（Co-activator/repressors、Chromatin Remodeling Factors、General Cofactors、Histone-modifying Enzymes、Cell Cycle和Other Cofactors）。

动物转录因子数据库AnimalTFDB3.0提供了网页工具进行转录因子分析。该网页工具一次仅支持上传不超过1000条蛋白序列。该网页工具采用HMM算法进行计算，分析速度比较快。数据库作者对大部分转录因子家族采用了PFam数据库的HMM模型，并自己构建了一些转录因子家族的HMM模型。但作者没有提供HMM数据库的下载，仅提供了97个动物转录因子和转录辅助因子蛋白序列的FASTA格式文件下载。

使用AnimalTFDB3.0数据库的网页工具仅能进行转录因子分析，且由于仅能支持不超过1000条序列进行分析而比较麻烦，不利于动物全基因组的转录因子分析。以下讲解下载AnimalTFDB3.0数据库FASTA文件，并自行编写程序进行转录因子和转录辅助因子注释。

2. 下载AnimalTFDB3.0数据库蛋白序列

首先，进入AnimalTFDB3.0数据库下载页面，采用复制粘贴方式获取97个物种的物种名，保存到文件species.txt中。

Ailuropoda melanoleuca;Anas platyrhynchos;Anolis carolinensis;Aotus nancymaae;Astyanax mexicanus;Bos taurus;Caenorhabditis elegans;Callithrix jacchus;Canis familiaris;Capra hircus;Carlito syrichta;Cavia aperea;Cavia porcellus;Cebus capucinus;Cercocebus atys;Chinchilla lanigera;Chlorocebus sabaeus;Choloepus hoffmanni;Ciona intestinalis;Ciona savignyi;Colobus angolensis palliatus;Cricetulus griseus chok1gshd;Cricetulus griseus crigri;Danio rerio;Dasypus novemcinctus;Dipodomys ordii;Drosophila melanogaster;Echinops telfairi;Equus caballus;Erinaceus europaeus;Felis catus;Ficedula albicollis;Fukomys damarensis;Gadus morhua;Gallus gallus;Gasterosteus aculeatus;Gorilla gorilla;Heterocephalus glaber female;Heterocephalus glaber male;Homo sapiens;Ictidomys tridecemlineatus;Jaculus jaculus;Latimeria chalumnae;Lepisosteus oculatus;Loxodonta africana;Macaca fascicularis;Macaca mulatta;Macaca nemestrina;Mandrillus leucophaeus;Meleagris gallopavo;Mesocricetus auratus;Microcebus murinus;Microtus ochrogaster;Monodelphis domestica;Mus caroli;Mus musculus;Mus pahari;Mus spretus;Mustela putorius furo;Myotis lucifugus;Nannospalax galili;Nomascus leucogenys;Notamacropus eugenii;Ochotona princeps;Octodon degus;Oreochromis niloticus;Ornithorhynchus anatinus;Oryctolagus cuniculus;Oryzias latipes;Otolemur garnettii;Ovis aries;Pan paniscus;Pan troglodytes;Papio anubis;Pelodiscus sinensis;Peromyscus maniculatus bairdii;Petromyzon marinus;Poecilia formosa;Pongo abelii;Procavia capensis;Propithecus coquereli;Pteropus vampyrus;Rattus norvegicus;Rhinopithecus bieti;Rhinopithecus roxellana;Saimiri boliviensis boliviensis;Sarcophilus harrisii;Sorex araneus;Sus scrofa;Taeniopygia guttata;Takifugu rubripes;Tetraodon nigroviridis;Tupaia belangeri;Tursiops truncatus;Vicugna pacos;Xenopus tropicalis;Xiphophorus maculatus

然后，根据物种名，编写下载数据库文件的脚本（每个物种下载4个数据文件）并批量下载数据。

$ mkdir /opt/biosoft/AnimalTFDB3.0
$ cd /opt/biosoft/AnimalTFDB3.0

生成文件species.list，包含97个物种名，每行一个物种名。根据http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/download网>页中的信息，使用vim编辑和perl单行命令抓取全部的物种名信息。
$ perl -p -i -e 's/;/\n/g;' species.txt

$ mkdir data_of_97_species
$ cd data_of_97_species
$ perl -e 'while (<>) { chomp; s/\s+/_/g; print "wget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF_protein.fa\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_Cof_protein.fa\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF_cofactors\n"; }' ../species.list > download.list
$ ParaFly -c download.list -CPU 2

Caenorhabditis_elegans的FASTA文件有问题，存在序列名一致的多条序列，推荐进行修改。
$ perl -p -i -e 's/>(\S+)\t(\S+)/>$2\t$1/' data_of_97_species/Caenorhabditis_elegans_*_protein.fa
$ cd ..

分别合并转录因子数据库和转录辅助因子数据库的FASTA文件。

得到转录因子序列数据库，FASTA文件头部仅包含序列名和基因家族信息。
$ cat data_of_97_species/*TF_protein.fa > AnimalTFDB3.0_TF_protein.fasta
$ perl -p -i -e 's/^(>\S+)\t.*?\t.*?\t/$1\t/' AnimalTFDB3.0_TF_protein.fasta


得到转录辅助因子数据库，FASTA文件头部仅包含序列名和基因家族信息。
$ cat data_of_97_species/*_Cof_protein.fa > AnimalTFDB3.0_Cof_protein.fasta
由于Cof_protein.fa文件中不包含转录辅助因子基因家族信息，需要编写程序读取TF_cofactors文件信息并将基因家族信息写入FASTA文件
$ perl -e 'while (<>) { @_ = split /\t/; $gf{$_[2]} = $_[3]; } open IN, "AnimalTFDB3.0_Cof_protein.fasta"; while (<IN>) { if (m/^>(\S+)\t(\S+)/) { if (exists $gf{$1}) {print ">$1\t$gf{$1}\n";} elsif (exists $gf{$2}) {print ">$2\t$gf{$2}\n";} else {print STDERR "No type: $_\n";} } else { print; } }' data_of_97_species/*_TF_cofactors > aa; mv aa AnimalTFDB3.0_Cof_protein.fasta

3. 使用DIAMOND将数据库所有蛋白和数据库自身进行比对，计算每个基因家族的BLAST阈值

使用ALL_VS_ALL的方式分别对转录因子数据库和转录辅助因子数据库进行DIAMOND比对（evalue 1e-5、identity 30%且coverage 30%）。根据DIAMOND比对结果，得到每个转录因子家族中每个基因的比对信息（evalue、hitNum），编写程序进而计算每个基因家族的BLAST阈值。

$ diamond --in AnimalTFDB3.0_TF_protein.fasta --db AnimalTFDB3.0_TF_protein
$ diamond blastp --db AnimalTFDB3.0_TF_protein --query AnimalTFDB3.0_TF_protein.fasta --out TF.tab --outfmt 6 --sensitive --max-target-seqs 50000 --evalue 1e-5 --id 30 --subject-cover 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1

编写程序geneFamily_BLAST_threshold_of_sensitivity.pl，输入DIAMOND结果文件和数据库FASTA文件，得到各个sensitivity级别上每个基因家族的evalue和hitNum阈值。

4. 将全基因组蛋白序列和数据库进行比对，再根据阈值文件进行准确注释。

1. FALCON软件简介

FALCON是PacBio公司开发的一款用于三代基因组De novo组装软件。相比于HGAP4软件，FALCON软件的基因组组装原理基本一致。但FALCON使用命令行运行，更适合于大基因组的组装，且能分析双倍体序列，并在基因组组装结果中给出包含变异位点信息的等位基因序列（alternative contigs / a-contigs）和主要的基因组序列（primary contig / p-contig）。每一条a-contig都有其对应的p-contig序列。因此，FALCON软件适合双倍体物种的基因组组装，能给出单倍的基因序列。其基因组组装结果中的p-contigs序列总长度要小于其它基因组组装软件（例如Canu和HGAP）的基因组序列。

FALCON-Unzip则是真正的单倍型组装软件，它能在FALCON或HGAP4软件的基因组组装结果基础上，利用较长的PacBio reads进行单倍型分析，对p-contigs序列向单倍型进行转换，同时输出单倍型序列（haplotig）区块。

2. FALCON软件下载与安装

FALCON软件在2018年在算法上有了较大的更改，将FALCON在内的多个软件整合到了pb-assembly软件（https://github.com/PacificBiosciences/pb-assembly）中。同时提供快捷的bioconda安装方法：

$ wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ sh Miniconda3-latest-Linux-x86_64.sh
    进入交互式界面：首先，输入yes，按Enter键，表示同意软件的license；然后，输入conda程序安装路径/opt/biosoft/miniconda3_for_pb-assembly，按enter键； 再设置是否添加Minicoda3到PATH变量，直接按Enter键，表示选择no，以免和系统自带的软件冲突。若需要切换到miniconda3环境，输入命令export PATH=/opt/biosoft/miniconda3_for_pb-assembly/bin:$PATH即可。
$ export PATH=/opt/biosoft/miniconda3_for_pb-assembly/bin:$PATH
$ conda config --add channels defaults
$ conda config --add channels bioconda
$ conda config --add channels conda-forge
    配置bioconda channel（https://bioconda.github.io/index.html）
$ conda install pb-assembly
    安装pb-assembly软件
$ source activate /opt/biosoft/miniconda3_for_pb-assembly
    激活FALCON软件的运行环境

3. FALCON软件的使用

FALCON软件运行示例：

首先，载入FALCON软件的运行环境。

然后，准备三种输入文件：测序数据FASTA格式文件、包含测序数据文件路径的FOFN列表文件、程序运行的参数配置文件。下载基因组大小为4.6Mb的E. coli基因组PacBio测序数据。
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139\
_c100713652400000001823152404301535_s1_p0.subreads.tar | tar xvf -
$ ls */ecoli.?.fasta > input.fofn
$ wget https://pb-falcon.readthedocs.io/en/latest/_downloads/fc_run_ecoli_local.cfg
当使用最新版的FALCON进行组装时，需要对此配置文件进行一些修改。

最后，对E. coli基因组进行组装
$ fc_run.py fc_run_ecoli_local.cfg

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

运行FALCON最关键的步骤是准备配置文件fc_run.cfg。一个配置文件示例及其参数讲解如下：

[General]
输入参数：
input_fofn = input.fofn
设置数据输入文件为input.fofn，该文本文件中每行是一个PacBio测序数据的FASTA文件路径。可以使用相对路径。
input_type = raw
设置输入数据类型为原始测序数据（raw），或者为修正后的数据（preads）。

对数据进行分块：
pa_DBsplit_option = -x500 –s200
设置预组装（pre-assembly）过程中对raw reads数据进行分块的参数。这些参数会传递给Dbsplit命令来执行数据分块。
-s<int>    default: 200
表示每份数据含有50Mb的数据量，该参数默认值为200。默认参数情况下，使用daligner进行比对需要消耗约16G内存。推荐设置该参数的值为预测的基因组大小。此外，软件作者推荐对小于10Mb的基因组设置该参数值为50，对于大基因组设置该参数值为200。
-x<int>
忽略长度小于指定阈值的reads。
ovlp_DBsplit_option=-x500 –s200
设置OLC组装过程中对preads数据进行分块的参数。

对重复序列进行屏蔽：
pa_HPCTANmask_option = -k18 -h480 -w8 -e.80
设置在预组装过程对串联重复进行屏蔽的参数。其参数会传递给HPC.TANmask命令。该命令是正常daligner命令的一个变种，适合对串联重复序列进行比对。
pa_HPCREPmask_option = -k18 -h480 -w8 -e.80
设置在预组装过程对散在重复进行屏蔽的参数。其参数会传递给HPC.REPmask命令。该分析步骤会迭代3次。HPC.TANmask命令进行数据分析时有两个关键参数-g和-c，用于设置分析的数据块数量和覆盖度阈值。pa_REPmask_code设置3次迭代过程中这两个参数的值。这两个参数值的设置与数据分块方法和测序数据量有关。
pa_REPmask_code = 1,666;4,266;8,53
设置三次HPC.REPmask命令中的-g和-c参数的值。以上参数值则表示不对散在重复进行屏蔽。该参数的设置可以参考网址https://dazzlerblog.wordpress.com/2016/04/01/detecting-and-soft-masking-repeats/。在该文章中，作者给出了一个示例：对30Gb的基因组测序30x并将数据分割成3600个大小为250Mb的数据块，使用参数值“1,20;10,15;100,10”。每个数据块的测序深度为0.008x，对1个数据块的数据进行比对后，对reads中覆盖度超过20x的位点进行屏蔽，即对高度重复（20 / 0.008 = 2500x）的基因组序列位点进行屏蔽；然后第二次迭代中，对10个数据块的数据进行比对，对reads中覆盖度超过15x的位点进行屏蔽，即对中度重复（15 / 0.08 = 187.5x）的基因组序列位点进行屏蔽；最后第三次迭代中，对100个数据块的数据进行比对，对reads中覆盖度超过10x的位点进行屏蔽，即对低度重复（10 / 0.8 = 12.5x）的基因组序列位点进行屏蔽。示例中基因组大小是30Gb，超级大，一般含有大量重复序列，使用上述参数对数据中的重复位点进行屏蔽。
对于常规的基因组组装，设置pa_REPmask_code参数后，在三次迭代过程中，推荐逐步对重复次数为1000x、100x和10x的序列位点进行屏蔽。例如：对300Mb大小的基因组测序100x，将数据分割成150个大小为200Mb的数据块，可以考虑设置pa_REPmask_code值为“1,666;4,266;8,53”。
ovlp_HPCTANmask_option = -k20 -h480 -w8 -e.80
设置在OLC组装过程对串联重复进行屏蔽的参数。其参数会传递给HPC.TANmask命令。

预组装（pre-assembly）参数：
length_cutoff = -1
设置对种子序列的长度筛选阈值。若该参数值设置为-1，则程序自动计算种子序列的长度筛选阈值，根据如下两个参数挑选最长的reads序列作为种子序列，直到其数据量达到基因组指定覆盖度为止。
genome_size = 4652500
seed_coverage = 30
设置基因组大小和种子序列覆盖度。推荐设置seed_coverage参数的值为20~40x。当然是用更多的种子序列有利于基因组的组装效果，但是会消耗更多运行时间。
pa_daligner_option = -k14 -w6 -h35 -e.70 -l1000
在新版本FALCON中该参数额外再分出了pa_HPCdaligner_option参数，专门用于对数据块的分配。pa_daligner_option参数专门用于调用daligner软件对raw reads进行比对，进行overlap分析。其常用参数：
-k<int>    default: 14
设置k-mer大小。推荐设置为14~18。设置较低的值，能提高sensitivity，得到更多的reads重叠结果，但增加磁盘、内存、计算和时间消耗，适合数据质量较差的情况；设置较高的值，能提高specificity，系统资源消耗更少，运行速度更快，仅适合数据质量较高的情况。
-w<int>    default: 6
daligner命令对两个reads进行比对，得到一个斜对角线的比对带。该带的宽度设置为2的w次方。
-h<int>    default: 35
在斜对角线带的比对区域上，需要至少有指定数目的碱基能被k-mers覆盖。-k、-w和-h是daligner命令进行比对的主要参数。默认设置下，寻找两两reads之间重叠，先将reads分割成长度为14bp的k-mers序列；对其中一条序列宽度为64bp的区域进行分析，要求至少有35个碱基位点能和被另外一条序列的k-mers覆盖，则找到重叠。
-k、-w和-h的默认参数适合于PacBio测序的raw reads。对于修正后的read，推荐使用参数“daligner -k20 -h60 -e.95”。
-e<float>    default: 0.70
设置identity阈值。推荐设置为0.70~0.80。设置为较高的值有利于单倍型序列的组装。
-l<int>    default: 1000
忽略长度小于低于设定值的reads。
-T<int>    default: 4
设置每个daligner使用的CPU线程数。
-M<int>
尝试仅使用指定大小GB的内存。程序通过忽略覆盖度过高的k-mers来降低内存使用。一般情况下，对
pa_HPCdaligner_option = -v -B4 -M16
程序调用HPC.daligner进行分析，相比于daligner命令，其增加的常用参数：
-v
将-v参数传递给LAsort和LAmerge命令。
-B<int>    default: 4
设置每个daligner任务对指定数目的数据块进行分析。在最新版本的FALCON中，该参数不能再设置到pa_daligner_option参数中，否则程序运行出错。
falcon_sense_option = --output-multi --min-idt 0.70 --min-cov 4 --max-n-read 200
FALCON使用fc_consensus命令根据daligner进行Overlapping分析的结果进行一致性分析，从而对种子序列进行校正。常用参数：
--output-multi
输出每个种子序列多个修正后的区域。
--min-cov    default: 6
设置当种子序列上某位点覆盖度低于指定阈值时，则在该位点对序列打断或截短。
--min-cov-aln    default: 10
设置当种子序列的平均覆盖度低于指定阈值时，过滤掉该序列。
--min-n-read    default: 10
--max-n-read    default: 500
对覆盖度在--min_n_read和--max_n_read两个参数设定范围内的位点进行一致性分析，从而对种子序列进行校正。
--min-idt    default: 0.7
设置比对结果中reads重叠部分的最小identity阈值。
--n-core    default: 24
设置运行的线程数，默认值是24。
falcon_sense_greedy = False
falcon_sense_skip_contained = False

OLC算法进行组装的参数：
length_cutoff_pr = 12000
对预组装结果中长度超过此阈值的reads使用OLC算法进行下一步组装。
ovlp_HPCdaligner_option = -v -B4 -M16
ovlp_daligner_option = -k20 -w6 -h60 -e.96 -l1000
调用daligner对校正后的preads进行overlapping分析时，参数相对要设置更加严格。
-k<int>    default: 14
对preads进行比对，推荐该参数的值为18~24。
-e<float>    default: 0.70
对preads进行比对，推荐该参数的值为0.93~0.96。
overlap_filtering_setting = --max-diff 100 --max-cov 100 --min-cov 2 --bestn 10
过滤overlaps。若reads首尾两端的覆盖度比平均覆盖度大很多，则表明reads首尾是重复序列；若reads首尾两端的覆盖度比平均覆盖度相差较小很多，则表明reads首尾出现错误的可能性很大。需要过滤掉这种reads的overlaps结果。
--bestn    default: 10
设置报告reads此数目的最优overlaps。
--min-cov和—max-cov
所允许的reads首尾的覆盖度范围。
--max-diff
设置所允许的首尾覆盖度值的最大差异。
fc_ovlp_to_graph_option = --min_len 4000 --min_idt 0.96

任务投递与计算资源消耗参数：
[job.defaults]
job_type = local
设置任务提交方式。local表示使用本地主机运行FALCON，运行相对稳定。也可以设置为集群任务提交方式，可以设置的值有：sge、lsf、pbs、torque和slurm。该参数设置不同的值，后面的submit也要设置对应的信息。
pwatcher_type=blocking
submit = bash -C ${CMD} >| ${STDOUT_FILE} 2>| ${STDERR_FILE}
MB=32768
NPROC=6
njobs=32
FALCON软件流程对数据进行了切割，再进行并行化运算，加快程序运行速度。MB设置任务实用的内存数；NPROC设置单个任务的CPU线程数；njobs设置并行任务数。后续分别对每个流程中的计算资源进行分配。da表示预组装过程中运行的daligner命令；la表示预组装过程中运行的Lasort和LAmerge命令；cns表示预组装过程中运行的fc_consensus命令；pda表示运行OLC组装过程中运行的的daligner命令；pla表示OLC组装过程中运行的Lasort和LAmerge命令；asm表示最终组装过程。
[job.step.da]
NPROC=4
MB=32768
njobs=240
[job.step.la]
NPROC=4
MB=32768
njobs=240
[job.step.cns]
NPROC=8
MB=65536
njobs=240
[job.step.pda]
NPROC=4
MB=32768
njobs=240
[job.step.pla]
NPROC=4
MB=32768
njobs=240
[job.step.asm]
NPROC=24
MB=196608
njobs=1

FALCON的结果文件：

0-rawreads/
该目录存放对raw subreads进行overlpping分析与校正的结果；
0-rawreads/cns-runs/cns_*/*/*.fasta存放校正后的序列信息。

1-preads_ovl/
该目录存放对校正后reads进行overlapping的结果；

2-asm-falcon/
该目录是最终结果目录，主要的结果文件是p_ctg.fa和a_ctg.fa。

1. 信号肽简介

信号肽是蛋白质N-末端一段编码长度为5-30的疏水性氨基酸序列，用于引导新合成蛋白质向通路转移的短肽链。信号肽存在于分泌蛋白、跨膜蛋白和真核生物细胞器内的蛋白中。

信号肽指引蛋白质转移的方式有两种：（1）常规的分泌（Sec/secretory）通路；（2）双精氨酸转移（Tat/twin-arginine）通路。前者存在于原核生物蛋白质转移到质膜过程中，以及真核生物蛋白质转移到内质网膜的过程中。后者存在于细菌、古菌、叶绿体和线粒体中，信号肽序列较长、疏水性较弱且尾部区含有两个连续精氨酸。相比于前者转运非折叠蛋白质，后者能转运折叠蛋白质跨越双层脂质膜 。

信号肽指引蛋白质转运后，将由信号肽酶进行切除。信号肽酶有三种：（1）一型信号肽酶（SPaseI）；（2）二型信号肽酶（SPaseII）；（3）三型信号肽酶（SPaseIII）。大部分信号肽由SPaseI进行移除，SPaseI存在古菌、细菌和真核生物中，且在真核生物的内质网膜上仅存在一型信号肽酶。细菌和古菌脂蛋白的信号肽C端含有一段称为 lipobox 的保守区域，由SPaseII切除其信号肽，且lipobox紧邻切除位点（CS/Cleavage Site）的氨基酸是半胱氨酸，这和锚定到膜的功能是相关的。细菌的四型菌毛蛋白信号肽由SPaseIII进行切除。此外：分泌通路（Sec）相关信号肽能由SPaseI、SPaseII和SPaseIII切除，但是双精氨酸转移（Tat）通路相关信号肽仅由 SPaseI和SPaseII切除。

使用SignalP 5.0能对原核生物的信号肽Sec/SPI、Sec/SPII和Tat/SPI，和对真核生物仅含有 Sec/SPI信号肽进行预测 。 SignalP 5.0目前不能对Tat/SPII进行预测。此外，由于没有足够大的数据进行训练，SignalP 5.0 也不能对Sec/SPIII进行分析。

SignalP是最为常用的信号肽分析软件，常用于分泌蛋白预测的第一步。到目前2019年4月，SignalP版本到了第5版。第一版本基于人工神经网络（Artificial Neural Network）；第二版本引入隐马尔科夫模型（Hidden Markov Models）；第三版本改进切除位点（Cleavage site）预测方法；第四版本改进对信号肽和跨膜螺旋的区分能力。这四个版本仅能对Sec/SPI类型的信号肽进行预测。而第五版本结合了深度神经网络（deep neural network）、条件随机分类（Conditional random field classification）和迁移学习（transfer learning）方法，能对信号肽进行更准确的预测。

可以使用SignalP网页工具进行分析。但一次仅支持最多5000条序列。以下讲解本地化运行SignalP软件。

2. SignalP软件下载和安装

# 需要填写edu邮箱和相关信息来获取下载地址
$ wget http://www.cbs.dtu.dk/download/6B91F6BC-5A05-11E9-8172-2ED6B9CD16B5/signalp-5.0.Linux.tar.gz -P ~/software/
$ tar zxf /home/train/software/signalp-5.0.Linux.tar.gz -C /opt/biosoft/
$ echo 'PATH=$PATH:/opt/biosoft/signalp-5.0/bin' >> ~/.bashrc
$ source ~/.bashrc

3. SignalP软件使用

对真核生物的全基因组蛋白序列进行信号肽预测：
$ signalp -batch 30000 -org euk -fasta proteins.fasta -gff3 -mature

signalp命令的参数：

-batch <int>    default: 10000
    程序并行化运行的序列条数。程序能多线程运行，速度很快。推荐设置该参数值大于FASTA文件的序列总条数，虽然增加内存消耗，但能加快程序运行。
-org <string>    default: euk
    设置分析的物种类型。该参数值有4种：arch，古菌；gram+，革兰氏阳性细菌，gram-，革兰氏阴性细菌；euk，真核生物。
-fasta <string>
    输入FASTA格式的蛋白序列文件。
-prefix <string>
    设置输出文件前缀。默认在程序运行目录下输出结果文件，和输入文件名的前缀相同，后缀为_summary.signalp5。
-gff3
    添加该参数，输出GFF3格式的信号肽预测结果。
-mature
    添加该参数，对含有信号肽的蛋白序列，切除信号肽后输出其成熟蛋白序列。可以用于下游的跨膜区分析。
-tmp <string>    default: /tmp
    设置临时文件夹路径。
-format <string>    default: short
    设置输出格式。该参数值有两个：short，仅输出一个信号肽预测的整合文本结果；long，额外输出每条序列的各位点预测文本结果和图片结果； 
-plot <string>    default: png
    设置输出图片结果的类型。当--format参数为long时，该参数生效。该参数值有三个：png；eps；none表示仅得到表格文件。
-version
    打印程序版本并退出。

我的wordpress站点总是有很多垃圾评论。于是需要对相应的IP地址禁止访问。

我的方法是：在/etc/httpd/conf.d目录下生成后缀为多个.conf的文件，每个文件分别对一个目标文件夹进行保护。比如生成文件/etc/httpd/conf.d/deny_chenlianfu.conf，其内容为：

<Directory "/home/chenlianfu/wordpress">
    AllowOverride None
    Options MultiViews FollowSymLinks ExecCGI
    Order allow,deny
    Allow from all
    Deny from 100.42.17.90
    Deny from 101.4.136.34
    Deny from 101.91.215.188
    Deny from 101.98.247.14
    Deny from 95.85.80.227
    Deny from 95.85.80.82
    Deny from 95.85.80.86
    Deny from 98.174.90.36
</Directory>

然后重启httpd服务，使配置文件生效，从而禁止配置文件中的IP访问我的网站。

/etc/init.d/httpd restart

为了对多个文件夹进行保护，则要生成多个配置文件。我编写名为wordpress_deny_ips.pl的程序，输入含有准备禁止访问IP地址的文本文件，自动生成多个配置文件：

#!/usr/bin/perl
use strict;

my $usage = <<USAGE;
USAGE:
    perl $0 web.txt

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

open IN, "/etc/httpd/conf.d/deny_chenlianfu.conf";
my %ip;
while (<IN>) {
    $ip{$1} = 1 if m/(\d+\.\d+\.\d+\.\d+)/;
}
close IN;

open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
my %ip_new;
while (<IN>) {
    $ip_new{$1} = 1 if m/(\d+\.\d+\.\d+\.\d+)/;
}
close IN;

my $number = 0;
foreach (keys %ip_new) {
    if (exists $ip{$_}) {
        print STDERR "Duplicate $_\n";
    }
    else {
        $ip{$_} = 1;
        $number ++;
    }
}

my @num = keys %ip;
my $num = @num;
print STDERR "$number ips were add\n$num ips in total\n";

my $out = '
    AllowOverride None
    Options MultiViews FollowSymLinks ExecCGI
    Order allow,deny
    Allow from all
';
foreach (sort keys %ip) {
    $out .= "    Deny from $_\n";
}
$out .= "</Directory>\n";


open OUT, ">", "/etc/httpd/conf.d/deny_chenlianfu.conf" or die "Can not create file /etc/httpd/conf.d/deny_chenlianfu.conf, $!\n";
print OUT '<Directory "/home/chenlianfu/wordpress">' . $out;
close OUT;

open OUT, ">", "/etc/httpd/conf.d/deny_zhengyue.conf" or die "Can not create file /etc/httpd/conf.d/deny_zhengyue.conf, $!\n";
print OUT '<Directory "/home/zhengyue/wordpress">' . $out;
close OUT;

open OUT, ">", "/etc/httpd/conf.d/deny_wuchangsong.conf" or die "Can not create file /etc/httpd/conf.d/deny_wuchangsong.conf, $!\n";
print OUT '<Directory "/home/wuchangsong/wordpress">' . $out;
close OUT;

很多垃圾评论是全英文或含有日文。可以在主体对应的functinons.php文件中添加一些设置来禁止全英文或包含日文的评论。我在文件./wp-content/themes/twentyeleven/functions.php的尾部添加以下代码：

function refused_english_comments( $incoming_comment ) {
    $pattern = '/[一-龥]/u';
    //禁止全英文评论
    if(!preg_match($pattern, $incoming_comment['comment_content'])) {
        wp_die( "您的评论中必须包含汉字!" );
    }

    $pattern = '/[あ-んア-ン]/u';
    //禁止日文评论
    if(preg_match($pattern, $incoming_comment['comment_content'])) {
        wp_die( "评论禁止包含日文" );
    }
    return( $incoming_comment );
}
add_filter('preprocess_comment', 'refused_english_comments');

然后重启网页服务，即可生效。