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