纤维素,半纤维素和果胶的成份及其降解酶

1. 纤维素

Cellulose is a dominant structural polysaccharide in plants composed ofβ -D-glucose units with β-1,4-linkages.

Cellulose decomposition requires multiple enzymes. In general, cellulose is degraded to cellodextrin or cellobiose by the synergistic action of two cellulases: endoglucanase (EC 3.2.1.4) and cellobiohydrolase (EC 3.2.1.91) (Tomme et al., 1995; Bayer et al., 1998). Degradation of cellodextrin or cellobiose into monomeric glucose units requires another enzyme, β-glucosidase (EC 3.2.1.21), that hydrolyzes non-reducing 1,4-linked-β-glucose (Henrissat et al., 1989).

2. 半纤维素

Cellulose fibers are cross-linked by other polysaccharides called `hemicelluloses’ to increase the physical strength of the cell wall. Hemicelluloses include xylan (β-D-xylose units with β-1,4-linkages), glucomannan (β-D-mannose units andβ -D-glucose units with β-1,4-linkages), xyloglucan (β-D-glucose units with β-1,4-linkages, andβ -D-xylose and β-D-glucose units withβ -1,6-linkages), 1,3-1,4-β-glucan (β-D-glucose units with β-1,3- and β-1,4-linkages), and a relatively small amount of other polysaccharides composed of β-D-glucose,β -D-xylose, β-D-mannose and other sugar units with various linkages (McNeill et al., 1984).

3. 果胶

The scaffold of cellulose and hemicelluloses is filled with pectin (α-D-galacturonic acid units with mainly α-1,4-linkages), which functions as a cement-like substance in the cell wall.

reference:
Sakamoto, Kentaro, and Haruhiko Toyohara. “A comparative study of cellulase and hemicellulase activities of brackish water clam Corbicula japonica with those of other marine Veneroida bivalves.” Journal of Experimental Biology 212.17 (2009): 2812-2818.

通过WIG格式将转录组数据展示到Gbrowse2中

1. WIG格式介绍

WIG格式(Wiggle Track Format),可用于将转录组数据进行可视化展示。bigWig格式则是WIG格式的二进制方式,可以使用wigToBigWig将WIG格式转换成BigWig格式。
一个 WIG 格式实例文件:

track type=wiggle_0 name="sampleA1" description="RNA-Seq read counts of species A"
variableStep chrom=chr01 span=10
10001    13
10011    15
10021    12
fixedStep chrom=chr01 start=100031 step=10 span=10
17
15
20

如上例子,有2个注意点:

1. 第一行必须如理示例中格式。只有name和description这两个参数的值可以随意填写。
2. 有两种方法进行数据描述。分别是variableStep和fixedStep。前者数据内容用2行表示,后者数据部分仅用1行表示。
3. 这两种方法的几个参素意义为:
    chrom    设置序列名
    start    fixStep中Locus的起始位置
    step     fixStep中Locus的步进
    span     一个数据对应碱基数目

2. 将Bam文件转换成WIG文件并进行压缩

使用bam2wig命令将bam文件转换成wig文件。bam2wig命令可以来自于Augustus软件。

$ bam2wig sampleA1.tophat.bam > sampleA1.wig

该wig文件的span参数值为1。因此,当基因组越大,则wig文件越大。可以考虑设置span的值为10,能有效减小wig文件的大小。例如编写如些perl程序进行压缩wig文件。

#!/usr/bin/perl
use strict;

my $usage = <<USAGE;
Usage:
    perl $0 RNA-Seq.wig > RNA-Seq.cutdown.wig
USAGE
if (@ARGV==0){die $usage}

open IN, $ARGV[0] or die $!;

$_ = <>;
print;

my $locus = 1;
my $count = 0;
while () {
    if (m/^variableStep/) {
        $count = int(($count + 0.5) / 10);
        print "$locus\t$count\n" if $count > 0;
        s/$/ span=10/;
        print;
        $locus = 1;
    }
    else {
        if (m/(\d+)\s+(\d+)/) {
            my ($num1, $num2) = ($1, $2);
            if ($num1 >= $locus + 10) {
                $count = int(($count + 0.5) / 10);
                print "$locus $count\n" if $count > 0;
                $locus = $num1;
                $count = 0;
            }
            $count += $num2;
        }
    }
}

3. 将wig文件转换成wig binary文件和一个gff3文件

使用Gbrowse2所带命令 wiggle2gff3.pl 将wig文件转换成wig binary文件和一个gff3文件。每个基因组序列得到一个二进制格式的wig文件。同时生成一个gff3文件。该gff3文件指向所有的wig binary文件。

$ mkdir $PWD/gbrowse_track_of_RNA_seq
$ wiggle2gff3.pl --source=sampleA1 --method=RNA_Seq --path=$PWD/gbrowse_track_of_RNA_seq --trackname=track_A1 sampleA1.wig > sampleA1.gff3

4. 导入gff3文件到数据库,并配置Gbrowse配置文件

导入gff3文件

$ bp_seqfeature_load.pl -a DBI::mysql -d gbrowse2_species -u train -p 123456 sampleA1.gff3

配置文件:

[RNA_Seq_sampleA1_xyplot]
feature        = RNA_Seq:sampleA1
glyph          = wiggle_xyplot
graph_type     = boxes
height         = 50
scale          = right
description    = 1
category       = RNA-Seq:sampleA1
key            = Transcriptional Profile

[RNA_Seq_sampleA1_density]
feature        = RNA_Seq:sampleA1
glyph          = wiggle_density
height         = 30
bgcolor        = blue
description    = 1
category       = RNA-Seq:sampleA1
key            = Transcriptional Profile

DISCOVAR的使用

1. DISCOVAR简介

DISCOVAR 是有 ALLPATHS-LG 软件开发团队做出来的软件。主要用于利用 PE 250bp 数据与参考基因组的比对结果,对基因组进行 Variants calling 的同时,进行基因组的组装。特别是近期公布的 DISCOVAR de novo (experimental) 还能进行基因组的 De novo 组装。

2. DISCOVAR的下载和安装

2.1 DISCOVAR的下载和安装

此软件的安装需要GCC 4.7或以上版本。

$ wget ftp://ftp.broadinstitute.org/pub/crd/Discovar/latest_source_code/LATEST_VERSION.tar.gz
$ tar zxf LATEST_VERSION.tar.gz
$ cd discov*
$ ./configure --prefix=/opt/biosoft/discovar && make -j 4 && make install
$ cd ..
$ rm -rf discov* LATEST_VERSION.tar.gz

2.2 DISCOVAR Denovo的下载和安装

此软件的安装需要GCC 4.7或以上版本,jemalloc 3.6.0或以上版本和samtools(如果使用bam文件,则需要)。

$ wget ftp://ftp.broadinstitute.org/pub/crd/DiscovarExp/LATEST_VERSION.tar.gz
$ tar zxf LATEST_VERSION.tar.gz
$ cd discov*
$ sudo yum install *malloc*
如果没有上一步,则在make过程中会提示错误“/usr/bin/ld: cannot find -ljemalloc”
$ ./configure --prefix=/opt/biosoft/discovarDenovo && make -j 4 && make install
$ echo 'export MALLOC_PER_THREAD=1' >> ~/.bashrc
上一步设置用于allowing per-threads memory management,能提高计算性能。
$ cd ..
$ rm -rf discov* LATEST_VERSION.tar.gz

3. 软件使用的注意事项

1. 强烈推荐使用 PCR-free protocol library 数据;数据量推荐为 ~60x,略大于或小于该值也是 OK 的。
2. 必须使用 Illumina MiSeq 或 HiSeq 2500 测序仪产生的 >=250 bp 长度的 Paired End 数据,并且首尾 reads 要有重叠。如果 PE 250bp 数据,则 Insert Size 长度要为 400-500 bp(需要注意的是软件的 manual 中可能写成 700bp,这是不对的)。
3. 只能使用一个文库的数据。,不支持输入 mate paired 数据。
4. DISCOVAR de novo (experimental) 能进行基因组的 de novo 组装,支持基因组大小可达 ~3 GB。

3. 软件的使用

3.1 DISCOVAR 的使用

软件的输入文件是 sort 过后的 Bam 文件,一个常用例子:

$ Discovar READS=sample-reads.bam REFERENCE=sample-genome.fasta              \
         REGIONS='10:30892106-30933760' OUT_HEAD=./discovar-variants/assembly\
         TMP=./discovar-variants/tmp

软件常用参数:

READS (String)
由逗号分割的一些 bam 文件,或内容为每行一个bam文件路径的 list 文件。
REGIONS (String)
对指定区域进行分析。多个区域则用逗号分割。区域的写法为 chr:start-sotp。如果 REGIONS=all,则对所有区域进行分析。
TMP (String)
指定临时文件路径
OUT_HEAD (String)
输出文件的前缀路径
NUM_THREADS (unsigned int) default: 0
使用的线程数。
REFERENCE (String)
参考序列 fasta 文件。若提供此文件,则能进行 variant calling,并给出 VCF 文件。

3.2 DISCOVAR de novo (experimental) 的使用

软件的输入文件是 sort 过后的 Bam 文件。程序在运行的时候会使用最大的线程数进行运算。

$ DiscovarExp --help special
上述命令用来查看软件的详细参数。
$ DiscovarExp READS=sample-reads.bam OUT_DIR=discovarexpOut
上述是软件的常用命令。同时,软件的参数非常少。
$ ls discovarexpOut/a.final/a.lines.fasta
查看主要结果。

4. DISCOVAR结果

4.1 结果表现形式

Understanding DISCOVAR output
图中,每个单独的箭头称为 edge,这些 edges 代表着序列;从起点到终点,有很多种不同的路径,称之为 lines;上图中有 4 个 cells,其中 3 个 cells 有 2 个 paths,有 1 个 cell 有 3 个 paths。
这种 multiple paths 可能表示:杂合位点;染色体变异;难以测序的位点等。

4.2 DISCOVAR 结果文件

生成的结果文件位于 discovar-variants/ 文件夹下,主要的结果文件是:

assembly.final.fasta 所有的 edges 序列 (edges overlap by K-1 bases)
assembly.final.fasta0 所有的 edges 序列 (without overlaps)
assembly.final.dot dot格式的组装图
assembly.final.variant VCF结果文件

4.3 DISCOVAR de novo 结果文件

生成的结果文件位于 discovarexpOut/a.final/ 文件夹下,主要结果文件有:

a.lines.fasta 多个 paths 中仅选择第一个 path,得到的 lines 序列的 fasta 文件。
a.lines.efasta 标准的 efasta 文件,有所有的 paths 结果。
a.fasta 所有的 edges 序列
a.lines 二进制文件
a.lines.src 上一个文件的文本形式结果

5. 总结

Discovar 能根据 Illumina 测序数据比对到基因组上的结果来进行基因组 de novo 组装,得到 edges 序列;若在提供了基因组序列的情况下,还能进行 Vaiants calling。
Discovar de novo (experimental) 能根据 Illumina 测序数据比对到基因组上的结果来进行基因组 de novo 组装,得到 edges 序列。相比与前者,还能得到 lines 序列,这是比较完整的序列文件。

RTL8188CE 802.11b无线网卡驱动的安装

lspci

03:00.0 Network controller: Realtek Semiconductor Co., Ltd. RTL8188CE 802.11b/g/n WiFi Adapter (rev 01)
$ lspci -nn|grep -i net
03:00.0 Network controller [0280]: Realtek Semiconductor Co., Ltd. RTL8188CE 802.11b/g/n WiFi Adapter [10ec:8176] (rev 01)
0c:00.0 Ethernet controller [0200]: Realtek Semiconductor Co., Ltd. RTL8111/8168/8411 PCI Express Gigabit Ethernet Controller [10ec:8168] (rev 07)

EL Repo

在网站http://elrepo.org/tiki/DeviceIDs上搜索10ec:8176

$ rpm --import https://www.elrepo.org/RPM-GPG-KEY-elrepo.org
$ rpm -Uvh http://www.elrepo.org/elrepo-release-6-6.el6.elrepo.noarch.rpm
$ yum search  kmod-r8192ce
$ yum install kmod-r8192ce

使用 RAxML 构建进化树

1. RAxML 简介

RAxML (Random Axelerated Maximum Likelikhood) 能使用多线程或并行化使用最大似然法构建进化树。
网页版工具:http://epa.h-its.org/raxml/submit_single_gene
参考文献:RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies

2. RAxML 下载与安装

$ wget https://github.com/stamatak/standard-RAxML/archive/v8.2.12.tar.gz -O ~/software/RAxML-v8.2.12.tar.gz
$ tar zxf ~/software/RAxML-v8.2.12.tar.gz -C /opt/biosoft/
$ mv /opt/biosoft/standard-RAxML-8.2.12/ /opt/biosoft/RAxML-8.2.12/
$ cd /opt/biosoft/RAxML-8.2.12/
$ make -f Makefile.SSE3.PTHREADS.gcc -j 4
$ rm *.o
$ make -f Makefile.AVX.PTHREADS.gcc -j 4
$ rm *.o
$ source ~/.bashrc.mpich
$ make -f Makefile.SSE3.HYBRID.gcc -j 4
$ rm *.o
$ make -f Makefile.AVX.HYBRID.gcc -j 4
$ rm *.o
$ chmod 755 /opt/biosoft/RAxML-8.2.12/usefulScripts/*
$ echo 'PATH=$PATH:/opt/biosoft/RAxML-8.2.12/' >> ~/.bashrc
$ source ~/.bashrc

2. RAxML 的使用

RaxML 软件包中带有一个 PDF 格式的 Manual 文档,介绍得非常详细。

2.1 RaxML 版本的选择

Sequential 版本适合于中小型的数据; PThreads 版本适合于长序列或多条序列;MPI 版本适合于较大(100~1000) bootstraps 的运行。

2.2 常用例子与参数

常用例子:

简单快速方式
$ raxmlHPC ­-f a ­-x 12345 ­-p 12345 ­-# 100 ­-m PROTGAMMALGX ­-s ex.phy ­-n ex -T 20

并行化软件支持,能最快速计算。并行化20个任务,每个任务使用8线程,能使用全部160线程计算资源:
$ /opt/sysoft/mpich2-1.5/bin/mpirun -np 20 raxmlHPC ­-f a ­-x 12345 ­-p 12345 ­-# 100 ­-m PROTGAMMALGX ­-s ex.phy ­-n ex -T 8

RAxML 的参数非常多,设置非常复杂,上述常用例子的参数为:

-f a
此参数用于选择 RAxML 运算的算法。可以设定的值非常之多。 a 表示执行快速 Bootstrap 分析并搜索最佳得分的 ML 树。
-x 12345
指定一个 int 数作为随机种子,以启用快速 Bootstrap 算法。
-p 12345
指定一个随机数作为 parsimony inferences 的种子。
-# 100
指定 bootstrap 的次数。
-m PROTGAMMALGX
指定核苷酸或氨基酸替代模型。PROTGAMMALGX 的解释: "PROT" 表示氨基酸替代模型; GAMMA 表示使用 GAMMA 模型; X 表示使用最大似然法估计碱基频率。
-s ex.phy
指定输入文件。phy 格式的多序列比对结果。软件包中包含一个程序来将 fasta 格式转换为 phy 格式。
-n ex
输出文件的后缀为 .ex 。
-T 20
指定多线程运行的 CPUs 。

2.3 结果文件

RAxML_bootstrap.ex           bootstrapped trees
RAxML_bestTree.ex            最佳得分 ML 树
RAxML_bipartitions.ex        有 bootstrap 分值支持的最佳得分树,分值在 node 上。
RAxML_bipartionsBranchLabels.ex 有 bootstrap 分值支持的最佳得分树, 分值在 branch 上。FigTree不能识别此文件。

使用 PhyML 构建进化树

1. PhyML 简介

使用 PhyML 构建最大似然树。
参考文献:New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0

2. PhyML 的下载和安装

$ wget http://www.atgc-montpellier.fr/download/binaries/phyml/PhyML-3.1.zip
$ unzip PhyML-3.1.zip
$ mv PhyML-3.1 /opt/biosoft/
$ ln -s /opt/biosoft/PhyML-3.1/PhyML-3.1_linux64 /opt/biosoft/PhyML-3.1/PhyML
$ echo 'PATH=$PATH:/opt/biosoft/PhyML-3.1/' >> ~/.
$ source ~/.bashrc

3. PhyML 的使用

PhyML 的输入文件为 phylip 格式。
常用例子:

$ PhyML -i proteins.phy -d aa -b 1000 -m LG -f m -v e -a e -o tlr

常用参数:

-i seq_file_name
输入文件,phylip 格式的多序列比对结果。
-d data_type default:nt
该参数的值为 nt, aa 或 generic。
-b int
设置 bootstrap 次数。
-m model
设置替代模型。 核酸的模型有: HKY85(默认的), JC69, K80, F81, TN93, GTR ; 氨基酸的模型有:LG (默认的), WAG, JTT, MtREV, Dayhoff, DCMut, RtREV, CpREV, VT, Blosum62, MtMam, HIVw, HIVb 。
-f e,m or fA,fC,fG,fT
设置频率计算的方法。 e 表示使用比对结果中不同氨基酸或碱基出现的频率来计算; m 表示使用最大似然法计算碱基频率,或使用替换模型计算氨基酸频率; fA,fC,fG,fT 则是 4 个浮点数,表示 4 中碱基的频率,仅适合核酸序列。
-v prop_invar
设置不变位点的比例,是一个[0,1]区间的值。或者使用 e 表示程序获得其最大似然估计值。
-a gamma
gamma 分布的参数。此参数值是个正数,或者使用 e 表示程序获得其最大似然估计值。在 ProtTest 软件给出的最优模型中含有 G 时,使用该参数。
-o params
参数优化的选项。t 表示对 tree topology 进行优化; l 表示对 branch length 进行优化; r 表示对 rate parameters 优化。
params=tlr 这表示对 3 者都进行优化。 params=n 表示不进行优化。

4. PhyML 结果

PhyML 的输出结果为:

proteins.phy_phyml_tree.txt        :    最大似然法构建的进化树
proteins.phy_phyml_boot_stats.txt  :    bootstrap 的统计信息
proteins.phy_phyml_boot_trees.txt  :    bootstrap 树
proteins.phy_phyml_stats.txt       :    程序运行的中的参数和结果统计

使用 ProtTest 来选择最优氨基酸替代模型

1. ProtTest 简介

ProtTest 用来进行最优氨基酸替代模型的选择。相应的,适用于核苷酸的软件是 jModeltest。
ProtTest 通过 PhyML 对进化树和模型参数的最大似然估计,通过 AIC, BIC 分值或 DT 来寻找最佳模型。分值越小越优。
ProTest 3.2 版本包含 15 种不同类型的 rate matrices;考虑到位点的 rate variation (+I: invariable sites; +G: gamma-distributed rates) 和 observed amino acid frequencies (+F), 共有 120 种不同的模型。
ProtTest 官网:https://code.google.com/p/prottest3/
从此处下载该软件。可能需要设置代理后下载。
参考文献:ProtTest 3: fast selection of best-fit models of protein evolution

2. ProtTest 下载和安装

$ tar zxf prottest-3.4-20140123.tar.gz -C /opt/biosoft
$ cd /opt/biosoft/prottest-3.4-20140123
$ echo 'export PROTTEST_HOME=/opt/biosoft/prottest-3.4-20140123' >> ~/.bashrc
查看说明文档:
$ less README

3. ProtTest 的使用

ProtTest 使用 JAVA 编写,有图形化和命令行两种运行模式。

3.1 图形化界面使用

必须要进入到程序的所在的目录运行程序以启动图形化界面
$ cd /opt/biosoft/prottest-3.4-20140123/runXProtTestHPC.sh
$ runXProtTestHPC.sh

启动 JAVA 界面后,点击 File–Load Alignment, 上传多序列比对结果;然后点击 Analysis–Compute likehood scores, 选择所使用的线程数,以及候选模型的选择,和计算 likelihood 的 topology;然后点击 Compute, 进行计算,所需要消耗的实际有点长;计算完毕后,点击 Selection–Results 来查看结果。通过 AIC, BIC, AICc 和 DT 来查看其得分,点击表格的第1行进行排序,寻找分值最小的模型作为最优氨基酸替代模型。

3.2 命令行运行

常用例子:

不加参数运行,则给出帮助文档:
java -jar /opt/biosoft/prottest-3.4-20140123/prottest-3.4.jar

常用的命令行:
java -jar /opt/biosoft/prottest-3.4-20140123/prottest-3.4.jar -i proteins.phy -all-distributions -F -AIC -BIC -tc 0.5 -threads 24 -o prottest.out

ProtTest 的常用参数:

-i alignment_filename
必须参数,输入多序列比对结果文件。
-o output_filename
输出的文件名。不设置,则默认输出到标准输出。
-[matrix]
指定需要分析的 matrix 。 该 matrix 可以被替换为 JTT LG DCMut MtREV MtMam MtArt Dayhoff WAG RtREV CpREV Blosum62 VT HIVb HIVw FLU 这 15 种 matrix。 若不指定,则默认全选。
-all-distributions
指定 matrix 模型结合 I 或 G 或 I+G
-F
指定 matrix 模型结合 empirical frenquency estimation
-AIC
输出结果中按 AIC (Akaike Information Criterion) 排序
-BIC
输出结果中按 BIC (Bayesian Information Criterion) 排序
-AICC
输出结果中按 AICc (Corrected Akaike Information Criterion) 排序
-DT
输出结果中按 DT (Decision Theory Criterion) 排序
-tc consensus_threshold
输出满足指定阈值的一致树。该值在 0.5~ 1.0 之间。[0.5 = majority rule consensus ; 1.0 = strict consensus]
-threads number_of_threads
使用的 CPU 数。

使用 Gblocks 提取保守序列

1. Gblocks 简介

Gblocks用于从多序列比对结果中提取保守位点,以利于下一步的进化分析。
在线说明文档:http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_documentation.html
在线服务网站:http://molevol.cmima.csic.es/castresana/Gblocks_server.html

2. Gblocks 安装

$ wget http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_Linux64_0.91b.tar.Z
$ sudo yum install -y ncompress
$ tar Zxf Gblocks_Linux64_0.91b.tar.Z -C /opt/biosoft/
$ echo 'PATH=$PATH:/opt/biosoft/Gblocks_0.91b/' >> ~/.bashrc
$ source ~/.bashrc

3. Gblocks 使用

Gblosk 有两种使用方式,第一种是交互式的方式(按提示输入文件改变参数),第二种是命令行式(在命令行中输入参数)。
命令行式的常用例子:

使用默认的设置:
$ Gblocks proteins.fasta -t=p
必须是 fasta 文件在前,参数在后。若没有参数,则进入交互式界面。

得到更长的序列:
$ Gblocks proteins.fasta -b4=5 -b5=h

命令行式的常用参数:

-t= default:p
设置序列的类型,可选的值是 p,d,c 分别代表 protein, DNA, Codons 。
-b1= default:( 序列条数的 50% + 1 )
设定保守性位点必须有 >= 该值的序列数。该参数后接一个 integer 数,默认下比序列条数的 50% 大 1.
-b2= default: 序列条数的 85%
确定保守位点的侧翼位点时,其位点必须有 >= 该值的序列数。该值必须要比 -b1 的值要大。
-b3= default: 8
最大连续非保守位点的长度。
-b4= default: 10
保守位点区块的最小长度。该值必须 >=2 。
-b5= default: n
设置允许含有 Gap 位点。可选的值有 n,h,a 分别代表 None, With Half, All 。 当为 h 时,表示
-e= default: -gb
设置输出结果的后缀。

使用 MAFFT 进行多序列比对

1. MAFFT 简介

最经典和广为熟知的多序列比对软件是 clustalw 。 但是现有的多序列比对软件较多,有文献报道:比对速度(Muscle>MAFFT>ClustalW>T-Coffee),比对准确性(MAFFT>Muscle>T-Coffee>ClustalW)。因此,推荐使用 MAFFT 软件进行多序列比对。

2. MAFFT 下载安装

$ wget http://mafft.cbrc.jp/alignment/software/mafft-7.158-without-extensions-src.tgz
$ tar zxf mafft-7.158-without-extensions-src.tgz
$ cd mafft-7.158-without-extensions/core
$ perl -p -i -e 's#PREFIX =.*#PREFIX = /opt/biosoft/mafft#' Makefile
$ perl -p -i -e 's#BINDIR =.*#BINDIR = /opt/biosoft/mafft/bin/#' Makefile
$ make
$ make install
$ echo 'PATH=$PATH:/opt/biosoft/mafft/bin/' >> ~/.bashrc
$ source ~/.bashrc

检测软件是否正确安装
$ cd ../test
$ rehash                                                   # if necessary
$ mafft sample > test.fftns2                               # FFT-NS-2
$ mafft --maxiterate 100  sample > test.fftnsi             # FFT-NS-i
$ mafft --globalpair sample > test.gins1                   # G-INS-1
$ mafft --globalpair --maxiterate 100  sample > test.ginsi # G-INS-i
$ mafft --localpair sample > test.lins1                    # L-INS-1
$ mafft --localpair --maxiterate 100  sample > test.linsi  # L-INS-i
$ diff test.fftns2 sample.fftns2
$ diff test.fftnsi sample.fftnsi
$ diff test.gins1 sample.gins1
$ diff test.ginsi sample.ginsi
$ diff test.lins1 sample.lins1
若 diff 的结果不换回异常,则正确安装。

3. MAFFT 使用

MAFFT 有一些参数,适合不同情况下的多序列比对。
软件输入为 FASTA 格式,默认输出到标准输出。

3.1 较为精确的方法

L-INS-i

最准确的方法。适合于 <200 条序列,且序列长度 <~2000 aa/nt 的比对。

$ mafft –localpair –maxiterate 1000 input [> output]
$ linsi input [> output]

G-INS-i

适合于序列长度相似的多序列比对。序列条数 <200, 序列长度 <~2000 aa/nt 。

$ mafft –globalpair –maxiterate 1000 input [> output]
$ ginsi input [> output]

E-INS-i

适合序列中包含较大的非匹配区域。序列条数 <200, 序列长度 <~2000 aa/nt 。

$ mafft –ep 0 –genafpair –maxiterate 1000 input [> output]
$ einsi input [> output]

3.2 节约速度的方法

FFT-NS-i

减少迭代次数,最大迭代次数减为 2 。

$ mafft --retree 2 --maxiterate 2 input [> output]
$ fftnsi input [> output]

FFT-NS-2

最大迭代次数减为 0 。

$ mafft --retree 2 --maxiterate 0 input [> output]
$ fftns input [> output]

FFT-NS-1

此方法非常快速,适合 >2000 条序列的多序列比对。

$ mafft --retree 1 --maxiterate 0 input [> output]

NW-NS-i

迭代过程中不进行 FFT aproximation

$ mafft --retree 2 --maxiterate 2 --nofft input [> output]
$ nwnsi input [> output]

NW-NS-2

$ mafft --retree 2 --maxiterate 0 --nofft input [> output]
$ nwns input [> output]

NW-NS-PartTree-1

3 个参数都设置为最不消耗时间的类型,适合于 ~10,000 到 ~50,000 条序列的比对。

$ mafft --retree 1 --maxiterate 0 --nofft --parttree input [> output]

antiSMASH 使用

1. antiSMASH 简介

antiSMASH 用于寻找次级代谢基因簇。一般情况下,参与次级代谢途径中生物合成酶的基因在染色体上成簇排列。基于指定类型的 profile hidden Markov models, antiSMASH 能准确鉴定所有已知的次级代谢基因簇。
antiSMASH 的使用说明:http://antismash.secondarymetabolites.org/help.html
antiSMASH 的参考文献:antiSMASH 2.0 — a versatile platform for genome mining of secondary metabolite producers.
Kai Blin, Marnix H. Medema, Daniyal Kazempour, Michael A. Fischbach, Rainer Breitling, Eriko Takano, & Tilmann Weber
Nucleic Acids Research (2013), doi: 10.1093/nar/gkt449.

2. 次级代谢基因簇简介

在 antiSMASH 中,将次级代谢基因簇分为了 24 类。
最常见的次级代谢基因簇是 type I, II and III polyketides synthase(PKS) 和 non-ribosomal peptides synthase(NRPS)。例如:四环素、大环内酯类、安莎类、聚醚类由 PKS 途径合成;beta-内酰胺类、多肽类、糖肽类由 NRPS 途径合成。
此外,还有 heterocyst glycolipid-like polyketides, terpenes, lantibiotics, bacteriocins, β-lactams, aminoglycosides/aminocyclitols, aminocoumarins, siderophores, ectoines, butyrolactones, indoles, nucleosides, phosphoglycolipids, melanins, oligosaccharide, furans, homoserine lactones, thiopeptides, phenazines, others.

3. antiSMASH 安装

首先要安装 ncbi-blast+, hmmer3.0, hmmer2.3.2(hmmpfam), glimmer3, GlimmerHMM 3.0.2, muscle.
安装 ncbi-blast+

$ wget ftp://ftp.ncbi.nih.gov/blast/executables/LATEST/ncbi-blast-2.2.29+-x64-linux.tar.gz
$ tar zxf ncbi-blast-2.2.29+-x64-linux.tar.gz -C /opt/biosoft/
$ echo 'PATH=$PATH:/opt/biosoft/ncbi-blast-2.2.29+/bin/' >> ~/.bashrc

安装 hmmer3.0

$ wget http://selab.janelia.org/software/hmmer3/3.0/hmmer-3.0.tar.gz
$ tar zxf hmmer-3.0.tar.gz
$ cd hmmer-3.0
$ ./configure --prefix=/opt/biosoft/hmmer-3.0/ && make -j 4 && make install
$ echo 'PATH=$PATH:/opt/biosoft/hmmer-3.0/bin/' >> ~/.bashrc

安装 hmmer2.3.2

$ wget http://selab.janelia.org/software/hmmer/2.3.2/hmmer-2.3.2.tar.gz
$ tar zxf hmmer-2.3.2.tar.gz
$ cd hmmer-2.3.2
$ ./configure --prefix=/opt/biosoft/hmmer-2.3.2/ && make -j 4 && make check
$ sed -e "s#\(cp src/\$\$file \$(BINDIR)/\);#\1\$\${file}2;#" -i Makefile
$ make install
$ echo 'PATH=$PATH:/opt/biosoft/hmmer-2.3.2/bin/' >> ~/.bashrc

安装 glimmer3

$ wget http://ccb.jhu.edu/software/glimmer/glimmer302.tar.gz
$ tar zxf glimmer302.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/glimmer3.02
$ wget https://bitbucket.org/antismash/antismash2/downloads/Allow-glimmer-to-compile-on-g-4.4.3.patch -O Allow-glimmer-to-compile-on-g-4.4.3.patch
$ patch -p1 < Allow-glimmer-to-compile-on-g-4.4.3.patch
$ cd src
$ make -j 4
$ echo 'PATH=$PATH:/opt/biosoft/glimmer3.02/bin/' >> ~/.bashrc

安装 GlimmerHMM 3.02

$ wget ftp://ccb.jhu.edu/pub/software/glimmerhmm/GlimmerHMM-3.0.2.tar.gz
$ tar zxf GlimmerHMM-3.0.2.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/GlimmerHMM/sources
$ make
$ cp glimmerhmm ../bin/
$ echo 'PATH=$PATH:/opt/biosoft/GlimmerHMM/bin/' >> ~/.bashrc

安装 Muscle

$ wget http://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_src.tar.gz
$ tar zxf muscle3.8.31_src.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/muscle3.8.31/src
$ make -j 4
$ mkdir ../bin/
$ cp muscle ../bin/
$ echo 'PATH=$PATH:/opt/biosoft/muscle3.8.31/bin/' >> ~/.bashrc
$ source ~/.bashrc

安装一些系统软件:

$ sudo yum install -y perl-Archive-Tar python-pip python-virtualenv git java-1.7.0-openjdk python-devel libxslt-devel libxml2-devel gcc-c++ patch glibc-static cairo
$ wget http://dl.fedoraproject.org/pub/epel/6/x86_64/epel-release-6-8.noarch.rpm
$ sudo rpm -Uvh epel-release-6-8.noarch.rpm

安装 antiSMASH。可以不需要以上所有步骤,直接进行下面的安装(需要联网)。以上步骤则是将所有的程序安装到 /opt/biosoft 目录下,以便于管理生物信息学软件。

$ mkdir /opt/biosoft/antiSMASH/
$ cd /opt/biosoft/antiSMASH/
$ wget https://bitbucket.org/antismash/antismash2/downloads/install_centos.sh -O install_centos.sh
此文件中有处错误,导致下载不了 antiSMASH 的软件包。
$ perl -p -i -e 's#\"\${ANTISMASH_BASE}\${ANTISMASH_TARBALL}\"#\${ANTISMASH_BASE}/\${ANTISMASH_TARBALL} -O \${ANTISMASH_BASE}#' install_centos.sh
此脚本适合于 CentOS 6.4 系统,如果是 CentOS 6.5 系统, 则要多进行下一步
$ perl -p -i -e 's/6_4/6_5/' install_centos.sh
$ sh install_centos.sh
$ echo 'PATH=$PATH:/opt/biosoft/antiSMASH/' >> ~/.bashrc
$ source ~/.bashrc
$ run_antismash

4. antiSMASH 使用

4.1 注意事项

antiSMASH 支持 Fasta/Genbank/EMBL (要分别以 .fasta .gbk .embl 作为后缀以利于程序识别)格式的文件作为输入。推荐使用 Genbank 格式文件作为输入。该文件包含了编码蛋白基因的注释信息。否则,以 fasta 文件作为输入,程序则需要调用 Glimmer3 和 GlimmerHMM 来进行基因预测后再进行次级代谢基因簇的鉴定。
使用 –clusterblast 和 –subclusterblast 参数,antiSMASH 使用 blastp 来将氨基酸序列比对到已知的次级代谢 clusters 或 subclusters 上,来寻找 query 序列中的基因簇。网页版中默认使用此参数。
使用 –smcogs 参数,antiSMASH 能分析次级代谢基因家族 (smCOGs),并使用其家族的基因(最多100个)构建系统发育树。网页版中默认使用此参数。
使用 –full-hmmer 参数, antiSMASH 将进行全基因组的 PFAM 分析,寻找次级代谢 domains 出现过于频繁的基因组区域。这样能找到一些 clusterblast 步骤中漏掉的基因簇。网页版中默认使用此参数。

4.2 常用例子

$ run_antismash --clusterblast --subclusterblast --smcogs --full-hmmer species.gbk

结果文件生成于 species 文件夹下。点击 index.html 进行结果的网页查看。
有关于基因组 genbank 文件的生成,可以使用 tbl2asn 软件进行生成。