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 安装

安装 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