使用 SSPACE 进行 scaffoding

SSPACE 能利用 paired reads 的比对结果,将 contigs 或 scaffolds 连接成 scaffolds。其参考文献:Boetzer M, Henkel C V, Jansen H J, et al. Scaffolding pre-assembled contigs using SSPACE[J]. Bioinformatics, 2011, 27(4): 578-579.

1. 安装 SSPACE

软件下载页面:http://www.baseclear.com/lab-products/bioinformatics-tools/sspace-standard/

$ tar zxf SSPACE-STANDARD-3.0_linux-x86_64.tar.gz
$ ./SSPACE-STANDARD-3.0_linux-x86_64/SSPACE_Standard_v3.0.pl

解压缩软件包后,运行软件文件夹中的 perl 程序即可运行 SSPACE。软件主目录下包含一些软件使用说明和示例等,其中 README 文件描述得非常详细。

2. SSPACE 使用方法

2.1 library 文件

首先要建立一个描述 library 信息的文本文件,例如:

Lib1 bwa file1.1.fasta file1.2.fasta 400 0.25 FR
Lib1 bowtie file2.1.fasta file2.2.fasta 400 0.25 FR
Lib2 bwasw file3.1.fastq file3.2.fastq 4000 0.5 RF
Lib2 TAB file4.tab 4000 0.5 RF
Lib3 TAB file5.tab 10000 0.5 RF
unpaired bowtie unpaired_reads1.fasta
unpaired bwasw unpaired_longreads1.gz

此 library 文件由多列组成,列与列之间由 1 个 空格 或 tab 分隔,各列意义如下:

第 1 列: library 名称。程序运行过程中产生的临时文件以此来命名; 多个行可以拥有同一个 library 名称,则其具有相同的 library 设置和不同的数据文件; 同时,libraries 必须按 insert size 来排序,inert size 最小的必须放到第一行,这是因为进行 scaffold 构建时,按此文件提供的 libraries 的顺序来输入数据的; unpaired reads, 则第一列是 ‘unpaired’。
第 2 列: 将 reads 比对到基因组上所使用的软件名, 可以为 bowtie 、 bwa 和 bwasw 等; 如果输入的数据是 reads 比对过后的 tab 格式结果,则此列为 “TAB”。
第 3,4 列: Fasta 或 Fastq 格式的双末端测序文件,并且文件中成对的 paired reads 必须在两个文件中并处于相同的行号上,同时,软件读取数据与序列的 headers 无关。如果是 unpaired reads,则仅需要第 3 列,为 tab 格式的 reads mapping 结果,过后详述。
第 5,6 列:第 5 列为 insert size 的期望值; 第 6 列为 insert size 允许的最小偏差。 比如,这两列值分别为 4000 和 0.5,则 insert size 在 2000-6000 之间的 pairs 才是有效 pairs。
第 7 列:paired-reads 的方向,有 FF,FR,RF 或 RR 几种选项。

2.2 程序参数

-l 输入的 library 文件
-s 输入的 Fasta 文件
-x 是否对 contigs 进行延长。其值可以为 0 或 1。 1 表示进行延伸,0 表示不延伸。默认值为 0。

延伸参数:
-m 进行延伸时,read 和基因组序列最小的 overlap。此值越大,则结果越准确,同时耗内存越少。推荐此值接近最长的 read 的长度。比如,对于 26 bp 长度的 reads, 该值适合设为 32~35。 默认此值为 32 。此值取值范围为 15~50 。软件运行时,将 unmapped reads 全部打断成 m+1 长度的序列,这些序列用于进行 contigs 的延伸。
-o 进行延伸时,延伸 1 个碱基需要的最小 reads 数。此值越大,则结果越准确。默认值为 20 。
-r 进行延伸时,延伸 1 个碱基,此碱基在所有匹配的 reads 中的最小比例。此值越大,则结果越准确。默认值为 0.9 。

Scaffolding 参数:
-k 将两个 contigs 连接成 scaffold 时,需要的最小的 reads pairs 数目。默认值为 5 。
-a 将两个 contigs 连接成 scaffold 时,这两个 contigs 之间的连接数 与 其和其它 contigs 的连接数之间的最小比值。此值越大,则结果越准确。默认值为 0.70
-n 在 scaffold 中,将两个邻近的 contigs 合并到一起需要的最小的 overlap。默认值为 15。
-z 进行 scaffolding 时,允许的最小的 contig 长度。低于此长度的 contig 将不能用于进行 scaffold 组装。默认值为 0 。较长的 contigs 产生的 scaffolds 比较可信; 而小于 100bp 的 contigs 容易是重复序列。

bowtie 比对参数:
-g 使用 bowtie 进行比对时,允许的最大 gaps 数。默认值为 0

其它参数:
-T 设定运行的线程数。默认值为 1。
-b 输出文件夹名及文件夹内的文件前缀。
-S 当程序正在运行时,跳过读取 reads 的阶段。和 -b 参数结合使用,则可以同时运行多个 SSPACE 程序,对每个程序设置不同的参数,这样能较快得到较好的结果。
-v verbose mode
-p 生成可供可视化的 .dot 文件。

2.3 其它工具

SSPACE 提供了一些其它比较有用的小工具:

estimate_insert_size.pl 用于计算 insert size。此程序计算的结果有些问题。
fastq_qualitytrim_pairs.pl 对 reads pairs 进行质量控制的程序。

sam_bam2tab.pl 将 bam sam 文件转换为 tab 格式的程序。

GAGE-B 进行微生物基因组组装软件的评估

参考文献:Magoc T, Pabinger S, Canzar S, et al. GAGE-B: an evaluation of genome assemblers for bacterial organisms[J]. Bioinformatics, 2013, 29(14): 1718-1725.

评价了 8 个基因组组装软件: ABySS CABOG MIRA MaSuRCA SGA SOAPdenovo SPAdes Velvet。
若这 8 个软件组装效果全部一致,则值全部为 1 。

对 7 个物种的组装结果, 根据 N50, contig 的组装性能:
ABySS 0.79
CABOG 0.55
MIRA 0.85
MaSuRCA 2.08
SGA 0.24
SOAPdenovo 1.10
SPAdes 1.91
Velvet 0.48

对 5 个物种的组装结果, 根据 N50, scaffold 的组装性能:
ABySS 0.76
CABOG 0.81
MIRA 0.85
MaSuRCA 1.53
SGA 0.31
SOAPdenovo 0.81
SPAdes 1.20
Velvet 0.72

可以看出,从基因组的连续性考虑,最佳的细菌基因组组装软件是:
MaSuRCA > SPAdes 。仅有这两个软件能高于平均水平, 特别是 MaSuRCA 显著比其它软件要好。

文中没有进行组装结果的准确性分析。但是提到 SOAPdenvo 可以产生最大的 scaffolds,但是比其它软件的结果有根多的 local errors。

由于数据问题,没有对 ALLPATHS-LG 进行分析。

使用 SPAdes 进行基因组组装

1. SPAdes 简介

SPAdes 主要用于进行单细胞测序的细菌基因组组装,当然也能用于非单细胞测序数据。输入数据可以是 Illumina、IonTorrent reads,或 PacBio、Sanger reads,也可以把一些 contigs 序列作为 long reads 进行输入。该软件可以同时接受多组 paired-end、mate-pairs 和 unpaired reads 数据的输入。同时该软件有一个独立的模块用于进行杂合基因组的组装。
文献:Bankevich A, Nurk S, Antipov D, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[J]. Journal of Computational Biology, 2012, 19(5): 455-477.
软件说明文档:http://spades.bioinf.spbau.ru/release3.0.0/manual.html#sec2.1
SPAdes 包含多个独立模块:

BayesHammer: 用于 Illumina reads 的修正
IonHammer: 用于 IonTorrent 数据的修正
SPAdes: 用于基因组组装;K 值是软件自动选择的。
MismatchCorrector: 对组装结果 contigs 或 scaffolds 的 mismatch 和 sort indels 进行修正;此模块需要使用 BWA;默认下此模块是关闭的,但是推荐开启它。
dipSPAdes: 用于组装双倍型高杂合基因组;更多说明:dipSPAdes manual

2. SPAdes 下载和安装

从网页http://bioinf.spbau.ru/spades下载需要填写一些个人资料。

$ wget http://spades.bioinf.spbau.ru/release3.1.0/SPAdes-3.1.0-Linux.tar.gz
$ tar -xzf SPAdes-3.1.0-Linux.tar.gz -C /opt/biosoft/
$ /opt/biosoft/SPAdes-3.1.0-Linux/bin/spades.py --test
测试运行

3. SPAdes 的使用

3.1 输入数据

SPAdes 3.0.0 的运行至少需要有以下一种数据:

Illumina paired-end/unpaired reads
IonTorrent paired-end/unpaired reads
PacBio CCS reads

并且,值得注意的是:Illumina 数据和 Ionorrent 数据不能同时用于组装; 仅有 mate-paired,PacBio CLR reads, Sanger reads 或 additional contigs 数据时,不应该使用 SPAdes 进行组装。
SPAdes 支持的 Paired-end 和 Mate-Paired 的数据,其数据需要为 fastq 格式,软件需要对其进行 reads 的 error correction 。同时, SPAdes 也支持使用 Sanger 或 PacBio CCS 的 reads 数据,但软件不能对此数据进行 error correction。

3.1.1 Read-pair 数据

Read-pair 数据输入到程序中有 3 种方式:
1. left 和 right 的 reads 分别在两个 fastq 文件中。
2. left 和 right 的 reads 交叉融合在一个 fastq 文件中。
3. 将所有的输入数据信息整合在一个 YAML 格式的文本文件中。

使用非 YAML 方式输入数据,这种方式最多能使用 5 组 paired-end 数据 和 5 组 mate-paired 数据。
仅有一个 library 数据时:

--12 file_name
12 表示后面接的文件是交叉融合的 paired 数据,下同。
-1 file_name
1 表示 left 数据
-2 file_name
2 表示 right 数据
-s file_name
s 表示 single 数据, 也用于输入 PacBio CCS reads。

有多个 paired-end library 数据时:

--pe{int}-12 编号为 int 的 library 的交叉融合后的 paired 数据。int 取值只能是 1,2,3,4,5 。下同。
--pe{int}-1  编号为 int 的 library 的 left 数据
--pe{int}-2  编号为 int 的 library 的 right 数据
--pe{int}-s  编号为 int 的 library 的 single 数据
--pe{int}-{fr|rf|ff} 编号为 int 的 library 的数据的方向,默认为 --pe{int}-fr 。

有多个 mate-paired library 数据时:

--mp{int}-12 编号为 int 的 library 的交叉融合后的 paired 数据
--mp{int}-1  编号为 int 的 library 的 left 数据
--mp{int}-2  编号为 int 的 library 的 right 数据
--mp{int}-s  编号为 int 的 library 的 single 数据
--mp{int}-{fr|rf|ff} 编号为 int 的 library 的数据的方向,默认为 --mp{int}-rf 。

3.1.2 PacBio 数据

PacBio 数据有两种: CCS (circular consensus sequence) 和 CLR (continuous long read)。PacBio CLR 数据有利于杂合基因组的组装。
数据输入参数:

--pacbio  输入 PacBio CLR reads
--sanger  输入 sanger reads

3.1.3 Additional contigs

--trusted-contigs
输入可信度高的 contigs,用于 graph construction, gap closure 和 repeat resolution。
--untrusted-contigs
输入可信度较低的 contigs, 用于gap closure 和 repeat resolution。

这两个参数不能输入其它邻近物种的基因组序列。仅用于输入同一个物种的基因组的 contigs 。

3.1.4 YAML 方式输入数据

--dataset
YAML 格式的文件
    [
      {
        orientation: "fr",
        type: "paired-end",
        right reads: [
          "/FULL_PATH_TO_DATASET/lib_pe1_right_1.fastq",
          "/FULL_PATH_TO_DATASET/lib_pe1_right_2.fastq"
        ],
        left reads: [
          "/FULL_PATH_TO_DATASET/lib_pe1_left_1.fastq",
          "/FULL_PATH_TO_DATASET/lib_pe1_left_2.fastq"
        ]
      },
      {
        orientation: "rf",
        type: "mate-pairs",
        right reads: [
          "/FULL_PATH_TO_DATASET/lib_mp1_right.fastq"
        ],
        left reads: [
          "/FULL_PATH_TO_DATASET/lib_mp1_left.fastq"
        ]
      },
      {
        type: "single",
        single reads: [
          "/FULL_PATH_TO_DATASET/pacbio_ccs.fastq"
        ]
      },
      {
        type: "pacbio",
        single reads: [
          "/FULL_PATH_TO_DATASET/pacbio_clr.fastq"
        ]
      }
    ]

3.2 参数

使用 spades.py 运行 SPAdes 程序:

$ /opt/biosoft/SPAdes-3.0.0-Linux/bin/spades.py [options] -o output_dir

参数:

-o output_dir
指定输出的文件夹
--sc
此 flag 用于 MDA (single-cell) 数据
--iontorrent
此 flag 用于 IonTorrent 数据的组装
--test
使用 test 数据运行 SPAdes,用于检测软件是否正确安装
-h | --help
打印帮助信息
--only-error-correction
仅仅执行 reads error correction 步骤
--only-assembler
仅仅运行组装模块
--careful
通过运行 MismatchCorrector 模块进行基因组上 mismatches 和 short indels 的修正。推荐使用此参数。
--continue
从上一次终止处继续运行程序。
--restart-from
从指定的位置重新开始运行程序。和上一个参数相比,此参数可以用于改变一些组装参数。可选的值有:

ec 从 error correction 处开始
as 从 assembly module 处开始
k{int} 从指定的 k 值处开始
mc 从 mismatch correction 处开始

--disable-gzip-output
使用此参数来设定不对 corrected reads 进行压缩。默认下 corrected reads 是 .fastq.gz 格式的
-t int
使用的线程数,默认为16
-m int
设定内存的限制,单位为 Gb。如果程序使用的内存达到此值,则程序会终止运行。默认值是 250 。
--tmp-dir dir_name
设置 reads error correction 的临时文件存放路径。默认为 output_dir/corrected/tmp 。
-k int,int,...
由逗号分隔的 k-mer sizes。这些数值必须为奇数,要小于 128,且按升序排列。如果使用了 --sc 参数,则默认值为 21,33,55 。 若没有 --sc 参数,则程序会根据 reads 长度自动选择 k-mer 参数。
--phred-offset
碱基质量格式, 33 或 64

3.3 常用例子

单个 illumina paired-end 文库:

$ spades.py -o output_dir -1 reads1.fastq -2 reads2.fastq

多个 illumina paired-end 和 mate-paired 文库,以及 Pcabio sanger contigs 数据:

$ spades.py -o output_dir\
 --pe1-1 pe1_1.fq --pe1-2 pe1_2.fq --pe2-1 pe2_1.fq --pe2-2 pe2_2.fq\
 --mp1-1 mp1_1.fq --mp1-2 mp1_2.fq --mp2-1 mp2_1.fq --mp2-2 mp2_2.fq\
 -s pacibo_ccs.fastq --pacbio pacbio_clr.fastq\
 --sanger sanger.fa\
 --trusted-contigs trusted_contig.fa\
 --untrusted-contigs untrusted_contig.fa\
 --careful -t 16 --phred-offset 33 -m 250 -k 21,33,55 [--sc]

使用 PASHA 进行基因组组装

1. PASHA 简介

并行化运算,速度快的基因组组装软件。

2. PASHA 的下载和安装

$ wget http://cznic.dl.sourceforge.net/project/pasha/Pasha-1.0.10.tar.bz2
$ tar jxf Pasha-1.0.10.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/Pasha-1.0.10/
安装 PASHA 前需要安装 TBB
$ wget www.threadingbuildingblocks.org/sites/default/files/software_releases/linux/tbb42_20140416oss_lin.tgz
$ tar zxf tbb42_20140416oss_lin.tgz
修改 Pasha 软件的 Makefile 文件中 TBB 文件夹路径的指向
$ perl -p -i -e 's#~/install#/opt/biosoft/Pasha-1.0.10/tbb42_20140416oss#' Makefile
$ make
$ echo 'PATH=$PATH:/opt/biosoft/Pasha-1.0.10/bin' >> ~/.bashrc
$ source ~/.bashrc

2. 程序的参数

2.1 pasha-kmergen

此命令用于生成 k-mers。输入文件可以为 fasta 或 fastq 文件。
参数:

-fasta str
输入 fasta 格式的数据
-fastq str
输入 fastq 格式的数据
-k int
输入 k-mer size,从 1 到 31 的奇数,默认值为 21

常用例子:

$ pasha-kmergen DIRECTORY -fasta in.fasta -k KMER_SIZE
$ pasha-kmergen DIRECTORY -fasta in.fasta -fasta in2.fasta -k KMER_SIZE
$ pasha-kmergen DIRECTORY -fasta in.fasta -fastq in2.fastq -k KMER_SIZE
$ mpirun -np 8 pasha-kmergen DIRECTORY -fasta in.fasta -k KMER_SIZE

每个 pasha-kmergen 使用 2 个线程。因此,推荐 -np 对应的值最好低于 CPU cores 的一半。

2.2 pasha-pregraph

此命令用于利用 k-mers 来构建初步的 de Bruijn 图,再去除 tips 和 low-coverage paths。
参数:

-fasta str
输入 fasta 格式的数据
-fastq str
输入 fastq 格式的数据

常用例子:

$ pasha-pregraph DIRECTORY -fasta in.fasta
$ pasha-pregraph DIRECTORY -fasta in.fasta -fasta in2.fasta
$ pasha-pregraph DIRECTORY -fasta in.fasta -fastq in2.fastq
$ mpirun -np 8 pasha-pregraph DIRECTORY -fasta in.fasta

当使用 MPI 进行并行运算是,pasha-pregraph 和 pasha-kmergen 使用 -np 指定的线程数要一致。

2.3 pasha-graph

此命令用于融合 bubbles,生成 contigs,并进行 scaffolding。
参数:

-fasta str
输入 fasta 格式的数据。程序识别此数据为单端数据。
-fastq str
输入 fastq 格式的数据。程序识别此数据为单端数据。
-fastaPaired str
输入 fasta 格式的数据。后接两个用空格分隔的文件。程序识别此数据为双端数据。
-fastqPaired str
输入 fastq 格式的数据。后接两个用空格分隔的文件。程序识别此数据为双端数据。
-fastaPairedFile str
输入 fasta 格式的数据。后接 1 个文件,此文件是双端数据交叉整合的结果。程序识别此数据为双端数据。
-fastqPairedFile str
输入 fastq 格式的数据。后接 1 个文件,此文件是双端数据交叉整合的结果。程序识别此数据为双端数据。
-numthreads int
线程数,默认值为 1。

常用例子:

不进行 scaffolding:
$ pasha-graph DIRECTORY -fasta in.fasta
$ pasha-graph DIRECTORY -fasta in.fasta -fasta in2.fasta
$ pasha-graph DIRECTORY -fasta in.fasta -fastq in2.fastq
$ pasha-graph DIRECTORY -fasta in.fasta -numthreads 4
进行 scaffolding:
$ pasha-graph DIRECTORY -fastaPaired in_1.fasta in_2.fasta
$ pasha-graph DIRECTORY -fastaPaired in_1.fasta in_2.fasta -fastaPaired in2_1.fasta in2_2.fasta
$ pasha-graph DIRECTORY -fastqPaired in_1.fastq in_2.fastq -numthreads 4

使用Platanus进行基因组组装

1. platanus 的安装

$ mkdir /opt/biosoft/platanus
$ wget http://platanus.bio.titech.ac.jp/Platanus_release/20130901010201/platanus -P /opt/biosoft/platanus
$ wget http://platanus.bio.titech.ac.jp/Platanus_release/20130901010201/README -P /opt/biosoft/platanus
$ chmod 755 /opt/biosoft/platanus

下载的两个文件分别是主程序和说明文件。

2. platanus 的使用

platanus 下包含三个命令,分别是 assemble, scaffold, gap_close 。其用法如下:个
这 3 个命令的共同参数为:

-t 使用的线程数,此值<=100,默认值为 1 。
-o 输出文件的前缀,默认值为 out 。

2.1 assemble

此命令基于 Bruign 图的算法来组装出 contig

-f FILE1 [File2 ...]
输入的文件,支持输入的文件总输入最大为 100 。文件可以为 fasta 或 fastq 格式。 软件会自动识别其格式。不会运用到碱基质量值,碱基质量值对组装无任何影响。
-k INT
初始的 k-mer 大小,默认值为 32 。数据覆盖度低时,该值要设小些。
-s INT
k-mer 值的步进。此值必须 >= 1,默认值为 10 。程序会使用多个 K-mer 值进行 contigs 组装。
-n INT
初始的 k-mer 覆盖度的 cutoff。 默认值为 0,即自动取值。自动取值依赖于 k-mer 的频率分布。如果其分布不正常,则应该手动设置。
-c INT
设置最小的 k-mer 覆盖度。默认值为 2 。在 k-mer 值越大的时候,则 k-mer 覆盖度越小,其 cutoff 值越小,但此 cutoff 值不能低于此参数设置的值。
-a FLOAT
K-mer 值增大的安全性水平,默认值为 10.0 。增大最终的 k-mer 值。如果牺牲准确性来延伸 contig,则设置较低的值,比如为 5.0 。
-u FLOAT
消除气泡所运行的最大差异,默认值为 0.1 。此值越大,则越容易消除气泡。特别是基因组杂合率高时,此值推荐设置更高,比如为 0.2 。
-d FLOAT
当分支的覆盖率超过此值时,则截断分支,默认值为 0.5 。此值越小,则准确率越高。如果碱基错误率较低,则适合设置较低的值,比如 0.3 。
-m INT
限制内存,单位为 GB,默认值为 16 。当程序需要消耗的内存超过此值,则会提示警告,但不会中断运行。

此程序输出的文件为

PREFIX_contig.fa  组装出的连续的序列
PREFIX_contigBubble.fa  融合并删除的气泡序列
PREFIX_kmerFrq.tsv  k-mers 频数的分布

2.2 scaffold

scaffold 用于将 paired reads 比对到 contigs 上,并确定 contigs 的顺序和方向,构建出 scaffolds 。

-c FILE1 [FILE2 ...]
contig 文件。 在此 fasta 文件的 header 中,程序识别字符 ‘cov’并将其后的数值作为覆盖度的值。即使没有 cov 信息,程序也能处理。
-b FILE1 [FILE2 ...]
Bubble_seq_file
-ip{INT} PAIR1 [PAIR2 ...]
INT 在参数的名称中,是 LIB 的名称,后接 Inward Paired 的数据文件。首尾 reads 在同一个文件中,此参数后接这个文件名。
-IP{INT} FWD1 REV1 [FWD2 REV2 ...]
INT 在参数的名称中,是 LIB 的名称, 后接 Inward Paired 的数据文件。首尾 reads 分别位于两个文件中,此参数后接这两个文件名。
-op{INT} PAIR1 [PAIR2 ...]
INT 在参数的名称中,是 LIB 的名称,后接 Outward Paired 的数据文件。首尾 reads 在同一个文件中,此参数后接这个文件名。
-OP{INT} FWD1 REV1 [FWD2 REV2 ...]
INT 在参数的名称中,是 LIB 的名称, 后接 Outward Paired 的数据文件。首尾 reads 分别位于两个文件中,此参数后接这两个文件名。
-n{INT1} INT2
INT1 在参数的名称中,是 LIB 的名称, INT2 是最小的 insert size 。 在进行 scaffolding 时,程序会自动估算各个文库的 insert size 值。如果文库中 paired reads 估算出的 insert size < IN2,则舍弃此 paired reads 信息。
-a{INT1} INT2
INT1 在参数的名称中,是 LIB 的名称, INT2 是设定的 insert size 的平均值 。设定此参数,则程序不进行自动估算 insert size 。
-d{INT1} INT2
INT1 在参数的名称中,是 LIB 的名称, INT2 是设定的 insert size 的标准差。设定此参数,则程序不进行自动估算 insert size 。
-s INT
Mapping seed length (default 32)。此值设定不能超过 reads 的长度。越小,则程序运行速度越低。
-v INT
最小的重叠长度(default 32)。 如果临近的 contigs 重叠的长度 >= INT,则这两个 contigs 则连接起来。
-l INT
连接两个 contig 所需要的最小的 paired reads 的连接数。
-u FLOAT
消除气泡所运行的最大差异,默认值为 0.1 。此值越大,则越容易消除气泡。特别是基因组杂合率高时,此值推荐设置更高,比如为 0.2 。

此程序的输出文件:

PREFIX_scaffold.fa  组装出来的 scaffold 序列
PREFIX_scaffoldBubble.fa  去除的 bubble 序列
PREFIX_scaffoldComponent.tsb  scaffold 由相应的 contigs 组成的信息

2.3 gap_close

程序将 paired reads 比对到 scaffolds,将 reads 定位到 gaps 上,并关闭一些 gap 。

-c FILE1 [FILE2 ...]
scaffold 文件
-ip{INT} PAIR1 [PAIR2 ...]
INT 在参数的名称中,是 LIB 的名称,后接 Inward Paired 的数据文件。首尾 reads 在同一个文件中,此参数后接这个文件名。
-IP{INT} FWD1 REV1 [FWD2 REV2 ...]
INT 在参数的名称中,是 LIB 的名称, 后接 Inward Paired 的数据文件。首尾 reads 分别位于两个文件中,此参数后接这两个文件名。
-op{INT} PAIR1 [PAIR2 ...]
INT 在参数的名称中,是 LIB 的名称,后接 Outward Paired 的数据文件。首尾 reads 在同一个文件中,此参数后接这个文件名。
-OP{INT} FWD1 REV1 [FWD2 REV2 ...]
INT 在参数的名称中,是 LIB 的名称, 后接 Outward Paired 的数据文件。首尾 reads 分别位于两个文件中,此参数后接这两个文件名。
-s INT
Mapping seed length (default 32)。此值设定不能超过 reads 的长度。越小,则程序运行速度越低。
-v INT
最小的重叠长度(default 32)。 此值越小,例如 20,则关闭的 gap 会越多。
-e FLOAT
重叠区所允许的最低错误率,默认值为 0.05 。 此值越大,例如 0.1 ,则关闭的 gap 会越多。

程序输出的文件:

PREFIX_gapClosed.fa  补洞后的序列文件

使用ABySS进行基因组组装

1. ABySS的安装

安装google-sparsehash
$ sudo rpm -ivh http://sparsehash.googlecode.com/files/sparsehash-2.0.2-1.noarch.rpm

$ wget http://www.bcgsc.ca/platform/bioinfo/software/abyss/releases/1.3.7/abyss-1.3.7.tar.gz
$ tar zxf abyss-1.3.7.tar.gz
$ cd abyss-1.3.7/
$ ./configure --prefix=/opt/biosoft/abyss-1.3.7/ && make -j 4 && make install
默认下ABySS支持的最大Kmer是64,在 configure 中使用 --enable-maxk=96 改变软件所能支持的最大 kmer 长度。
$ make
$ make install

$ cd ..
$ rm abyss-1.3.7/ -rf
$ echo 'PATH=$PATH:/opt/biosoft/abyss-1.3.7/bin/' >> ~/.bashrc
$ source ~/.bashrc

2. ABySS的使用

使用 ABYSS 命令能将 short reads 组装成 contigs。而如果需要组装成 scaffolds,则需要使用 abyss-pe 命令。

2.1 使用 ABYSS 将 short reads 组装成 contigs

使用如下命令查阅 ABYSS 命令的说明:

$ ABYSS --help
或者
$ less /opt/biosoft/abyss-1.3.7/share/man/man1/ABYSS.1

简要使用方法:

$ ABYSS -k 31 -o 31_contigs.fa read1.fastq reads2.fastq
注意事项: -k 和 -o 参数是必须参数; 输入文件可以有多个。

主要使用参数:

--chastity
    去除污染的reads,这是默认选项。
--no-chastity
    不去除污染的reads。
--trim-masked
    从序列某端去除低质量的碱基,这是默认选项。
--no-trim-masked
    不从序列末端去除低质量碱基。
-q | --trim-quality=N
    从序列尾端去除碱基质量低于此值的碱基。
--standard-quality
    碱基质量格式为 phred33 ,这是默认选项。
--illumina-quality
    碱基质量格式为 phred64 。
-o | --out=FILE
    输出的 contigs 文件的文件名。
-k | --kmer=N
    k-mer 长度。
-t | --trim-length=N
    maximum length of dangling edges to trim
-c | --coverage=FLOAT
    去除 k-mer 覆盖读低于此值的 contigs 。
-b | --bubbles=N
    pop bubbles shorter than N bp [3*k]
-b0 | --no-bubbles
    do not pop bubbles

2.1 使用 abyss-pe 将 short reads 组装成 scaffolds

使用如下命令查阅 ABYSS 命令的说明:

$ abyss-pe --help
或者
$ less /opt/biosoft/abyss-1.3.7/share/man/man1/abyss-pe.1

几种使用示例:

1 个 paired-end 文库:
$ abyss-pe k=64 name=ecoli in='reads1.fa reads2.fa'

多个 paired-end 文库:
$ abyss-pe k=64 name=ecoli lib='lib1 lib2' lib1='lib1_1.fa lib1_2.fa' lib2='lib2_1.fa lib2_2.fa' se='se1.fa se2.fa'

paired-end 和 mate-pair 文库:
$ abyss-pe k=64 name=ecoli lib='pe1 pe2' mp='mp1 mp2' pe1='pe1_1.fa pe1_2.fa' pe2='pe2_1.fa pe2_2.fa' mp1='mp1_1.fa mp1_2.fa' mp2='mp2_1.fa mp2_2.fa' se='se1.fa se2.fa'

使用 RNA-Seq 的组装结果进行 rescaffolding :
$ abyss-pe k=64 name=ecoli lib=pe1 mp=mp1 long=long1 pe1='pe1_1.fa pe1_2.fa' mp1='mp1_1.fa mp1_2.fa' long1=long1.fa

使用 MPI :
abyss-pe np=8 k=64 name=ecoli in='reads1.fa reads2.fa'

使用集群:
qsub -N ecoli -t 64 -pe openmpi 8 abyss-pe n=10 in='reads1.fa reads2.fa'

使用多个 k 值进行基因组组装,再寻找最佳 k 值:

$ export k
$ for k in {20..40}; do
$    mkdir k$k
$    abyss-pe -C k$k name=ecoli in=../reads.fa
$ done
$ abyss-fac k*/ecoli-contigs.fa

使用MaSuRCA进行基因组组装

1. MaSuRCA 简介

MaSuRCA(Maryland Super Read Cabog Assembler)基因组组装软件集合了 de Bruijn 和 Overlap-Layout-Consensus 的优点。
文献:Zimin A V, Marçais G, Puiu D, et al. The MaSuRCA genome assembler[J]. Bioinformatics, 2013, 29(21): 2669-2677.

2. MaSuRCA 下载和安装

$ wget wget ftp://ftp.genome.umd.edu/pub/MaSuRCA/MaSuRCA-2.2.1.tar.gz
$ tar zxf MaSuRCA-2.2.1.tar.gz -C /opt/biosoft
$ cd /opt/biosoft/MaSuRCA-2.2.1
$ ./install.sh

3. MaSuRCA 使用

3.1 配置文件准备

将模板配置文件 “/opt/biosoft/MaSuRCA-2.2.1/sr_config_example.txt” 拷贝到当前工作目录,并修改之。此配置文件含有输入文件和参数 的一些信息。内容如下:

# 测序数据的信息。分为 3 种类型:PE JUMP OTHER。每种类型的数据后接 5 列:1)2 个字符的前缀;2)平均插入片段长度;3)插入片段长度标准差;4)fastq(.gz)格式的 reads1; 5)fastq(.gz)格式的 reads2。如果有 jump 数据是 FR 类型,则,则使用 JUMP,但是平均插入片段长度为负数。其它的数据,则必须要转换成 Celera 兼容的 .frg 文件。
DATA
PE= p1 180 20 180_1.fastq 180_2.fastq
PE= p2 500 50 500_1.fastq 500_2.fastq
JUMP= j1 2000 200 2000_1.fastq 2000_2.fastq
JUMP= j2 5000 500 5000_1.fastq 5000_2.fastq
OTHER= file.frg
END

PARAMETERS
# 设置 k-mer size,大小为 25~101,或者为 auto,表示自动计算最优值。
GRAPH_KMER_SIZE=auto
# 如果仅分析 Illumina 数据,则值为 1;如果有 1x 及以上的 454 数据,则设置为 0。
USE_LINKING_MATES=1
# 如果 jumping library 的数据过多,可能会 confuse the assembler,设置此值为 60,则仅使用 60x 左右的 jumping 数据用于基因组组                  装。对于细菌基因组,一般设置为 60。如果基因组较大,则设置此值大些。对于一些较大的真核基因组,可以大至 1000。
LIMIT_JUMP_COVERAGE = 60
# Celera Assembler 的参数。如果是 mammals 的基因组,cgwErrorRate的值不能高于 0.15。
CA_PARAMETERS = ovlMerSize=30 cgwErrorRate=0.25 ovlMemory=4GB
# 舍弃频数低于此值的 k-mer。如果覆盖度大于 100,可以设置此值为 2。
KMER_COUNT_THRESHOLD = 1
# 设置使用的线程数。
NUM_THREADS= $NUM_THREADS
# 设置 jellyfish 的 hash size。此值可以设置为 "基因组大小+reads的数目"。
JF_SIZE=100000000
# 设置是否 trim long reads 的 3' homopolymers(e.g. GGGGGGG)。适合于高 GC 含量的基因组。
DO_HOMOPOLYMER_TRIM=0
END

3.2 运行 masurca 和 assemble.sh 进行基因组组装

运行程序 masurca,生成 assemble.sh; 然后运行 assemble.sh 进行组装。

$ /opt/biosoft/MaSuRCA-2.2.1/bin/masurca config.txt
$ ./assemble.sh

3.3 运行中断后继续运行

由于程序出错,或手动终止后,可以终止步骤所生成的文件,在继续运行 masurca ,生成含有后续步骤的 assemble.sh,再继续运行程序。

4. 结果文件

最终的结果文件为 CA/10-gapclose/genome.ctg.fasta 。

bwa的使用

1. 简介

最近需要使用软件GAM-NGS软件,但是该软件貌似不支持bowtie2的结果。可能需要用到BWA,于是参考BWA的网站其Manual
bwa,即Burrows-Wheeler-Alignment Tool。
BWA 包含 3 种算法:

BWA-bactrack: 用于进行 Illumina reads 的比对。reads 的长度最大为 100bp。
BWA-SW: 用于比对 long-read ,支持的长度为 70bp-1Mbp;同时支持剪接性比对。
BWA-MEM: 最新的算法。和 BWA-SW 的适用性一致,但是更加快速和准确;同时与 BWA-bactrack 相比,在对 70-100bp reads 的比对上有更优的性能。

2. BWA 的下载和安装

$ wget http://jaist.dl.sourceforge.net/project/bio-bwa/bwa-0.7.9a.tar.bz2
$ tar jxf bwa-0.7.9a.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/bwa-0.7.9a/
$ make
$ echo 'PATH=$PATH:/opt/biosoft/bwa-0.7.9a' >> ~/.bashrc
$ source ~/.bashrc

3. BWA用法

2.1 index

在进行 reads 的比对前,需要对 fasta 文件构建 FM-index。
常用例子和参数如下:

$ bwa index ref.fa -p genome

-p STR
输出的数据库前缀。 默认和输入的文件名一致,则输出的数据库在其输入文件所在的文件夹,并以该文件名为前缀。
-a [is|bwtsw]
输入构建 index 的算法。 is 算法简单快速,是默认的选项,但是不能用于基因组大于 2GB 的数据库; bwtsw 适合于大基因组,比如人类基因组。

2.2 mem

该算法先使用 MEM(maximal exact matches) 进行 seeding alignments,再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。
BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配,导致 reads 有多个最佳匹配位点。这对 long reads 的比对时比较重要的结果。但是却会和 Picard 的 markDuplicates 程序部兼容。
常用例子和参数如下:

$ bwa mem genome reads.fq > aln-se.sam
$ bwa mem genome read1.fq read2.fq > aln-pe.sam

-t INT
使用的线程数
-p
若无此参数:输入文件只有1个,则进行单端比对;若输入文件有2个,则作为paired reads进行比对。若加入此参数:则仅以第1个文件作为输入(输入的文件若有2个,则忽略之),该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-M
加入此参数用于将 shorter split hits 标记为次优,有利于兼容 Picard。

2.3 bwasw

对输入的第1个文件的所有序列进行比对。如果输如有 2 个文件,则进行 paired-end 比对,此模式仅对 Illumina 的 short-insert 数据进行比对。在 Paired-end 模式下,BWA-SW依然输出剪接性比对结果,但是这些结果会标记为 not properly paired; 同时如果有多个匹配位点,则不会写入 mate 的匹配位置。
常用例子和参数如下:

$ bwa bwasw genome long_read.fq > aln.sam
$ bwa bwasw genome read1.fq read2.fq > aln-pe.sam

-t INT
使用的线程数

2.4 backtrack

经典的 bwa 先使用 aln 命令将单独的 reads 比对到参考序列,再使用 samse 或 sampe 生成 sam 文件。
常用例子:

$ bwa aln genome read1.fq > aln_sa1.sai
$ bwa aln genome read2.fq > aln_sa2.sai
$ bwa samse genome aln_sa1.sai read1.fq > aln_se.sam
$ bwa sampe genome aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam