1. DBG2OLC软件简介

DBG2OLC能利用二代和三代混合数据组装大基因组。其文章于2016年发表在Scientific Reports上。

2. DBG2OLC软件下载与安装

使用git下载DBG2OLC软件

$ cd /opt/biosoft/
$ git clone https://github.com/yechengxi/DBG2OLC.git
$ cd /opt/biosoft/DBG2OLC
按照说明中对软件进行编译，编译出的3个可执行程序全部都是DBG2OLC命令
$ g++ -O3 -o SparseAssebmler DBG2OLC.cpp
$ g++ -O3 -o DBG2OLC *.cpp
$ g++ -O3 -o Sparc *.cpp
直接拷贝作者编译好的程序即可
$ chmod 755 compiled/*
$ cp compiled/* .
$ echo 'PATH=$PATH:/opt/biosoft/DBG2OLC' >> ~/.bashrc
$ source ~/.bashrc

DBG2OLC程序第三步需要blasr, sparc/pbdagcon软件支持。其中sparc在DBG2OLC安装文件夹下。
安装blasr

下载BLASR
$ git clone https://github.com/PacificBiosciences/blasr.git /opt/biosoft/blasr
$ cd /opt/biosoft/blasr/
下载libcpp和pbbam两个submodules
$ make update-submodule

blasr编译需要hdf5支持,从hdf5官网下载适合centos6的二进制包并安装 
$ wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0-patch1/bin/linux-centos6-x86_64-gcc447/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz
$ tar zxf hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz -C /opt/sysoft/

可以使用cmake, pitchfork和make三种方式对BLASR进行编译，以下使用常规的make方法进行编译，需要高版本gcc支持
对BLASR进行configure
$ ./configure.py --shared --sub --no-pbbam HDF5_INCLUDE=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/include/ HDF5_LIB=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/
对submodules进行configure
$ make configure-submodule
对submodules进行make
$ make build-submodule -j 4
对BLASR进行make
$ make blasr -j 4
对其它工具，例如pls2fasta, loadPulses, sawriter等进行编译，其结果文件在utils文件夹中
$ make -j 4

可选步骤：手动将有用的命令和库文件放置到指定的地方
blasr的正常运行需要依赖libcpp里面三个库文件和hdf5软件中的库文件
$ mkdir bin lib
$ cp blasr bin/
$ find utils -maxdepth 1 -perm 775 -exec cp {} bin/ \;
$ cp ./libcpp/pbdata/libpbdata.so ./libcpp/hdf/libpbihdf.so ./libcpp/alignment/libblasr.so lib/
$ echo 'export LD_LIBRARY_PATH=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/:/opt/biosoft/blasr/lib/:$LD_LIBRARY_PATH
PATH=/opt/biosoft/blasr/bin/:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

若DBG2OLC流程第三步选择使用pbdagcon进行运算，则需要安装pbdagcon软件

pbdagcon软件的编译需要高版本gcc支持
$ git clone https://github.com/PacificBiosciences/pbdagcon.git /opt/biosoft/pbdagcon
$ cd /opt/biosoft/pbdagcon
$ ./configure.py --boost --gtest --sub --no-pbbam
$ make init-submodule
$ make -j 4
$ make check
$ mkdir bin
$ cp src/cpp/dazcon src/cpp/pbdagcon bin/
$ echo 'PATH=/opt/biosoft/pbdagcon/bin:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

3. 程序运行

使用DBG2OLC软件利用二代和三代数据混合的基因组组装，其运行流程分3步。

3.1 使用SparseAssembler利用二代数据进行DBG组装

首先，利用Illumina小片段文库数据使用SparseAssembler命令组装出contigs序列。此外，也可以使用其他基因组组装软件组装出contigs序列后，直接跳到DBG2OLC的第二个步骤。值得注意的是：输入到第二步骤的contigs必须是没有经过repeat resolving的原始序列；绝大部分基因组组装软件为了获得更完整连续的contigs序列，都牺牲了contigs的准确性，其结果不能用于DBG2OLC软件的第二步，否则最终结果会很差；作者推荐可以直接用于第二步的其它contig组装软件有Platanus和Meraculous。
一般情况下，输入到SparseAssembler命令中~50x的Illumina数据，能获得较好的contigs结果。
SparseAssembler命令参数：

常用参数：
LD 
    是否载入k-mer graph。第一次运行SparseAssembler命令的时候，该参数的值必须是0；若为了使用SparseAssembler得到更好的contigs结果，则需要调整NodeCovTh和EdgeCovTh参数；调整这些参数的时候，不需要再次计算k-mer graph，设置该参数为1来跳过这个步骤，从而节约很多时间。
k 
  设置使用DBG方法计算时的Kmer大小，支持的Kmer大小为15-127。
g 
  number of skipped intermediate k-mers, support 1-25.该参数软件示例中使用的值是15。
NodeCovTh 
  设置k-mers覆盖度阈值，去除覆盖度较低的k-mers。该值设定范围为0-16，默认值为1。
EdgeCovTh 
  设置link覆盖度阈值，去除覆盖度较低的links。该值设定范围为0-16，默认值为0。
GS 
  设置一个基因组大小的值。该参数用于决定预先占用的内存量。推荐设置得比基因组大，例如设置为2倍基因组大小。
f 
  输入单端测序数据的路径。输入文件可以是fasta或fastq文件。若有多个输入文件，则使用多个f参数。
i1  i2 
  输入inward paired-end数据。若有多组paired-end数据，则多次使用i1/i2参数。

其它参数：
o1  o2 
  输入outward paired-end数据。
i1_mp  i2_mp 
  输入插入片段长度>10kb的inward paired-end数据。
o1_mp  o2_mp 
  输入插入片段长度>10kb的outward paired-end数据。
PathCovTh 
  设置path覆盖度阈值，去除覆盖度较低的paths。该值设定范围为0-100。根据经验，不推荐添加该参数。
TrimLen 
  将所有过长的序列截短到此指定的长度。
TrimN 
  若read中的碱基N数目超过此设定的值，则去除该read数据。
TrimQual 
  从尾部截短质量低于此值的碱基。
QualBase 
  设置Fastq文件中最低碱基质量对应的ASCII码符号。默认值是'!'，等同于Pred33。
Denoise 
  设置是否对reads进行修正。默认值是0，表示不对reads进行修正。
H 
  混合模式。默认值是0，表示对reads的尾部进行截短，来保证高质量的数据进行reads修正。
CovTh 
  覆盖度 < 此设定值的k-mer会被检测，从而被校正。若该参数值设置为0，则软件会自动计算该值。
CorrTh 
  覆盖度 >= 此设定值的k-mer可以用来对reads做校正。若该参数值设置为0，则软件会自动计算该值。

SparseAssembler运行示例：

对某物种Illumina小片段文库测序的PE150bp数据使用trimmomatic质控，再使用FindErrors进行修正，再运行SparseAssembler：
$ SparseAssembler LD 0 k 95 g 15 NodeCovTh 1 EdgeCovTh 0 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq
$ cp Contigs.txt Contigs.txt.00
增大NodeCovTh和EdgeCovTh参数后，再次运行SparseAssembler，并比较两次结果。第二次运行较第一次运行，耗时少了很多很多。
$ SparseAssembler LD 1 k 95 g 15 NodeCovTh 2 EdgeCovTh 1 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq

SparseAssembler在当前目录下生成了18个文件结果，其中Contigs.txt文件是Fasta格式的Contigs序列。
运行SparseAssembler的注意事项：

1. SparseAssembler只可以简单地对Fastq文件进行质控和错误修正。推荐使用其它软件进行reads质控和修正，以获得更好的结果。
2. 参数k设置了k-mer的大小，该参数的值对结果的影响较大。若基因组较小，推荐设置多个k-mer值进行多次计算，从而选择最优k-mer值。个人经验，PE150bp数据的最优的k-mer值约为91~99。
3. 选定了k-mer大小后，使用默认的NodeCovTh和EdgeCovTh参数（默认参数一般能得到很好的结果）运行一遍SparseAssembler。然后尝试增大NodeCovTh和EdgeCovTh参数，设置LD 1参数再次运行SparseAssembler，以获得最优的Contigs结果。
4. 可能是先使用了最小的NodeCovTh和EdgeCovTh参数做运算后，才能再次使用更大的参数进行运算。
5. SparseAssembler虽然也有输入大片段文库数据的参数和Scaffolding参数，但是不推荐输入大片段文库数据进行Scaffolding操作，没太大意义。
6. 虽然SparseAssembler命令的文件输入方式有多种，若是仅进行contigs组装，没有利用到paired信息，因此使用i1 i2参数输入文件和使用f参数输入文件的结果是一模一样的。

3.2 使用DBG2OLC找Contigs序列和Pacbio reads的Overlap并进行Layout

DBG2OLC通过比较contigs和Pacbio reads之间的overlap，将contigs序列定位到Pacbio reads上，将DBG的contigs结果运用到OLC算法中。
DBG2OLC命令参数：

主要参数：
LD 
  是否载入compressed reads information。第一次运行DBG2OLC命令的时候，该参数的值必须是0；若为了得到更好的结果，则需要调整其它参数；调整这些参数的时候，设置该参数为1来跳过这个步骤，从而节约很多时间。
k 
  设置k-mer大小。k-mer用来比较contig和pacbio read之间的重叠，而不是用于基因组组装，推荐设置为 17 即可。
AdaptiveTh 
  若contig和pacbio read之间匹配的k-mers数目 < AdaptiveTh * contig长度，则认为contig和pacbio read没有重叠。推荐设置为0.001-0.02。
KmerCovTh 
  若contig和pacbio read之间匹配k-mers的覆盖度 < KmerCovTh，则认为contig和pacbio read没有重叠。推荐设置为2-10。
MinOverlap 
  两条Pacbio read之间匹配的k-mers数目 < MinOverlap，则认为它们之间没有重叠。推荐设置为10-150。
RemoveChimera 
  去除嵌合体Pacbio reads。若Pacbio数据覆盖度大于10x，推荐设置该参数为 1 。
Contigs 
  输入contigs序列文件路径
f 
  输入Pacbio测序Fasta/Fastq文件路径。

其它参数：
MinLen 
  设置能用于计算的最小Pacbio reads长度。
ChimeraTh 
  该参数默认值是 1 ；若Pacbio数据覆盖度大于100x，则推荐加入该参数并设置为 2 。
ContigTh 
  该参数默认值是 1 ；若Pacbio数据覆盖度大于100x，则推荐加入该参数并设置为 2 。

DBG2OLC运行示例：

$ DBG2OLC LD 0 k 17 AdaptiveTh 0.001 KmerCovTh 2 MinOverlap 20 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq
$ DBG2OLC LD 1 k 17 AdaptiveTh 0.005 KmerCovTh 3 MinOverlap 30 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq

DBG2OLC的结果文件很多，其主要结果文件是backbone_raw.fasta和DBG2OLC_Consensus_info.txt，是第三步的输入文件。
运行DBG2OLC的注意事项：

1. AdaptiveTh, KmerCovTh和minOverlap这3个计算Overlap的参数对结果的影响最大。对于10x/20x PacBio数据：KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01；对于50x-100x PacBio数据：KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02。
2. 不推荐对Pacbio数据就行修正后再运行DBG2OLC。可以比较使用修正前的数据用于DBG2OLC的结果，一般情况下使用未修正的Pacbio数据能获得更好的结果。此外，DBG2OLC运行过程中有一步多序列比对模块来进行错误修正。
3. 可能是先使用了最小的AdaptiveTh, KmerCovTh和minOverlap参数做运算后，才能再次使用更大的参数进行运算。

3.3 Call consensus

本步骤需要使用/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本，来调用blasr和consensus模块Sparc（可以考虑使用pbdagcon作为consensus模块，但DBG2OLC没有提供相应的脚本）进行运算。

先将/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本拷贝到当前目录
$ cp /opt/biosoft/DBG2OLC/utility/*.sh /opt/biosoft/DBG2OLC/utility/*.py ./
若使用了最新版本的blasr软件，其参数书写方法有一个中划线变成了两个中划线，因此需要修改.sh文件中blasr命令的参数书写方法。
此外，也需要修改.sh文件中Sparc命令路径，或者将Sparc命令也拷贝到当前目录。

将Contigs序列和Pacbio reads数据合并成一个文件ctg_pb.fasta
$ cp Contigs.txt ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell01.fastq >> ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell02.fastq >> ctg_pb.fasta

运行脚本程序split_and_run_sparc.sh
$ ./split_and_run_sparc.sh backbone_raw.fasta DBG2OLC_Consensus_info.txt ctg_pb.fasta ./ 2 > cns_log.txt
结果会输出到 ./ 目录下。最后的结果文件是final_assembly.fasta。

AGP文件为NCBI数据上传要求的标准格式，用来描述小片段序列（比如contig）如何构成大片段序列（比如scaffold和chromosome）。详细的说明文档请见：http://www.ncbi.nlm.nih.gov/projects/genome/assembly/agp/AGP_Specification.shtml
AGP文件有9列，分别是：

1. 大片段的序列名(object)
2. 大片段起始(object_begin)
3. 大片段结束(object_end)
4. 该段序列在大片段上的编号(part_number)
    一般一个大片段由多个小片段和gap组成。此处则为这些小片段和gap在大片段上的编号。
5. 该段序列的类型(component_type)
    常用的是W、N和U。W表示WGS contig；N表示指定大小的gap；U表示不明确长度的gap，一般用100bp长度。
6. 小片段的ID或gap长度(component_id or gap_length)
    如果第5列不为N或U，则此列为小片段的ID。
    如果第5列是N或U，则此列为gap的长度。如果第5列为U，则此列值必须为100。
7. 小片段起始或gap类型(component_begin or gap_type)
    如果第5列是N或U，则此列表示gap的类型。常用的值是scaffold，表示是scaffold内2个contigs之间的gap。其它值有：contig，2个contig序列之间的unspanned gap，这样的gap由于没有证据表明有gap，应该要打断大片段序列；centromere，表示中心粒的gap；short_arm，a gap inserted at the start of an acrocentric chromosome；heterochromatin，a gap inserted for an especially large region of heterochromatic sequence；telomere，a gap inserted for the telomere；repeat，an unresolvable repeat。
8. 小片段结束或gap是否被连接(component_end or linkage)
    如果第5列是N或U，则此列一般的值为yes，表示有证据表明临近的2个小片段是相连的。
9. 小片段方向或gap的连接方法(orientation or linkage_evidence)
    如果第5列不为N或U，则此列为小片段的方向。其常见的值为 +、-或？。
    如果第5列是N或U，则此列表明临近的2个小片段能连接的证据类型。其用的值是paired-ends，表明成对的reads将小片段连接起来。其它值有：na，第8列值为no的时候使用；align_genus，比对到同属的参考基因组而连接；align_xgenus，比对到其它属的参考基因组而连接；align_trnscpt，比对到同样物种的转录子序列上；within_clone，gap两边的序列来自与同一个clone，但是gap没有paired-ends跨越，因此这种连接两边小片段无法确定方向和顺序；clone_contig，linkage is provided by a clone contig in the tiling path (TPF)；map，根据连锁图，光学图等方法确定的连接；strobe，根据PacBio序列得到的连接；unspecified。如果有多中证据，则可以写上多种证据，之间用分号分割。

例子：
Scaffold from component (WGS)
Chromosome from scaffold (WGS)

1. ABySS 进行 scaffolding 的目的与优点

目的： 对其它基因组 denovo 的 assembly 结果，使用 abyss 再进行一次 scaffolding。
优点： 可以使用 RNA-seq 的转录子数据进行基因组的辅助组装。

2. ABySS 进行 scaffolding 的命令行

输入文件： assembly.fasta, 2000.1.fastq, 2000.2.fastq, 5000.1.fastq, 5000.2.fastq。
输入文件是基因组的组装结果，和 3 对 mate-paired Illumina 数据。

2.1 对 assembly.fasta 进行序列改名

去除序列之间的换行
fasta_no_blank.pl assembly.fasta > 11; mv 11 assembly.fasta
给序列按顺序重命名
perl -e '$num = 0; while (<>) {if (/^>/) { s/>(.*)/>$num/; print; $num ++; } else { print } }' assembly.fasta > ledodes-6.fa

2.2 将 mate-paired 数据比对到基因组序列上

根据比对结果，得到 mate-paired library 的 insertSize 信息（以 .hist 为后缀的文件）和 序列之间的连接、距离与顺序信息 (以 .dist.dot 为后缀的 graph 文件）。

abyss-map -j24 -l87 3000.1.fastq 3000.2.fastq ledodes-6.fa \
  |abyss-fixmate -l87 -h mp1-6.hist \
  |sort -snk3 -k4 \
  |DistanceEst --dot -j24 -k87 -l87 -s200 -n10 -o mp1-6.dist.dot mp1-6.hist
abyss-map -j24 -l87 8000.1.fastq 8000.2.fastq ledodes-6.fa \
  |abyss-fixmate -l87 -h mp2-6.hist \
  |sort -snk3 -k4 \
  |DistanceEst --dot -j24 -k87 -l87 -s200 -n10 -o mp2-6.dist.dot mp2-6.hist

以上命令行中参数：

-j24
    使用 24 个线程运行
-l87
    使用的 kmer值 为 87
-s200
    sedd contigs的最小长度为 200bp
-n10
    所允许连接两条序列的最小的pairs的数目

2.3 进行 scaffolding

abyss-scaffold -k87 -s200 -n5 -g ledodes-6.path.dot ledodes-6.fa mp1-6.dist.dot mp2-6.dist.dot > ledodes-6.path
PathConsensus -k87 -p0.9 -s ledodes-7.fa -g ledodes-7.adj -o ledodes-7.path ledodes-6.fa ledodes-6.fa ledodes-6.path
cat ledodes-6.fa ledodes-7.fa \
  | MergeContigs -k87 -o ledodes-8.fa - ledodes-7.adj ledodes-7.path
ln -sf ledodes-8.fa ledodes-scaffolds.fa
PathOverlap --overlap --dot -k87 ledodes-7.adj ledodes-7.path > ledodes-8.dot

2.4 使用转录子序列进行 rescaffolding

bwa index ledodes-8.fa
bwa mem -a -t2 -S -P -k87 ledodes-8.fa transcripts.fasta \
  |gzip > long1-8.sam.gz
abyss-longseqdist -k87 long1-8.sam.gz \
  |grep -v "l=" >long1-8.dist.dot
abyss-scaffold -k87 -s200 -n1 -g ledodes-8.path.dot ledodes-8.dot long1-8.dist.dot > ledodes-8.path
PathConsensus -k87 -p0.9 -s ledodes-9.fa -g ledodes-9.adj -o ledodes-9.path ledodes-8.fa ledodes-8.dot ledodes-8.path
cat ledodes-8.fa ledodes-9.fa \
  | MergeContigs -k87 -o ledodes-10.fa - ledodes-9.adj ledodes-9.path
ln -sf ledodes-10.fa ledodes-long-scaffs.fa

1. GCE 简介

GCE(Genome Characteristics Estimation) 是华大基因用于基因组评估的软件，其参考文献为：Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects。下载地址：ftp://ftp.genomics.org.cn/pub/gce。

GCE 软件包中主要包含 kmer_freq_hash 和 gce 两支程序。前者用于进行 kmer 的频数统计，后者在前者的结果上进行基因组大小的准确估算。

2. GCE 下载和安装

$ wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
$ tar zxf gce-1.0.0.tar.gz -C /opt/biosoft

3. kmer_freq_hash 的使用

kmer_freq_hash 的常用例子：

$ /opt/biosoft/gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash \
  -k 21 -l reads.list -a 10 -d 10 -t 24 -i 50000000 -o 0 -p species &> kmer_freq.log

kmer_freq_hash 的常用参数：

-k <17>
    设置 kmer 的大小。该值为 9~27，默认值为 17 。
-l string
    list文本文件，其中每行为一个fastq文件的路径。
-t int
    使用的线程数，默认为 1 。
-i int
    初始的 hash 表大小，默认为 1048576。该值最好设置为 （kmer 的种类数 / 0.75）/ 线程数。如果基因组大小为 100M，测序了 40M 个 reads，reads 的长度为 100bp，测序错误率为 1%，kmer的大小为 21，则kmer的种类数为100M+40M*100*1%*21=940M，若使用24线程，则该参数设置为 i=940M/0.75/24=52222222。
-p string
    设置输出文件的前缀。
-o int
    是否输出 k-mer 序列。1: yes, 0: no，默认为 1 。推荐选 0 以节约运行时间。
-q int
    设置fastq文件的phred格式，默认为 64。该值可以为 33 或 63。
-c double
    设置k-mer最小的精度，该值位于 0~0.99，或为 -1。 -1 表示不对 kmer进行过滤。设置较高的精度，可以用于过滤低质量 kmer。精度是由 phred 格式的碱基质量计算得来的。
-r int
    设置获取 k-mer 使用到的 reads 长度。默认使用 reads 的全长。
-a int
    忽略read前面该长度的碱基。
-d int
    忽略read后面该长度的碱基。
-g int
    设置使用该数目的碱基来获取 k-mers，默认是使用所有的碱基来获取 k-mer。

kmer_freq_hash 的主要结果文件为 species.freq.stat。该文件有 2 列：第1列是kmer重复的次数，第二列是kmer的种类数。该文件有255行，第225行表示kmer重复次数>=255的kmer的总的种类数。该文件作为 gce 的输入文件。
kmer_freq_hash 的输出到屏幕上的信息结果保存到文件 kmer_freq.log 文件中。该文件中有粗略估计基因组的大小。其中的 Kmer_individual_num 数据作为 gce 的输入参数。

4. gce 的使用

gce 的常用例子：

$ /opt/biosoft/gce-1.0.0/gce \
  -f species.freq.stat -c 85 -g 4112118028 -m 1 -D 8 -b 1 > species.table 2> species.log

gce 的常用参数：

-f string
    kmer depth frequency file
-c int
    kmer depth frequency 的主峰对应的 depth。gce 会在该值附近找主峰。
-g int
    总共的 kmer 数。一定要设定该值，否则 gce 会直接使用 -f 指定的文件计算 kmer 的总数。由于默认下该文件中最大的 depth 为 255，因此，软件自己计算的值比真实的值偏小。同时注意该值包含低覆盖度的 kmer。
-M int
    支持最大的 depth 值，默认为 256 。
-m int
    估算模型的选择，离散型（0），连续型（1）。默认为 0，对真实数据推荐选择 1 。
-D int
    precision of expect value，默认为 1。如果选择了 -m 1，推荐设置该值为 8。
-H int
    使用杂合模式（1），不使用杂合模式（0）。默认值为 0。只有明显存在杂合峰的时候，才选择该值为 1 。
-b int
    数据是（1）否（0）有 bias。当 K > 19时，需要设置 -b 1 。

gce 的结果文件为 species.table 和 species.log 。species.log 文件中的主要内容：

raw_peak	now_node	low_kmer	now_kmer	cvg	genome_size	a[1]	b[1]
84	35834245	22073804	4044916750	84.6637	4.83093e+07	0.928318	0.637648

raw_peak： 覆盖度为 84 的 kmer 的种类数最多，为主峰。
now_node： kmer的种类数。
low_kmer： 低覆盖度的 kmer 数。
now_kmer： 去除低覆盖度的 kmer 数，此值 = （-g 参数指定的总 kmer 数） - low_kmer 。
cvg：估算出的平均覆盖度
genome_size：基因组大小，该值 = now_kmer / cvg 。
a[1]： 在基因组上仅出现 1 次的 kmer 之 种类数比例。
b[1]： 在基因组上仅出现 1 次的 kmer 之 数量比例。该值代表着基因组上拷贝数为 1 的序列比例。

如果使用 -H 1 参数，则会得额外得到如下信息：

for hybrid: a[1/2]=0.223671 a1=0.49108
kmer-species heterozygous ratio is about 0.125918

上面结果中，0.125918 是由 a[1/2] 计算出来的。 0.125918 = a[1/2] / （ 2- a[1/2] ) 。
a[1/2]=0.223671 表示在所有的 uniqe kmer 中，有 0.223671 比例的 kmer 属于杂合 kmer 。

此外，有 a[1/2] 和 b[1/2] 的值在最后的统计结果中。重复序列的含量 = 1 - b[1/2] - b[1] 。

则杂合率 = 0.125918 / kmer_size 。 若计算出的杂合率低于 0.2%，个人认为测序数据应该是纯合的。这时候，应该不使用 -H 1 参数。使用 -H 1 参数会对基因组的大小和重复序列含量估算造成影响。

5. 不同杂合率，有无重复序列的 kmer species 和 kmer individuals 图

下图中 a 和 b 是对理想中无重复的基因组在不同杂合率下的曲线图；
下图中 c 和 d 是对有重复的基因组(human)在不同杂合率下的曲线图。
从下图可以参考不同杂合率下的曲线状况。

1. FGAP简介

FGAP利用BLAST将contigs序列比对到基因组草图序列上，寻找重叠到gap区间的最优序列，从而进行补洞。其参考文献：Piro, Vitor C., et al. “FGAP: an automated gap closing tool.” BMC research notes 7.1 (2014): 371.

2. FGAP下载和安装

FGAP官网：http://www.bioinfo.ufpr.br/fgap/。

$ wget http://sourceforge.net/projects/fgap/files/MCR_LINUX64b.tar.gz/download
$ tar zxf MCR_LINUX64b.tar.gz
$ cd MCR_LINUX64b
$ ./installMCR.sh /opt/biosoft/MCR

$ wget http://sourceforge.net/projects/fgap/files/FGAP_1_7_LINUX64b.tar.gz/download
$ tar zxf FGAP_1_7_LINUX64b.tar.gz -C /opt/biosoft/

2. FGAP的使用

FGAP的简单使用示例：

$ ln -s /opt/biosoft/FGAP_1_7_LINUX64b/sample_data/* .
$ /opt/biosoft/FGAP_1_7_LINUX64b/run_fgap.sh /opt/biosoft/MCR/v717/ \
 -d DRAFT_ecoli_hiseq454.fasta \
 -a "DATASET_ecoli_454.fasta,DATASET_ecoli_hiseq.fasta"

FGAP的使用参数：

-d /--draft-file	Draft genome file [fasta format - Ex: "draft.fasta"]
-a /--datasets-files	List of datasets files to close gaps [fasta format - Ex: "dataset1.fasta,dataset2.fasta"]

-s /--min-score		Min Score (raw) to return results from BLAST (integer) - Default: 25
-e /--max-evalue	Max E-Value to return results from BLAST (float) - Default: 1e-7
-i /--min-identity	Min identity (%) to return results from BLAST (integer [0-100]) - Default: 70

-C /--contig-end-length	Length (bp) of contig ends to perform BLAST alignment (integer) - Default: 300
-T /--edge-trim-length	Length of ignored bases (bp) upstream and downstrem of the gap (integer) - Default: 0
-R /--max-remove-length	Max number of bases (bp) that can be removed (integer) - Default: 500
-I /--max-insert-length	Max number of bases (bp) that can be inserted (integer) - Default: 500

-p /--positive-gap	Enable closing of positive gaps (with insertion) (integer [0-1]) - Default: 1
-z /--zero-gap		Enable closing of zero gaps (without insert any base) (integer [0-1]) - Default: 0
-g /--negative-gap	Enable closing of negative gaps (overlapping contig ends) (integer [0-1]) - Default: 0

-c /--gap-char				Base that represents the gap (char) - Default: "N"
-b /--blast-path			Blast+ package path (only makeblastdb and blastn are needed, version 2.2.28+ or higher) - Default: ""
-l /--blast-alignment-parameters	BLAST alignment parameters (opengap,extendgap,match,mismatch,wordsize) - Default: "1,1,1,-3,15"
-r /--blast-max-results			Max results from BLAST for each query (integer) - Default: 200
-t /--threads				Number of threads (integer) - Default: 1

-m /--more-output	More output files with gap regions after and before gap closing (integer [0-1]) - Default: 0
-o /--output-prefix	Output prefix [File or folder - Ex: "out" or "out/" ] - Default: "output_fgap"
-h /--help		This help message

1. ABACAS简介

ABACAS 用于在有参考序列的情况下，对contigs序列进行连接。软件官网：http://abacas.sourceforge.net/

2. ABACAS 下载和安装

ABACAS是一支perl程序，其运行需要安装有MUMmer。

$ wget http://sourceforge.net/projects/abacas/files/abacas.1.3.1.pl
$ chmod 755 abacas.1.3.1.pl
$ mv abacas.1.3.1.pl ~/bin/
$ ln -s ~/bin/abacas.1.3.1.pl ~/bin/abacas.pl
$ perl -p -i -e 's#/usr/local/bin/perl#/usr/bin/perl#' ~/bin/abacas.1.3.1.pl

3. ABACAS 的使用

ABACAS 的简单例子：

$ abacas.pl -r ref.fasta -q query.fasta -p nucmer -o out_prefix

ABACAS 的使用参数：

-r string
    输入参考序列。参考序列为fasta格式并仅有1条序列。
-q string
    输入query序列。query序列为fasta格式。
-p nucmer|promer
    该参数可以设为nucmer或promer，表示使用相应的程序进行contig ordering。
-o string
    设定输出文件的前缀。默认为query.fasta_ref.fastq。
-c
    参考序列是一个环。
-m
    使用该参数后，则将定位并定向的query序列按顺序打印到文件 out_prefix.MULTIFASTA.fa 中。
-b
    使用该参数后，则将未能定位的 congtig 序列内容打印到文件 ou_prefix.contigsInbin.fas 中。
-N
    使用该参数后，则将由定位定向的conitg序列连接而成的pseudomolecule的序列打印到文件 out_prefix.NoNs.fasta 中。默认下，输出文件 out_prefix.fasta 中的 pseudomolecule 序列中 congtig 之间的 gap 使用 N 表示，有重叠则加 100个N；该参数则是生成额外的fasta文件，其中contig之间序列无N。
-t
    使用该参数后，则对 ou_prefix.contigsInbin.fas 中的序列运行 blastall，额外生成 blast.out 文件。使用该参数则必须有 -b 参数。
-g out_prefix
    使用该参数后，则将 pseudomolecule 的gap部分在ref上的序列答应到文件 out_prefix.Gaps_onRef 中。
-P
    使用该参数后，则设计用于 close gaps 的引物。需要安装 prime3_core，但我个人加该参数运行失败。
-R
    使用该参数，则不需要运行 mumer，适合于多次运行程序，节约时间。但我个人加该参数运行失败。
-e
    不进行 contig ordering，有利于直接进行引物设计。
-i int
    允许的最小的identity的百分比。默认为 40。
-v int
    允许的最小的 contig 覆盖率。默认为 40。
-V 0|1
    contig比对到ref序列上的位置的数目。默认为1。使用 -V 0 表示将contig序列随机放到1个比对上的位点。
-l int
    不使用低于此长度的contig序列。默认为 100。 
-d
    ABACAS在调用MUmer的时候，默认下使用了--maxmatch参数来增加匹配，使用 -d 参数则会使用MUmer的默认参数，即不会使用--maxmatch参数。

1. REAPR 简介

REAPR(Recognition of Errors in Assemblies using Paired Reads)能利用成对的reads来识别基因组序列中的错误。从而，能将基因组序列从错误的gap处断开或将错误序列使用 Ns 代替。同时，对错误信息进行统计。

2. REAPR下载和安装

REAPR官网：http://www.sanger.ac.uk/resources/software/reapr/
安装 REAPR 需要先安装 R 和 Perl 模块： File::Basename, File::Copy, File::Spec, File::Spec::Link, Getopt::Long, List::Util。

$ wget ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.tar.gz
$ tar zxf Reapr_1.0.17.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/Reapr_1.0.17
$ wget ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.manual.pdf
$ ./install.sh
$ echo 'PATH=$PATH:/opt/biosoft/Reapr_1.0.17' >> ~/.bashrc
$ source ~/.bashrc

测试Reapr能否正常运行：
$ wget ftp://ftp.sanger.ac.uk/pub4/resources/software/reapr/Reapr_1.0.17.test_data.tar.gz
$ tar zxf Reapr_1.0.17.test_data.tar.gz 
$ ./test.sh

3. REAPR 使用

reapr 最少需要输入基因组fasta文件和mate-pair数据（需要将RF方向转变成FR方向）文件。常用的运行步骤如下：

使用 facheck 命令检查基因组fasta文件，以防fasta文件序列名含有怪异字符。
$ reapr facheck genome.fasta
使用下列命令对fasta文件进行修改。
$ reapr facheck genome.fasta new_genome.fasta

使用SMALT将 mate-pair reads比对到基因组上。
$ reapr smaltmap genome.fasta jumping.1.fastq jumping.2.fastq jumping.bam

运行reapr pipeline对基因组进行评估。将结果输出到 outdir 文件夹中。
$ reapr pipeline genome.fasta jumping.bam outdir

reapr 也支持输入paired-end数据(可选)。通过将paired-end数据uniquely并perfecly比对到基因组，从而能计算出基因组上error-free的碱基位点并进行统计。该项目分析对error calling没有影响。值得注意的是： 默认设置下call一个error-free位点需要至少5x的覆盖度；同时输入的paired-end数据的方向要为inward(一般paired-end数据方向即为inward)，此外reapr仅支持inward方向，因此输入的jumping数据也要先转换为inward方向。

使用paired-end数据的命令行：

$ reapr smaltmap genome.fasta jumping.1.fastq jumping.2.fastq jumping.bam
$ reapr perfectmap genome.fasta frag.1.fastq frag.2.fastq insert_size perfect_prefix
$ reapr pipeline genome.fasta jumping.bam outdir perfect_prefix

reapr的运行流程：

4. reapr的结果

REAPR 会对基因组每个位点的 FCD(fragment coverage distribution)进行分析。期望的FCD和实际上的上FCD之差为”FCD error”。 FCD error往往代表不正确的scaffold连接、大的插入缺失或contig序列中不正确的连接。为了能更好地计算FCD，特别是gap位点的FCD，mate-pair数据的insert size越大越好。

REAPR对输入的paired reads数据进行检测，并将结果放入到文件夹 outdir/00.Sample/ 下。

gc_vs_cov.lowess.pdf: GC偏好性检测图
insert.in.pdf: FR(inner/inward)类型插入片段的长度检测图
insert.out.pdf: RF类型插入片段的长度检测图
insert.stats.txt: 插入片段长度统计

从错误位点打断后的基因组序列： outdir/04.break.broken assembly.fa 。
如果FCD error区有gap，则从gap处将基因组序列打断；若FCD error区没有gap，则将其用Ns代替。同时将被代替的序列写入到文件 outdir/04.break.broken assembly bin.fa中。

基因组错误统计文件： outdir/05.summary.report.txt。该文件统计位于 contig 中的 FCD error 的数目，和含有 gap 的 FCD error 的数目。