使用SignalP对蛋白序列进行信号肽预测

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
打印程序版本并退出。

网页服务器禁止IP访问

我的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');

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

检测全基因组序列中的端粒序列

使用三代测序数据能获得较好的、甚至完整的基因组序列。通过检测基因组序列两端的属于端粒(Telomere)的特定的重复序列,可以知道基因组组装是否得到完整的染色体水平的序列。若在序列中间检测到端粒序列,可以知道基因组组装过程中对Contigs有错误的连接。

在染色体序列首尾存在端粒序列。在人类中,端粒序列由重复单元TTAGGG,串联重复约2500次组成。 不同的物种,端粒的重复单元可能不一样,可以在端粒数据库中查询。 我对一种子囊菌使用PacBio测序数据进行了基因组组装,得到12条序列,发现大部分序列首尾均出现了TTAGGG/CCCTAA的串联重复序列。我对一种昆虫物种PacBio组装基因组序列进行分析,发现序列首尾出现一些长度较长微卫星重复序列,而没有固定的重复单元,这和端粒数据库中的结果一致。

编写Perl程序对全基因组序列的端粒序列进行搜索,查看基因组序列的完整情况:

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

my $usage = <<USAGE;
Usage:
    $0 genome.fasta > telomere_info.txt

    大部分物种端粒序列的重复单元是TTAGGG/CCCTAA。本程序能在基因组中寻找端粒重复单元的串联重复序列,并给出位点信息。

    --split-length <int>    default: 100000
    --overlap-length <int>    default: 10000
    程序会将每条序列打断后进行重复单元搜索。这两个参数设置打断的序列长度和相邻两序列之间的重叠长度。

    --repeat-unit <string>    default: TTAGGG
    设置重复单元碱基序列,该重复单元的反向互补也将作为重复单元进行搜索。可以在端粒数据库(http://telomerase.asu.edu/sequences_telomere.html)中寻找目标端粒重复单元。
    vertebrate sp.      TTAGGG
    plants sp.          TTTAGGG
    Pezizomycotina      TTAGGG

    --min-repeat-num <int>    default: 4
    设置重复单元最小重复次数。


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

my ($splitLength, $overlapLength, $repeatunit, $minRepeatNum);
GetOptions(
    "split-length:i" => \$splitLength,
    "overlap-length:i" => \$overlapLength,
    "repeat-unit:s" => \$repeatunit,
    "min-repeat-num:s" => \$minRepeatNum,
);
$splitLength ||= 100000;
$overlapLength ||= 10000;
$repeatunit ||= "TTAGGG";
$repeatunit = uc($repeatunit);
my $repeatunit_reverse = reverse $repeatunit;
$repeatunit_reverse =~ tr/ATCG/TAGC/;
$minRepeatNum ||= 4;

# 读取基因组序列
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
my (%seq, $seq_id);
while (<IN>) {
    chomp;
    if (m/^>(\S+)/) {
        $seq_id = $1;
    }
    else {
        $_ = uc($_);
        $seq{$seq_id} .= $_;
    }
}
close IN;

# 将基因组序列打断
my (%seq_split, %seq_length);
foreach my $id (keys %seq) {
    my $seq = $seq{$id};
    my $length = length($seq);
    $seq_length{$id} = $length;
    my $pos = 0;
    while ($pos < $length) {
        $seq_split{$id}{$pos} = substr($seq, $pos, $splitLength + $overlapLength);
        $pos += $splitLength;
    }
}

print "SeqID\tSeqLength\tStart\tEnd\tLength\tType\n";
foreach my $id (sort keys %seq_split) {
    foreach my $pos (sort {$a <=> $b} keys %{$seq_split{$id}}) {
        my $seq = $seq_split{$id}{$pos};
        while ($seq =~ m/(($repeatunit){$minRepeatNum,})/g) {
            my $length = length($1);
            my $end = pos($seq);
            $end = $end + $pos;
            my $start = $end - $length + 1;
            print "$id\t$seq_length{$id}\t$start\t$end\t$length\t$repeatunit\n";
        }
        while ($seq =~ m/(($repeatunit_reverse){$minRepeatNum,})/g) {
            my $length = length($1);
            my $end = pos($seq);
            $end = $end + $pos;
            my $start = $end - $length + 1;
            print "$id\t$seq_length{$id}\t$start\t$end\t$length\t$repeatunit_reverse\n";
        }
    }
}

使用DIAMOND将全基因组蛋白序列比对到Nr数据库

1. DIAMOND简介

对全基因组的基因进行Nr注释是必不可少的一步。由于Nr数据库非常大,导致使用BLAST会消耗巨大的计算资源和时间。使用DIAMOND则能快500-20000倍,而获得和BLAST比较一致的结果。特别是对于长度为100-150bp,数量超过1M的核酸序列,DIAMOND的速度比BLASTX快20000倍;当e-value设置为1e-5时,DIAMOND能报告80~90%的BLASTX匹配结果;若设置sensitivie模式,DIAMOND速度是BLASTX的2500倍,DIAMOND能报告94%的BLASTX匹配结果。因此,对大批量的序列进行Nr注释或规模较大的蛋白数据库的比对,使用DIAMOND是优先选择的方式。

虽然DIAMOND软件的运行速度快,但是其结果中没有匹配的Subject序列的描述信息。若需要蛋白名称和物种的描述信息,需要自行编写软件从数据库的FASTA文件中提取并整合到DIAMOND结果中。

2. DIAMOND下载和安装

$ wget https://github.com/bbuchfink/diamond/releases/download/v0.9.24/diamond-linux64.tar.gz -P ~/software/
$ mkdir /opt/biosoft/diamond
$ tar zxf ~/software/diamond-linux64.tar.gz -C /opt/biosoft/diamond
$ echo 'PATH=$PATH:/opt/biosoft/diamond' >> ~/.bashrc
$ source ~/.bashrc

3. DIAMOND使用

DIAMOND软件的主命令是diamond,它的使用包含几个子命令。DIAMOND最常用的使用方法:

使用DIAMOND软件的子命令makedb将FASTA格式的蛋白序列创建成后缀为dmnd的数据库文件:
$ diamond makedb --in nr_eukaryon.fasta -d nr_eukaryon_20190405
    将Nr数据库(版本20190405)的真核子集创建成DIAMOND数据库,数据库包含37.4M条蛋白序列,服务器是24线程2.0GHz的CPU资源,耗时17min。

将Illumina测序数据比对到数据库:
$ diamond blastx --db nr_eukaryon_20190405 --query reads.fq.gz --out reads.tab

将全基因组蛋白序列比对到数据库:
$ diamond blastp --db nr_eukaryon_20190405 --query proteins.fasta --out nr.tab --outfmt 6 --sensitive --max-target-seqs 20 --evalue 1e-5 --id 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1
    对18700个蛋白序列进行Nr注释,服务器160线程2.4GHz的CPU资源,耗时91min(不加--sensitive参数只需要16min)。

此外,DIAMOND还包含其它几个子命令:

blastx
    将核酸序列比对到蛋白序列数据库
view
    从DAA文件生成格式化的输出结果
getseq
    从DIAMOND数据库文件获取FASTA序列
version
    打印DIAMOND软件的版本号
dbinfo
    打印数据库的信息
help
    打印DIAMOND软件的详细帮助信息

DIAMOND的参数:

数据库参数:
--in <string>    default: STDIN
    输入FASTA格式的蛋白序列数据库文件。
--db | -d <string>
    设置数据库文件路径和前缀。创建数据库时,会生成一个后缀为.dmnd的数据库文件。比对时,则是输入相应的数据库文件。
--taxonmap <string>
--taxonnodes <string>
    这两个参数输入NCBI的物种分类信息。前者输入Nr数据库中蛋白编号和物种编号的对应关系信息( ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz);后者输入物种编号的层次信息(ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdmp.zip)。输入这两个文件有利于DIAMOND和Nr数据库进行比对时,可以选择和数据库中指定的物种大类的数据进行比对。

输入参数:
--db | -d <string>
    设置数据库文件路径和前缀。创建数据库时,会生成一个后缀为.dmnd的数据库文件。比对时,则是输入相应的数据库文件。
--query | -q <string>    default: STDIN
    输入需要注释的FASTA或FASTQ格式的序列文件。可以是带.gz后缀的压缩文件。
--taxonlist <string>
    输入NCBI分类编号,仅对数据库中的目标子集进行比对。可以输入多个使用逗号分隔的编号ID。例如,古菌(2157)、细菌(2)、病毒(10239)、真菌(4751)、动物(33208)和植物(3193)。使用该参数,必须要使用--taxonmap和--taxonnodes参数构建数据库。
--query-gencode <int>    default:1
    设置在进行BLASTX模式进行比对时所使用的遗传密码。
--strand <string>    default: both
    设置对query序列的那条链进行比对。可以设置的值为:boh、plus和minus。
--min-orf | -l <int>
    在使用BLASTX模式进行比对时,若序列的某个ORF低于此值,则忽略该ORF。默认设置下:若核酸序列长度低于30,则值为1;若核酸序列长度低于100,则值为20;若核酸序列长度不低于100,则值为40。若该值设为1,则表示在BLASTX模式时,识别所有长度的ORF并进行比对。

比对参数:
--sensitive
    添加该参数,则能得到更多比对结果。该模式适合比对较长的序列。默认模式主要适用于比对short reads序列(Illumina reads),搜寻比对长度为30~40aa且bit得分大于50的匹配结果。
--more-sensitive
    相比于sensitive,能得到更全的比对结果。
--frameshifit | -F <int>
    在进行BLASTX模式时,设置对移码匹配的罚分,且一般推荐设置该参数值为15左右。设置该参数后,有利于比对较长的有INDEL错误的reads序列。默认情况是是禁用了此功能。
--gapopen <int>    default: BLOSUM62矩阵的默认值为1
    设置打开Gap罚分。
--gapextend <int>    default: BLOSUM62矩阵的默认值为1
    设置延长Gap罚分。
--matrix <string>    default: BLOSUM62
    设置计分矩阵。可以设置的值有:BLOSUM45,BLOSUM50,BLOSUM62,BLOSUM80,BLOSUM90,PAM250,PAM70和PAM30。
--comp-based-stats <bool>    default: 1
    设置是否对比对得分进行校正。
--algo <bool>
    设置种子序列搜索的算法。0表示使用double-indexed算法,1表示query-indexed算法。第一种算法是软件的主要算法,但不适合query数据较少的情况,这时需要使用第二种算法。默认设置下,软件会根据query和数据库的数据量自动选择相应的算法。两种算法的结果基本一致。
--freq-sd <int>
    设置种子序列最大的匹配频率。若种子序列在query或数据库种出现的频率过高,则忽略之。在一般模式下,该参数值为50;在sensitive模式下,该参数值为100,在more-sensitive模式下,该参数值为200。


输出参数:
--out | -o <string>    default:STDOUT
    设置输出文件。
--outfmt | -f <int>    default: 6
    设置输出格式。支持的格式有:0,两两比对格式;5,XML格式;6,BLAST表格格式;100,DIAMOND匹配存档(DDA)格式,该格式可以使用diamond view命令转换成其它格式;101,SAM格式;102,分类鉴定结果,该模式下结果文件分三列,QueryID、NCBI物种分类ID、最佳匹配evaule;103,PAF格式。
--salltitles
    在DDA格式结果中包含all subject titles。
--sallseqid
    在DDA格式结果中包含all subject ids。由于id比title要短,相比于--salltitles,本参数输出的结果文件更小。
--compress <bool>    default: 0
    设置是否对暑促文件进行压缩。0表示不压缩,1表示进行gzip压缩。
--max-target-seqs | -k <int>    default: 25
    设置每个query比对结果的最大匹配数目。
--top <int>
    若设置该参数的值,则报告一定数目的top匹配结果。例如设置该参数值为10,则报告匹配得分为top10的结果,即报告所有得分和最高得分相差10%以内的匹配结果。该参数会取代--max-target-seqs参数设置的值。
--range-culling
    添加该参数能剔除匹配到query序列相同位置的hit结果。该参数适合于较长的query DNA序列。当DNA序列较长时,它可能跨越多个基因,而一般情况下仅报告25个最佳匹配结果,若这25个hits结果可能evaule值很好但得分较低,且一起不能覆盖query序列的50%,会导致阵势的比对结果较差。因此,添加该参数能剔除匹配到query序列相同位置的hit结果,以使最终的比对结果对query序列的覆盖率超过50%。
--range-cover <float>    default: 50.0
    和参数--range-culling配合使用,表示剔除hits以使所有的比对结果能对query序列的覆盖率超过此参数设定的百分比。
--long-reads
    对于long read DNA测序数据,推荐添加该参数,等同参数--range-culling --top 10 -F 15。
--evalue | -e <float>    default: 0.001
    设置比对的evalue阈值。
--min-score <int>
    设置比对的bit score阈值。设置该参数会替代--evalue参数。
--id <int>
    设置identity阈值。
--query-cover <int>
    设置对query序列的覆盖度阈值。该阈值是HSP的阈值。
--subject-cover <int>
    设置对subject序列的覆盖度阈值。该阈值是HSP的阈值。
--max-hsps <int>
    设置报告的HSPs数量。默认是报告所有的HSPs,若设置成1,则仅报告最优HSP。
--unal <bool>
    设置是否报告没有匹配的比对结果。默认设置下,--outfmt为1、5和101时,该参数值为1,其它输出方式该参数值为0。
--no-self-hits
    添加该参数,则不报告序列和自己的比对结果。

性能参数:
--threads | -p <int>    default: Max
    设置程序运行所使用的CPU线程数。默认是服务器可用的最大CPU线程数。
--block-size | -b <float>    default: 2.0
    设置每次处理多少G的序列字符数。该参数控制程序消耗的内存量,一般内存消耗该值的6倍。设置更大的值,会消耗更多的内存和临时磁盘空间,但能提高性能。默认设置下,程序一次处理2G个序列字符,消耗内存12G。
--tmpdir | -t <string>    default: directory of --out
    设置临时文件夹路径。推荐该文件所在剩余磁盘空间有100G及以上。若将该参数设置为/dev/shm,则会将临时文件存放在内存,会增加内存消耗和计算性能。
--index-chunks | -c <int>    default: 4
    将seed index分成指定的份数。推荐将该参数值设置成1,能增加计算性能和内存使用量。

DIAMOND软件运行过程中,若使用参数 –block-size 20.0 –tmpdir /dev/shm –index-chunks 1 能使用较大内存,不会对数据进行分块,程序运行较快。程序运行过程中可能完全不会生成过程文件。但在程序运行过程中给出一些进度信息。当shape到15的时候,则表示程序快运行完毕。通过查看shape值来确定程序运行的进度。

CPU threads: 160
Scoring parameters: (Matrix=BLOSUM62 Lambda=0.267 K=0.041 Penalties=11/1)
Temporary directory: /dev/shm
Opening the database…
Target sequences to report alignments for: 20
Opening the input file…
Opening the output file…
Loading query sequences…
Masking queries…
Building query seed set…
Algorithm: Double-indexed
Building query histograms…
Allocating buffers…
Loading reference sequences…
Building reference histograms…
Allocating buffers…
Initializing temporary storage…
Processing query chunk 0, reference chunk 0, shape 0, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
Processing query chunk 0, reference chunk 0, shape 1, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
...
Processing query chunk 0, reference chunk 0, shape 15, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
Deallocating buffers…
Computing alignments…
Deallocating reference…
Loading reference sequences…
Deallocating buffers…
Deallocating queries…
Loading query sequences…
Closing the input file…
Closing the output file…
Closing the database file…
Deallocating taxonomy…
Total time = 659.515s
Reported 33888205 pairwise alignments, 33890340 HSPs.
122518 queries aligned.

DIAMOND默认设置下输出表格格式的结果。结果分12列,其结果信息和BLAST默认设置-outfmt 6输出的格式完全一致。

1. query序列ID
2. subject序列ID
3. Identity百分比
4. 匹配长度
5. 错配长度
6. 打开Gap的次数
7. query序列起始位点
8. query序列结束位点
9. subject序列起始位点
10. subject序列结束位点
11. E-vaule值
12. bitscore得分

CentOS7系统配置firewalld使用NAT方法将网络共享给内部局域网机器

学校实验室有台服务器申请了固定的公网IP,能连接外部网络,同时该机器和其它几台内部服务器连接在一个路由器上。需要将该服务器的网络共享给其它内网服务器。进行如下设置即可。

首先,外网服务器有两根网线连接,一根和外网相连,一根和内网路由器相连。在外网服务器进行NAT转发设置:

第一种方法:
开启NAT转发
# firewall-cmd --permanent --zone=public --add-masquerade
开放DNS使用的53端口,否则可能导致内网服务器虽然设置正确的DNS,但是依然无法进行域名解析。
# firewall-cmd --zone=public --add-port=53/tcp --permanent
重启防火墙
# systemctl restart firewalld.service
# 检查是否允许NAT转发
firewall-cmd --query-masquerade
# 关闭NAT转发
firewall-cmd --remove-masquerade

第二种方法:
开启ip_forward转发
在/etc/sysctl.conf配置文件尾部添加
net.ipv4.ip_forward=1
然后让其生效
# sysctl -p
执行firewalld命令进行转发:
firewall-cmd --permanent --direct --passthrough ipv4 -t nat -I POSTROUTING -o enoxxxxxx -j MASQUERADE -s 192.168.1.0/24
注意enoxxxxxx对应外网网口名称。
重启防火墙
# systemctl restart firewalld.service

设置外网服务器的IP地址为192.168.1.1,然后,在内网服务器中对配置文件/etc/sysconfig/network-scripts/ifcfg-enxxxxxx进行修改,设置网卡的网关为外网服务器的IP地址:

ONBOOT=yes
IPADDR=192.168.1.2
PREFIX=24
GATEWAY=192.168.1.1
DNS1=119.29.29.29

重启内网服务器的网络服务,即可访问外网了。

# systemctl restart network.service

创建NCBI的Nr数据库的子集数据库

1. 为什么要做Nr子集数据库

NCBI官网仅提供Nr全数据库。该数据库太大,将物种的蛋白序列使用Blastp比对到Nr数据库非常消耗计算和时间。对1个蛋白序列可能需要1个CPU计算半个小时。若对全基因组2万个基因分析,普通台式机8个CPU要计算2000/(2*8*24)=52天。这主要是由于Nr数据库太大导致的。为了能尽快得到Nr注释结果,可以按物种分类将Nr数据库分割成子集数据库,能得到更快的比对速度。以下是创建Nr子集数据库的步骤。

2. 创建Nr子集数据库的步骤

从NCBI下载Nr数据库和分类数据库文件

cd /opt/biosoft/Nr_database
# 下载Nr数据库(FASTA文件)
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /blast/db/FASTA/nr.gz ./

# 下载NCBI的分类数据库文件
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /pub/taxonomy/taxdump.tar.gz ./
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /pub/taxonomy/accession2taxid/prot.accession2taxid.gz ./

# 解压缩两个NCBI的分类数据库文件
gzip -dc prot.accession2taxid.gz > prot.accession2taxid
mkdir ~/.taxonkit
tar zxf taxdump.tar.gz -C ~/.taxonkit
# 其主要有效文件有两个:
# names.dmp 记录物种名及其分类编号
# nodes.dmp 记录分类编号的节点信息
# 查看~/.taxonkit/names.dmp文件,使用关键词检索得到目标类的分类编号,例如:
# fungi 4751             # grep -P "\|\s+[fF]ungi\w*\s*\|" ~/.taxonkit/names.dmp
# plants 3193            # grep -P "\|\s+[pP]lant\w*\s*\|" ~/.taxonkit/names.dmp
# animals 33208          # grep -P "\|\s+[aA]nimal\w*\s*\|" ~/.taxonkit/names.dmp

下载并安装NCBI分类数据库解析软件TaxonKitTaxonKit,并解析nodes.dmp文件的物种节点信息,得到指定类的所有物种列表信息。

# 下载并安装NCBI分类数据库解析软件TaxonKit
wget https://github.com/shenwei356/taxonkit/releases/download/v0.2.4/taxonkit_linux_amd64.tar.gz
tar zxvf taxonkit_linux_amd64.tar.gz

# 提取古菌(2157)、细菌(2)和病毒(10239)这几个大类下的所有物种编号。
./taxonkit list -j 8 --ids 2,2157,10239 > sub.meta.list

# 再编写程序extract_sub_data_from_Nr.pl获得列表中物种在Nr数据库中的序列信息。
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.meta.list --acc2taxid prot.accession2taxid - > nr_meta.fasta

提取fungi/plants/animals子集

./taxonkit list -j 8 --ids 4751 > sub.fungi.list
./taxonkit list -j 8 --ids 3193 > sub.plants.list
./taxonkit list -j 8 --ids 33208 > sub.animals.list
./taxonkit list -j 8 --ids 10239 > sub.virus.list
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.fungi.list --acc2taxid prot.accession2taxid - > nr_fungi.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.plants.list --acc2taxid prot.accession2taxid - > nr_plants.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.animals.list --acc2taxid prot.accession2taxid - > nr_animals.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.virus.list --acc2taxid prot.accession2taxid - > nr_virus.fasta

使用makeblastdb创建blast本地数据库

makeblastdb -in nr_fungi.fasta -dbtype prot -title nr_fungi -parse_seqids -out nr_fungi_`date +%Y%m%d` -logfile nr_fungi_`date +%Y%m%d`.log
makeblastdb -in nr_plants.fasta -dbtype prot -title nr_plants -parse_seqids -out nr_plants_`date +%Y%m%d` -logfile nr_plants_`date +%Y%m%d`.log
makeblastdb -in nr_animals.fasta -dbtype prot -title nr_animals -parse_seqids -out nr_animals_`date +%Y%m%d` -logfile nr_animals_`date +%Y%m%d`.log
makeblastdb -in nr_virus.fasta -dbtype prot -title nr_virus -parse_seqids -out nr_virus_`date +%Y%m%d` -logfile nr_virus_`date +%Y%m%d`.log

CentOS7的服务开启方法和防火墙管理方法

CentOS6升级到CentOS7后,最明细的改变就是服务的启动方式和防火墙管理方式改变了。

1. CentOS的服务开启方法

在CentOS6中的服务开启方式(以httpd服务为例):

# /etc/init.d/httpd start      启动服务
# /etc/init.d/httpd stop       停止服务
# /etc/init.d/httpd restart    重启服务

# chkconfig httpd on           设置服务开机启动
# chkconfig httpd off          设置服务开机不启动
# chkconfig --list httpd       检测服务是否开机启动

在CentOS7中的服务器开启方式(以httpd服务为例):

# systemctl start httpd.service      启动服务
# systemctl stop httpd.service       停止服务
# systemctl restart httpd.service    重启服务

# systemctl enable httpd.service     设置服务开机启动
# systemctl disable httpd.service    设置服务开机不启动
# systemctl is-enabled httpd.service 检测服务是否开机启动

我常用的服务有:

防火墙服务: iptables firewalld.service
ssh远程联机:sshd     sshd.service    22端口
Apache网站: httpd    httpd.service   80端口
Mysql数据库:mysqld   mariadb.service 3306端口
网络服务:   network  network.service

CentOS7的服务都采用systemctl命令来进行操作。仅有network服务依旧可以和CentOS6一样进行操作。

2. CentOS的防火墙端口开放方法

在CentOS6中仅需要修改/etc/sysconfig/iptables配置文件内容,即可开放相应的端口:

-A INPUT -m state --state NEW -m tcp -p tcp --dport 22 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 80 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 3306 -j ACCEPT

在CentOS7中需要使用firewall-cmd命令对端口进行管理:

查看已经开放的端口
# firewall-cmd --list-ports
添加端口
# firewall-cmd --add-port=80/tcp --permanent
加入--permanent参数,使永久生效。
重启防火墙服务
# systemctl restart firewalld.service
重启后,再查看端口,则生效了。
设置服务开机启动
# systemctl enable firewalld.service

在CentOS7下安装chrome谷歌浏览器

CentOS7可以支持chrome浏览器的安装,而CentOS6不行。生信软件中HGAP4必须要使用chrome浏览器进行操作。

安装chrome浏览器方法:

首先添加google的yum源。
# echo '[google-chrome]
name=google-chrome
baseurl=http://dl.google.com/linux/chrome/rpm/stable/$basearch
enabled=1
gpgcheck=0
#gpgkey=https://dl-ssl.google.com/linux/linux_signing_key.pub' > /etc/yum.repos.d/google-chrome.repo

然后yum安装chrome软件
# yum install google-chrome-stable

安装完毕,即可在Applications – Internet – Google Chrome打开谷歌浏览器。

在CentOS7下使用VNC提供远程桌面服务

1. 为什么我要使用VNC

我使用VNC,是为了能远程连接到服务器的桌面。有些特殊的操作必须在服务器的桌面上完成。比如我需要通过百度网盘上传和下载高通量测序数据以及分析结果。这个时候,必须要安装windows虚拟机,而使用虚拟的windows系统必须要在服务器CentOS系统的桌面界面进行操作。

此外,我也使用了简单易用的teamviewer软件,但是该软件存在问题是,整个服务器系统仅能一个用户使用teamviewer软件。当我正在使用teamviewer的时候,某天,另外一个管理员人直接在服务器上登录了另外的账户桌面,我可能就无法再使用teamviewer连接服务器了。而VNC应该可以多用户同时远程访问各自的桌面。

以下讲解如何在家里通过一台公网服务器作为跳板,访问内网服务器的CentOS7系统桌面。

1. 在内网服务器安装VNC软件并进行设置

使用root用户联网安装软件
# yum install vnc vnc-server

开启5901、5902和5903端口。每增加一个用户使用VNC,则端口号加1,需要额外开放端口。
# firewall-cmd --zone=public --add-port=5901/tcp --permanent
# firewall-cmd --zone=public --add-port=5902/tcp --permanent
# firewall-cmd --zone=public --add-port=5903/tcp --permanent
重启防火墙,使开放的端口生效
# systemctl restart firewalld.service
查看被开启的端口
# firewall-cmd --list-ports

使用普通用户启动VNC软件
$ vncserver
Password:   设置使用VNC连接服务器的密码
Verify:     再次输入一遍密码
Would you like to enter a view-only password (y/n)? 选择n,以使连接到服务器后可以操作服务器
启动VNC后,在用户家目录下会生成 .vnc 配置文件夹,查看其中的log文件,得到启用的端口号,比如我启动了5901端口。

再使用普通用户在内网服务器上构建反向隧道
$ ssh -N -f -R 4499:localhost:5901 chenlianfu@115.29.105.12
我在aliyun上有一台服务器,以之为跳板,以上命令将内网服务器的5092端口映射到了跳板服务器的4499端口。

查看vncserver的服务情况
$ vncserver -list
关闭编号为:1的vncserver
$ vncserver -kill :1
修改vncserver的配置文件~/.vnc/config,尾部添加分辨率信息
geometry=1400x900
再次开启vncserver以获得更好的分辨率支持
$ vncserver

2. 在跳板服务器上的设置

在跳板服务器的防火墙设置中开放4499端口
# firewall-cmd --zone=public --add-port=4499/tcp --permanent
# systemctl restart firewalld.service

在跳板服务器的sshd配置中设置GatewayPorts参数值为yes
# perl -p -i -e 's/.*GatewayPorts.*/GatewayPorts yes/' /etc/ssh/sshd_config
# systemctl restart sshd.service

3. 在家里的笔记本电脑上连接

下载VNC viewer软件,安装完毕后,输入115.29.105.12:4499地址,再输入密码,即可连接到内网服务器的CentOS7桌面。默认画面质量可能较差,可以点击软件的设置图标,在Option选项卡的General栏下的Picture quality菜单中选择High,即可让画面显示清晰。

使用sequenceserver搭建本地化blast网页服务

1. 简介

现在NCBI不再提供wwwblast软件下载。且NCBI的wwwblast软件很久没更新,界面简陋,因此不推荐使用该软件进行本地化blast分析。使用sequenceserver部署本地化blast网页服务,界面简单易用,还可以同时勾选多个数据库进行比对分析。以下讲解使用sequenceserver在CentOS7上搭建本地化blast网页服务。

2. 安装软件

使用root用户安装sequenceserver软件。该软件使用Ruby编程,使用Rubygem即可安装软件。在CentOS7系统中需要额外安装ruby-devel软件。

# yum install ruby-devel
# gem install sequenceserver

3. 使用sequenceserver

使用sequenceserver要求安装有NCBI-blast+程序。

# sequenceserver -d /opt/biosoft/wwwblast/db -n 16 -p 8080
程序运行后,在浏览器中打开 localhost:8080 即可访问本地化blast网页服务。

常用参数:
-b | --bin 
    输入ncbi-blast+的bin文件路径
-d | --database_dir
    输入blast数据库路径
-n | --num_threads
    设置使用的CPU线程数
-p | --port
    设置使用的端口号,默认为4567。