使用DBG2OLC对二、三代混合数据进行基因组组装

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 (<>) {print; $_ = <>; print; <>; <>;}' Pacbio_Cell01.fastq >> ctg_pb.fasta
$ perl -e 'while (<>) {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格式简单说明

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)

使用 abyss 进行 scaffolding

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

使用 GCE 进行基因组大小评估

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)在不同杂合率下的曲线图。
从下图可以参考不同杂合率下的曲线状况。
kerm_pictures

使用FGAP进行补洞

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

使用ABACAS进行有参考基因组的contigs连接

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参数。

使用 REAPR 进行基因组错误评估

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的运行流程:
reapr_pipeline

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 的数目。

使用SOPRA进行scaffold连接

1. SOPRA 简介

SOPRA(Statistical Optimization of Paired Read Assembly),主要用于利用成对的reads信息来进行scaffold连接,其准确性非产高。其参考文献:Dayarian, Adel, Todd P. Michael, and Anirvan M. Sengupta. “SOPRA: Scaffolding algorithm for paired reads via statistical optimization.” BMC bioinformatics 11.1 (2010): 345.

2. SOPRA的下载和安装

$ wget http://www.physics.rutgers.edu/~anirvans/SOPRA/SOPRA_v1.4.6.zip
$ mkdir /opt/biosoft/SOPRA_v1.4.6/
$ unzip SOPRA_v1.4.6.zip -d /opt/biosoft/SOPRA_v1.4.6/
安装目录/opt/biosoft/SOPRA_v1.4.6/下有SOPRA的manual文档
$ chmod 755 /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/*.pl
$ perl -p -i -e 's\s*$/\n/' /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/*.pl
SOPRA的主要运行程序就是这4支perl程序。

3. SOPRA的使用流程

SOPRA 可以支持 Illumina 和 SOLID 的二代测序数据,但此处仅讲解使用 Illumina 的测序数据。SOPRA 也支持同时利用多个paired测序文库的数据,同时支持paired-end和mate-paired数据。

3.1 准备输入数据

不管是paired-end或mate-paired数据,首先需要将数据的方向转换为FR方向。即需要将mate-paired的2端reads都进行反向互补。然后,需要将fatstq文件2端的reads文件合并为一个文件,并转换为fasta文件作为SOPRA的输入数据。
此外,SOPRA还需要contig序列作为输入的基因组序列。
以1个paired-end数据和1个mate-paired数据为例:

先将 jumping 文库的数据进行反向互补
$ fastq_rc.pl jump.1.fastq > 11; mv 11 jump.1.fastq
$ fastq_rc.pl jump.2.fastq > 22; mv 22 jump.2.fastq

然后,合并两端的reads序列
$ /opt/biosoft/velvet_1.2.10/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  frag.1.fastq frag.2.fastq frag.fastq
$ /opt/biosoft/velvet_1.2.10/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  jump.1.fastq jump.1.fastq jump.fastq

再将fastq转换为fasta
$ perl -e '$num = 1; while (<>) { print ">$num\n"; $_=<>; print; $_=<>; $_=<>; $num ++; }' frag.fastq > frag.fasta
$ perl -e '$num = 1; while (<>) { print ">$num\n"; $_=<>; print; $_=<>; $_=<>; $num ++; }' jump.fastq > jump.fasta

3.2 SOPRA 对contigs序列和测序数据进行格式处理

使用s_prep_contigAseq_v1.4.6.pl程序对基因组的contig序列以及illunia测序数据进行序列名的修改,以利于后续程序运行。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_prep_contigAseq_v1.4.6.pl \
  -contig contig.fasta -mate frag.fasta jump.fasta -a SOPRA_OUT/

程序参数:
-contig string
    输入基因组 contigs 序列。
-mate string
    输入 illumina 数据。如果有多个数据,则多个fasta文件之间使用空格分开。
-a string
    设置输出文件夹目录。

在 SOPRA_OUT 目录中生成了 contigs_sopra.fasta、frag_sopra.fasta 和 jump_sopra.fasta 文件。查看这3个文件内容,发现是对fasta文件中的序列名进行了修改。

3.3 将Illumina数据比对到基因组

可以利用bowite,bwa等序列比对软件将Illumina数据按单端序列比对到基因组序列上。若一条序列比对到多个位置,推荐报告多个结果。因为后续的s_parse_sam_v1.4.6.pl程序会去除掉比对到多个位置的reads。若仅报告最优结果,我不知道程序会如何处理。要是不去除这样的比对结果,可能会对scaffolding的准确性有影响。因此,还是要注意报告多个比对结果,特别是bowtie2默认下仅报告最优结果。
此外,去除reads尾部碱基质量低的碱基后,illumina reads 的匹配率若是提高很多,则推荐去除reads尾部10~20个碱基,再用于比对。由于mate-pair数据进行了反向互补,则其原本的尾部成了首部。

$ cd SOPRA_OUT/
$ bowtie2-build contigs_sopra.fasta contig
$ bowtie2 -p 20 -x contig -k 10 -f -3 15 -U frag_sopra.fasta -S frag_sopra.sam
$ bowtie2 -p 20 -x contig -k 10 -f -5 15 -U jump_sopra.fasta -S jump_sopra.sam

程序参数:
-k 10
    若read比对到多个位点,则最多报告 10 个结果。
-3 15
    去除read尾部15bp碱基后再用于比对。
-5 15
    去除read首部15bp碱基后再用于比对。

3.4 对SAM文件进行分析

使用s_parse_sam_v1.4.6.pl对sam文件进行分析并简化其信息。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_parse_sam_v1.4.6.pl \
  -sam frag_sopra.sam jump_sopra.sam -a ./

程序参数:
-sam string
    输入 sam 文件。多个illumina 文库的sam文件用空格隔开。
-a string
    设置输出文件的路径。

得到sam文件分析后的结果frag_sopra.sam_parsed和jump_sopra.sam_parsed。

3.5 进行序列方向和距离分析

读取上一步简化后的sam文件信息,进行分析。得到contigs序列的覆盖度、library文库的插入片段长度,序列长度等。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_read_parsed_sam_v1.4.6.pl \
  -parsed frag_sopra.sam_parsed -d 500 -parsed jump_sopra.sam_parsed -d 3000 -a ./

程序参数:
-parsed string
    输入 .sam_parsed 文件。若有多个文件,则用空格分割。
-d int
    illumina数据的插入片段长度。
-c int default: 5
    如果read及其反向互补序列在library中出现的次数>=该值,则不计算其成对信息。
    在运行 s_scaf_v1.4.6.pl 程序后,输入日志中得到如下信息:
    Average number of links between two contigs using minlength 150 and minlink 2 is 63.38, ...
    Starting cycle 1 of orientation assignment
    Average number of links between two contigs using minlength 150 and minlink 4 is 184.84, ...
    如果该值(184.84)大于100,则需要考虑降低 -c 参数的值。
    此外,如果数据量比较少,平均覆盖度较低,比如低于 10,则考虑增加 -c 参数的值。

-e int default:0
    如果设置其值为 1, 则使用输入的插入片段长度值。默认下软件会计算插入片段长度值,并使用软件计算出来的值。
-a string
    设置输出文件的路径。

输出文件为 orientdistinfo_c5 。 此外,在 div 文件夹中有contigs序列的覆盖度、library文库的插入片段长度,序列长度等统计信息。

3.6 进行scaffold连接

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_scaf_v1.4.6.pl \
  -o orientdistinfo_c5 -a ./

程序常用参数:
-o string
    输入文件 orientdistinfo_c。
-w int default: 4
    连接2个contigs所需要的最小的pair links。
    在运行 s_scaf_v1.4.6.pl 程序后,输入日志中得到如下信息:
    Average number of links between two contigs using minlength 150 and minlink 2 is 63.38, ...
    Starting cycle 1 of orientation assignment
    Average number of links between two contigs using minlength 150 and minlink 4 is 184.84, ...
    如果该值(184.84)大于100,则需要考虑提高 -w 参数的值。一般 -w 取该值的 4~5%。如果该值较低,比如低于10的时候,则需要降低 -w 参数的值。
-L int default: 150
    用于连接scaffold的最小长度的contigs序列。
-h float default: 2.2
    设置一个系数,用于排除对高覆盖度contigs的scaffold连接。大于 mean_coverage + h * std_coverage 该覆盖度的 contigs ,不能用于scaffold连接
-a string
    设置输出文件的路径。

程序输出文件为 scaffolds_h2.2_L150_w4.fasta 。该文件为最终的结果文件。

使用 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

只根据两个转录组数据分析差异表达基因

1. 数据

手头有两个转录组数据,分别是某一真菌物种(该物种没有基因组数据)的双核菌丝阶段和子实体阶段的转录组数据。测序平台是Illumina Hiseq2000;插入片段长度200bp,测序的reads长度90bp。得到的数据文件为:

/home/user/RNA-seq/mycelium_reads1.fastq
/home/user/RNA-seq/mycelium_reads2.fastq
/home/user/RNA-seq/fruitingbody_reads1.fastq
/home/user/RNA-seq/fruitingbody_reads2.fastq

2. 数据的预处理

2.1. 去除N的比例大于5%的reads;去除低质量reads(质量值Q≤20的碱基数占整个read的50%以上);
2.2. 根据FastQC对上述过滤后的reads的质量检测,去除reads首尾各10bp的碱基,得到的预处理数据为:

/home/user/RNA-seq/clean_reads/mycelium_reads1.fastq
/home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
/home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq
/home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

3. 使用Trinity和TGICL进行转录组的组装

3.1. 对mycelium数据和fruitingbody数据分别进行转录组的组装

$ pwd
/home/user/RNA-seq/

$ mkdir assembly
$ cd assembly

$ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output mycelium_contig --left /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq --right /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
$ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output fruitingbody_contig --left /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq --right /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

当然,需要将生成的两个转录组的Trinity.fasta序列按长度进行排序;统一更改的fasta
头名称;更改fasta文件名

$ cp mycelium_contig/Trinity.fasta mycelium_contigs.fasta
$ cp fruitingbody_contig/Trinity.fasta fruitingbody_contigs.fasta

3.2. 使用TGICL将两个转录组序列合并

$ pwd
/home/user/RNA-seq/assembly

$ mkdir all_tissue_contig
$ cd all_tissue_contig

$ cat ../mycelium_contigs.fasta ../fruitingbody_contigs.fasta > all.contigs.fasta
$ tgicl -F all.contigs.fasta

tgicl生成了两个有用的文件: asm_1/contigs  和 all.contigs.fasta.single
tons。其中前者是聚类后的contigs结果;后者是没有聚类的单独的contigs的序列名,需
要分别到../mycelium_contigs.fasta 和 ../fruitingbody_contigs.fasta文
件中提取出相应的序列: all.contigs.fasta.singletons.mycelium.fasta 和 
all.contigs.fasta.singletons.fruitingbody.fasta

$ cat asm_1/contigs all.contigs.fasta.singletons.mycelium.fasta all.contigs.fasta.singletons.fruitingbody.fasta > all_contigs.fasta

在对all_contigs.fasta进行fasta头的重命名,并将序列按长度排序;同时要得到all_
contigs.fasta中的序列和mycelium_contigs.fasta,fruitingbody_contigs.
fasta中序列的对应关系。

$ cp all_contigs.fasta ../

3.3. 至此得出转录组的组装结果:

/home/user/RNA-seq/assembly/all_contigs.fasta

4. 使用cufflinks来分析差异表达基因

4.1 使用tophat来将预处理后的reads比对到转录组序列上

$ pwd
/home/user/RNA-seq/

$ mkdir cufflinks
$ cd cufflinks

$ bowtie2-build ../all_contigs.fasta all_contigs

$ tophat -o tophat_out_mycelium -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
$ tophat -o tophat_out_fruitingbody -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

4.2 自制一个transcripts.gtf文件
由于cuffdiff运行需要一个参考序列的transcripts.gtf文件: 该文件有9列,使用tab分隔;使exon的范围为整个contig。其格式如:

AA.auricula_all_contig_1     chenlianfu  exon  1  9401  .  .  .  gene_id "A.auricula_all_contig_1"; transcript_id "A.auricula_all_contig_1";
A.auricula_all_contig_10     chenlianfu  exon  1  4464  .  .  .  gene_id "A.auricula_all_contig_10"; transcript_id "A.auricula_all_contig_10";
A.auricula_all_contig_100    chenlianfu  exon  1  3090  .  .  .  gene_id "A.auricula_all_contig_100"; transcript_id "A.auricula_all_contig_100";
A.auricula_all_contig_1000   chenlianfu  exon  1  1768  .  .  .  gene_id "A.auricula_all_contig_1000"; transcript_id "A.auricula_all_contig_1000";
A.auricula_all_contig_10000  chenlianfu  exon  1  586   .  .  .  gene_id "A.auricula_all_contig_10000"; transcript_id "A.auricula_all_contig_10000"; 
A.auricula_all_contig_10001  chenlianfu  exon  1  586   .  .  .  gene_id "A.auricula_all_contig_10001"; transcript_id "A.auricula_all_contig_10001";

4.3 使用cuffdiff来分析转录子的表达量和差异表达基因

$ pwd
/home/user/RNA-seq/cufflinks

$ cuffdiff -L mycelium,fruitingbody --library-type fr-unstranded -p 8 -o cuffdiff ./transcriptome.gtf ./tophat_out_mycelium/accepted_hits.bam ./tophat_out_fruitingbody/accepted_hits.bam

至此得出转录子的表达量数据和差异表达分析

/home/user/RNA-seq/cufflinks/cuffdiff/gene_exp.diff