使用Virtualbox安装CentOS7系统，在安装增强包工具后，鼠标左键在桌面窗口中失灵。此时，需要在root权限下安装kernel软件，重启后恢复正常。
# yum install -y kernel*
# shutdown -r now
也有人说CentOS7安装完毕默认安装了open-vm-tools软件。将其默认安装的该软件删除后，再安装增强包工具后重启，能解决该问题。
解决Virtualbox安装CentOS7鼠标左键失效情况
回复
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 ~/.bashrc2. 示例数据下载和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 *.typecut -f 2 snoRNA_1.list > aa
paste snoRNA_2.list aa > snoRNA.list
rm aa snoRNA_1.list snoRNA_2.list snoRNA_type.pl3. 编写程序统计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
EOF2. 服务器CentOS系统安装完毕后的基本配置
使用root用户在管理节点和计算节点上进行以下操作：
- 修改/etc/profile.d/perl-homedir.sh配置文件，以免每次登录用户，自动在家目录下生成perl5文件夹
- 修改/etc/sudoers配置文件，将自己的用户（例如 train）变成超级管理员用户
- 修改/etc/selinux/config配置文件，永久关闭linux的一个安全机制，开启该安全机制会对很多操作造成阻碍。
- 修改/etc/ssh/sshd_config配置文件，使openssh远程登录更安全，更快速
- 增加系统资源权限
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
EOF3. 配置集群中各服务器的主机名和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
EOF4. 将控制节点的以太网共享给计算节点
控制节点通过电信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 -a6. 在集群计算机上创建新用户
首先，生成文件/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
donecreate_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使用wol对局域网内服务器通过网络唤醒方法进行远程开机
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值。特别是当控制服务器有多个网口连接时，不设置该参数的值，或设置错误的值，可能会导致不能唤醒目标服务器。
infiniband网卡设置
在CentOS7部署SGE
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.service2. 从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
删除任务使用eggNOG数据库进行基因功能注释
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 eggNOGemapper.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结果解读和分类图绘制
使用PHI-base数据库进行致病基因分析
1. PHI-base数据库简介
PHI-base数据库从文献中收集经过实验验证了的致病基因和效应基因的序列。目前（20190411）数据库版本为4.6版本，从3011篇文献中收集了263种致病菌（细菌、真菌、原生动物和线虫）的6438个基因和194种宿主（植物占~70%、脊椎动物、昆虫、线虫和真菌）的11340种相互关系，其中包含510中疾病。PHI-base将收集到的参考文献信息、 基因信息、病原和宿主信息、疾病信息、表型和相互关系等记录到数据库中，并提供关键词进行搜索。同时，PHI-base还提供BLAST网页工具，用于将序列比对到数据库中，从而根据BLAST结果进行致病基因分析。此外，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