使用 SGA 进行基因组组装

1. SGA 简介

SGA (String Graph Assembler), 采用 String Grapher 的方法来进行基因组的组装。软件通过创建 FM-index/Burrows-Wheeler 索引来进行查找 short reads 之间的 overlaps, 从而进行基因组组装。
使用 SGA 的注意事项:

1. SGA 的输入文件为 fastq 文件,序列长度推荐为 100bp 及以上。较短的序列长度,则使用 De Bruijn 的方法进行组装能或得更好的基因组序列。
2. 使用 SGA 进行基因组组装需要最少 40X

2. SGA 安装

SGA 的安装与使用需要依赖 google-sparsehash, bamtools, zlib 和 jemalloc 等软件。

安装 goolge-sparsehash,由于国内屏蔽了google,无法下载该软件。
$ wget www.chenlianfu.com/public/software/sparsehash-2.0.2.tar.gz
$ tar zxf sparsehash-2.0.2.tar.gz
$ cd sparsehash-2.0.2
$ ./configure && make -j 8 && sudo make install

安装 bamtools
$ wget https://github.com/pezmaster31/bamtools/archive/master.zip -O bamtools.zip
$ unzip bamtools.zip
$ cd bamtools-master/
$ mkdir build; cd build
$ cmake ../
$ make -j 8
$ cd ../../
$ mv bamtools-master/ /opt/biosoft/bamtools

安装 jemalloc,可选,有利于减少运算时间。
$ wget https://github.com/jemalloc/jemalloc/archive/dev.zip -O jemalloc.zip
$ unzip jemalloc.zip
$ ./autogen.sh 
$ ./configure && make -j 8 && sudo make install
$ cd ../
$ mkdir /opt/biosoft/sga/
$ mv jemalloc-dev /opt/biosoft/sga/

安装 SGA
$ wget https://github.com/jts/sga/archive/master.zip -O sga.zip
$ unzip sga.zip
$ cd sga-master/src
$ ./autogen.sh
$ ./configure --with-bamtools=/opt/biosoft/bamtools --with-jemalloc=/opt/biosoft/sga/jemalloc-dev/lib/ --prefix=/opt/biosoft/sga/
$ make -j 8
$ make install
$ cd ..
$ mv * /opt/biosoft/sga/
$ echo 'PATH=/opt/biosoft/sga/bin/' >> ~/.bashrc
$ source ~/.bashrc

此外,如果需要使用程序附带的一些 python 程序,则需要安装一些 python 模块。
$ wget https://pypi.python.org/packages/source/r/ruffus/ruffus-2.6.3.tar.gz
$ tar zxf ruffus-2.6.3.tar.gz
$ cd ruffus-2.6.3
$ sudo python setup.py install

3. SGA 的使用

SGA 包含很多个命令,直接输入命令,即可得到帮助信息。
以 2 个 PE 文库和 2 个 MP 文库的测序数据为输入,使用 SGA 进行基因组组装,来介绍运行流程。
输入文件:

pe180.1.fastq    pe180.2.fastq
pe500.1.fastq    pe500.2.fastq
mp2000.1.fastq   mp2000.2.fastq
mp5000.1.fastq   mp5000.2.fastq

3.1 数据预处理

使用 preprocess 程序进行低质量数据的截短和删除,示例:

去除 PE 文库数据中含有 N 的reads。此处仅对需要进行 contig 组装的 PE 数据进行处理。
$ sga preprocess -o pe180.pp.fastq --pe-mode 1 pe180.1.fastq pe180.2.fastq
$ sga preprocess -o pe500.pp.fastq --pe-mode 1 pe500.1.fastq pe500.2.fastq

结果中 pp 表示 preprocess 的意思。
一般用其它软件进行 Fastq 文件的 QC,此处使用 preprocess 去除含有 N 的 reads 即可。由于软件后续有对 reads 的校正过程,因此不推荐对 reads 进行截短处理。

preprocess 命令参数:

sga preprocess [OPTION] READS1 READS2 ...

-o | --output FILE
    将输出数据写入到该文件中。默认输出到标准输出。程序仅得到一个输出文件。成对的数据用 /1 /2 标识出来并在输出文件中交叉排列。
-p | --pe-mode INT
    paired end 模式。 有 3 中模式,用数字 0, 1, 2 表示。
    0:表示认为输入数据都是单端数据;
    1:表示 READS1 中的数据和 READS2 中的数据是成对的,这种情况下若有一条 read 被去除,则对应的另外一条也会被去除;
    2:表示 reads 是成对的,可能该值适合于3个及以上文件作为输入(我没有用过)。
--pe-orphans FILE
    如果 read pair 中有一条序列被去除,则另外一条序列写入到此文件中。
--phred64
    加入该参数,表明输入的 fastq 文件为 phred64 格式,输出结果会转换成 phred33 格式。
--discard-quality
    不输出碱基质量。
-q | --quality-trim INT
    对低质量碱基进行 trim。
-f | --quality-filter INT
    对低质量 reads 进行去除。若 read 中碱基质量 <=3 的低质量碱基多余该值,则去除整条 read。
-m | --min-length INT
    去除长度低于此值的 read。
-h | --hard-clip INT
    将所有的 reads 长度截短到该值。
--permute-ambiguous
    将 reads 中的 N 随机替换为 ATCG 中的一种。
--dust
    去除低复杂度 reads 。
-dust-threshold FLOAT
    去除 dust score 高于此值的 reads。默认为 4.0 。
--suffix STRING
    在每个 read ID 后添加后缀。
--no-primer-check
    不进行 primer 序列检测
-r | --remove-adapter-fwd  STRING
    去除 adapter 序列,提供 adapter 序列的第一条链。
-c | --remove-adapter-rev STRING
    去除 adapter 序列,提供 adapter 序列的第二条链。

2. FM 索引构建与合并

使用 index 命令对 reads 数据构建索引。当使用多组数据进行 contig 组装的时候,使用 merge 对多组数据的索引进行合并。示例:

$ sga index -a ropebwt -t 8 --no-reverse pe180.pp.fastq
$ sga index -a ropebwt -t 8 --no-reverse pe500.pp.fastq

$ sga merge -t 8 -p pe -r pe180.pp.fastq pe500.pp.fastq

index 命令常用参数:

sga index [OPTION] ... READSFILE

-a | --algorithm STRING
    BWT 构建算法,有两种 sais (默认) 和 ropebwt。前者适合与较长的 reads;后者适合与较短的(<200bp) reads,但是速度快,消耗内存少。
-t | --threads INT
    运算的线程数。默认为 1 。
-p | --prefix
    指定输出文件前缀,默认为输入文件的前缀。
--no-reverse
    不生成反向序列的 BWT 文件。由于对测序数据进行校正的时候不许要该文件,一般加入该参数以节约资源。
--no-forward
    不生成正向序列的 BWT 文件。
--no-sai
    不生成 SAI 文件。

merge 命令常用参数:

sga merge [OPTION] ... READS1 READS2

-t | --threads INT
    运算的线程数。默认为 1 。
-p | --prefix 
    指定输出文件前缀,默认为输入文件的前缀。
-r | --remove
    生成合并后的 index 文件后,删除旧的 index 文件和 fastq 文件。

3. reads 的校正

使用 correct 命令对 reads 进行校正,示例:

$ sga correct -t 8 -k 41 -x 2 --learn -o pe.ec.fastq pe.fastq

correct 命令参数:

sga correct [OPTION] ... READSFILE

-p | --prefix STRING
    输入的 index 文件的前缀,默认为输入文件的前缀。
-o | --outfile STRING
    将校正后的数据写入到此文件中,默认输出的文件名为 READSFILE.ec.fa。虽然后缀是 fa,但其实是 fastq 文件。
-t | --threads INT
    运算的线程数。默认为 1 。
--discard
    检测并去除低质量 reads 。
-a | --algorithm STRING
    进行 reads correction 的算法。有 3 种:
    kmer: 使用 kmer 的方法进行校正,默认是该值;
    overlap:基于重叠的方法进行校正;
    hybrid:应该是同时上面 2 中方法进行校正。
--metrics FILE
    校正后,将 read 中每个位点的碱基错误率写入到该文件。

使用 kmer 方法进行校正的参数:
-k | --kerm-size INT
    选取的 kemr 长度,默认为 31 。
-x | --kmer-threshold INT
    对出现次数低于此值的 kmer 进行校正。默认为 3 。
-i | --kmer-rounds INT
    执行 N 轮校正,最多校正 N 个碱基。默认的 N 值为 10。
--learn
    程序自行计算 -x 的参数。若测序覆盖度低于 20x 或测序数据在基因组上的覆盖不均匀,则不适合选择该参数。

使用 overlap 方法进行校正的参数:
-e | --error-rate    default:0.04
    2 条 reads 有 overlap 时,允许的最大错误率。
-m | --min-overlap    default:45
    2 条 reads 重叠的最小碱基数。
-c | --conflict    default:5
    检测出一个错误碱基所需要的最小重叠数目。
-r | --rounds    default:1
    进行迭代校正的指定的次数。

4. reads 的过滤

去除数据中完全一模一样的 duplicate;去除含有低频 kmer 的 reads。示例:

需要先对校正过后的数据重新构建 index
$ sga index -a ropebwt -t 8 pe.ec.fastq

$ sga filter -x 2 -t 8 -o pe.ec.f.fastq --homopolymer-check --low-complexity-check pe.ec.fastq

filter 命令常用参数:

sga filter [OPTION] ... READSFILE

-o | --outfile STRING
    将校正后的数据写入到此文件中,默认输出的文件名为 READSFILE.filter.pass.fa 。
--no-duplicate-check
    不进行 duplicate removal
--no-kmer-check
    不进行 kmer check
--kmer-both-strand
    对 reads 的正负链都进行检测
--homopolymer-check
    检测单碱基连续情况
--low-complexity-check
    去除低重复 reads
-k | --kmer-size    default:27
    kmer 检测时所用的 kmer 值。
-x | --kmer-threshold     default:3
    保留的 read 中所有的 kemr 的覆盖度必须大于该值。

5. 整合出简单的序列

使用 fm-merge 整合出一些简单序列(simple unbranched chains),示例:

$ sga fm-merge -m 55 -t 8 -o merged.fa pe.ec.f.fastq

fm-merge 命令参数:

sga fm-merge [OPTION] ... READSFILE

-m | --min-overlap    default:45
    2 个 reads 之间最小的 ovlerlap 长度。此值设置过小,则消耗更多的计算时间。

6. 去除重复序列

在上一步命令之后,生成的序列中,有很多的包含与被包含的关系,使用 rmdup 命令去除 substrings 序列。示例:

$ sga index -d 1000000 merged.fa
$ sga rmdup -t 8 -o rmdup.fa merged.fa

rmdup 命令参数:

sga rmdup [OPTION] ... READFILE

-o | --out FILE    default:READFILE.rmdup.fa
    设置输出文件名。
-e | --error-rate    default: exact matches required
    设定 2 条序列一致,这 2 条之间的允许的最大错误率。默认 2 条序列必须完全一致。

7. 构建 string grahp

使用 overlap 命令来构建 string graph。如果组装小物种的基因组,可以跳过 fm-merge 和 rmdup 这 2 步。示例:

$ sga overlap -m 55 -t 8 rmdup.fa

overlap 命令常用参数:

-m | --min-overlap    default:45
    2 个 reads 之间最小的 ovlerlap 长度。此值设置过小,则消耗更多的计算时间。

7. 组装 contigs

使用命令 assemble 进行 contigs 组装,示例:

$ sga assemble -m 75 -g 0 -r 10 -o assemble.m75 rmdup.asqg.gz

assemble 命令常用参数:

-m | --min-overlap    default:45
    2 个 reads 之间最小的 ovlerlap 长度。当参数在 fm-merge, overlap 和 assemble 这 3 个命令中都有出现。当测序 reads 长度为 100bp 时, 该参数的值在 3 个命令中比较合适的值分别是 65, 65, 75。当然, reads 长度越长,则该参数的值设置越大。此外,一般通过调节 overlap 中该参数的值来优化组装结果。
-d | --max-divergence    default:0.05
    如果序列中的 SNP 比率 < 0.05,则消除变异。当基因组杂合率非常高时,需要增大该值。
-g | --max-gap-divergence    default:0.01
    如果序列中的 INDEL 比率 < 0.01,则消除变异。当基因组杂合率非常高时,需要增大该值。如果设置改值为 0, 则表示不消除变异。
--max-indel    default:20
    当 indel 长度大于该值的时候,则不消除变异。
-l | --min-branch-length    default:150
    去除 tip 序列,如果其序列长度低于该值。当 reads 长度为 100 bp 时候,推荐设置为 150;若 reads 长度越大,该值也需要越大。
-r | --resolve-small    default: not perfomed
    对 samll repeat 序列进行打断,如果该跨过该区域的 overlaps 之间的最大长度差超过该值时。可以设置该值为 10 。

7. 组装 scaffold

首先,需要将 PE 和 MP 数据比对到 contigs 序列上

$ /opt/biosoft/sga/src/bin/sga-align --name pe180 assemble.m75-contigs.fa pe180.1.fastq pe180.2.fastq
$ /opt/biosoft/sga/src/bin/sga-align --name pe250 assemble.m75-contigs.fa pe250.1.fastq pe250.2.fastq
$ /opt/biosoft/sga/src/bin/sga-align --name mp2000 assemble.m75-contigs.fa mp2000.1.fastq mp2000.2.fastq
$ /opt/biosoft/sga/src/bin/sga-align --name mp5000 assemble.m75-contigs.fa mp5000.1.fastq mp5000.2.fastq

然后,计算 contig 之间的距离

$ sga-bam2de.pl -n 5 --prefix pe180 pe180.bam
$ sga-bam2de.pl -n 5 --prefix pe500 pe500.bam
$ sga-bam2de.pl -n 5 --prefix mp2000 mp2000.bam
$ sga-bam2de.pl -n 5 --prefix mp5000 mp5000.bam

再计算 contig 的拷贝数

$ sga-astat.py -m 200 pe180.refsort.bam > pe180.astat
$ sga-astat.py -m 200 pe500.refsort.bam > pe500.astat

最后,进行 scaffold 组装

$ sga scaffold -m 200 -a pe180.astat -a pe500.astat -pe pe180.de --pe pe500.de --mp mp2000.de --mp mp5000.de -o genome.scaf assemble.m75-contigs.fa
$ sga scaffold2fasta -m 200 -a assemble.m75-graph.asqg.gz -o genomeScaffold.n5.fa -d 1 --use-overlap --write-unplaced genome.scaf

发表评论

电子邮件地址不会被公开。 必填项已用*标注