基因组Repeat Sequence的寻找

发表于由
回复

1. 重复序列简介和相关软件

参考自文献：Review：Identifying repeats and transposable elements in sequenced genomes: how to find your way through the dense forest of programs。E Lerat。Heredity (2010) 104, 520–533; doi:10.1038/hdy.2009.165; published online 25 November 2009

1.1 Repeats的分类

基因组中的repeats依据其序列特征分成2类：串联重复(tandem repeats) 和 散在分布在基因组中的重复序列(interspersed repeats).其中第二类主要是transposable elements(TEs).

第一类串联重复包含：microsatellites 或 simple sequence repeats（1-6个碱基为一个重复单元） 和 minisatellites（10-60个碱基的长序列为一个重复单元）.

TEs包含2种类型：class-I TEs通过RNA介导的(copy and paste）机制进行转座；class-II TEs通过DNA介导的(cut and paste）机制来转座. 前者称为retroelements，后者称为DNA transposons。

class-I TEs中主要由LTR(long terminal repeat)构成。LTR的部分序列可能具有编码功能。而non-LTR则包含2个子类：LINEs(long interspersed nuclear elements)和SINEs(short interspersed elements），其中前者可能具有编码功能，后者则没有。

class-II TEs中加入了一个子类 MITEs(miniature inverted repeat transposable elements),基于DNA的转座因子，但是确通过”copy and paste”的机制来转座(Wicker et al., 2007)。

1.2 鉴定重复序列的软件

对于不同的重复序列，需要使用不同的软件来进行鉴定。而鉴定的方法可以分为：基于library，基于重复序列的特定结构 或 重头预测。文献中给出的软件很多：http://www.nature.com/hdy/journal/v104/n6/fig_tab/hdy2009165t1.html#figure-title。

1.2.1 基于library-based的软件

library-based的软件，需要构建library，该library中包含很多来自不同物种某一重复序列的一致性序列，然后通过相似性比对来鉴定repeats。这种方法能对所有的种类的重复序列进行鉴定。此方法最经典最流行的软件是RepeatMasker；此方法中CENSOR和MASKERAID两个软件可以用于改良RepeatMasker的结果；此外，用于基因组的重复序列鉴定的还有GREEDIER（Li et al.,2008），该软件在其文章中表明该软件性能还不错，在repeats鉴定的敏感性上稍微比RepeatMasker高一点，但是repeats的鉴定率只有RepeatMasker的一半左右.

1.2.2 基于signature的软件

基于signature的方法主要用于TEs的鉴定。
1. 鉴定LTR逆转座子: LTR_STRUC (McCarthy and McDonald, 2003), LTR_PAR (Kalyanaraman and Aluru, 2006), FIND_LTR (Rho et al., 2007), RETROTECTOR (Sperber et al., 2007), LTR_FINDER (Xu and Wang, 2007) and LTRHARVEST (Ellinghaus et al., 2008)。文献中，这些软件中LTR_STRUC的敏感性最高(96%)，但是LTR的鉴定率只有30%；而LTRharvest的鉴定率最高(42%),敏感性67%.总体上，作者依次推荐的软件是LTRHARVEST和FIND_LTR(敏感性83%，鉴定率37%)。
2. 鉴定non-LTR retrotransposons: TSDFINDER (Szak et al., 2002), SINEDR (Tu et al., 2004) and RTANALYZER (Lucier et al., 2007)。其中，第一个软件用于验证RepeatMasker检测到的L1 insertions；第二个软件用于检测侧翼有TSDs（target site duplications 当重复序列插入到基因组上时，其两侧会带入短核酸序列的重复）的SINEs；第三个软件通过一些特征，比如TSDs，polyA尾和5’端核酸内切酶位点等来通过打分来检测L1逆转座子。
3. 鉴定MITEs:FINDMITE (Tu, 2001), TRANSPO (Santiago et al., 2002), MITE Analysis Kit (MAK) (Yang and Hall, 2003) and MITE Uncovering SysTem (MUST) (Chen et al., 2009)。文献中作者使用第一个软件报错，使用第三个软件却下载不到。第二个软件不能寻找新的MITEs，看来最好是使用最新的第四个软件。
4. 鉴定helitrons: HELITRONFINDER。该软件(Du et al., 2008)用来寻找玉米基因组中的helitrons（在动物和植物中有发现）。

1.2.3 重头预测的软件

1.2.3.1 自我比较的方法

通过BLAST、PALS等方法，将序列和自身进行比较，从而找出重复序列。软件有：REPEAT PATTERN TOOLKIT (Agarwal and States, 1994), RECON (Bao and Eddy, 2002), PILER (Edgar and Myers, 2005) and the BLASTER suite (used in Quesneville et al., 2005).其中RECON软件使用最广泛。

1.2.3.2 k-mer and spaced seed approaches

一定长度的k-mer出现了多次，可以鉴定为重复序列；spaced seed则是k-mer的一种衍生，seed上允许有一定的差异。

软件有：REPUTER (Kurtz and Schleiermacher, 1999), VMATCH (Kurtz, unpublished), REPEAT-MATCH (Delcher et al., 1999), MER-ENGINE (Healy et al., 2003), FORREPEATS (Lefebvre et al., 2003), REAS (Li et al., 2005), REPEATSCOUT (Price et al., 2005), RAP (Campagna et al., 2005), REPSEEK (Achaz et al., 2007), TALLYMER (Kurtz et al., 2008) and P-CLOUDS (Gu et al., 2008).

1.2.4 其它重复序列鉴定软件

其它一些用于鉴定非TEs的软件：Tandem Repeats Finder (TRF) (Benson, 1999), Tandem Repeat Occurrence Locator (TROLL) (Castelo et al., 2002), MREPS (Kolpakov et al., 2003), TRAP (Sobreira et al., 2006) and Optimized Moving Window Spectral Analysis (OMWSA) (Du et al., 2007) have been developed specifically to detect tandem repeats. The Inverted Repeat Finder (IRF) program (Warburton et al., 2004) was designed to search for inverted repeats.

1.3 多个软件整合的pipeline程序

REPEATMODELER pipeline (Smit, unpublished http://www.repeatmasker.org/RepeatModeler.html) includes the programs RECON, REPEATSCOUT, REPEATMASKER and TRF. It uses the output of the RECON and REPEATSCOUT programs to build, refine and classify consensus models of putative interspersed repeats.

当然，在此文献中，也讲述了其它很多专门用途的其它pipeline软件。而REPEATMODELER pipeline是现在运用于基因组的重复序列鉴定最主流的软件。以下将讲述该软件运用。

2. RepeatModeler的安装与使用

2.1 软件的安装

RepeatMasker和RepeatModeler是ISB（Institute for Systems Biology)的软件。ISB is located in the South Lake Union neighborhood。

根据RepeatMasker说明，其安装与使用需要：Perl 5.8.0以上版本，序列比对Engine，TRF和Repeat Database。其中序列比对engine可以安装多个，但每次只能使用其一，可以使用Cross_match,RMBlast,HMMER和ABBlaast/WUBlast等。

根据RepeatMideker说明，其安装与使用需要：Perl 5.8.8以上版本，RepeatMasker,Repeat Database,RECON,RepeatScout,TRF和序列比对engine。其中序列比对engine可以安装多个，但每次只能使用其一，可以使用RMBlast和ABBlaast/WUBlast。

再安装完毕需要的软件后，对RepeatMasker和RepeatModeler进行configure的时候填入相应软件的路径即可安装。

2.2 RepeatModeler的使用

2.2.1 使用RepeatModeler来通过基因组序列构建library

$ $RepeatModelerHome/BuildDatabase -name species \
  -engine ncbi species.genome.fasta
$ $RepeatModelerHome/RepeatModeler -database species

结果生成了一个文件夹，名称为RM_[PID].[DATE] ie. “RM_5098.MonMar141305172005″。该文件夹中的”consensi.fa.classified”即为library，用于RepeatMasker的输入。值得注意的是，可能fasta文件需要序列间不能有空（空格和换行等），否则会程序出错。

2.2.2 使用RepeatMasker来进行重复序列掩盖和重复序列计算

cp RM_5098.MonMar141305172005/consensi.fa.classified .
mkdir Repeat_result
$ $RepeatMaskerHome/RepeatMasker -pa 8 \
  -e ncbi -lib consensi.fa.classified \
  -dir Repeat_result -gff species.genome.fasta

则生成的结果文件位于Repeat_result文件夹中。

基因组草图的gap closer软件：GapFiller

发表于由
1

1. GapFiller简介

组装出来的基因组草图的scaffold需要进一步进行gaps的关闭。进行这样功能的软件有：SOAPdenovo GapCloser v1.12r6; IMAGE; GapFiller.

GapFiller文章发表在Genome Biology上：Boetzer M,Pirovano W. 2012. Toward almost closed genomes with GapFiller. Genome Biol.13:R56。从文章可以完全明白该软件closing gap的原理。

GapFiller需要输入scaffold序列(FASTA)和NGS paired-read数据(FASTA or FASTAQ)，输出FASTA格式文件。该软件的获得需要填写一些邮箱和单位信息。商业license需要花钱；学术性需要引用其文章。

2. GapFiller安装

下载GapFiller的安装包，解压缩后，里面包含bowtie、bwa和example共3个文件，其最重要的是GapFiller.pl文件，为主程序。还有2个PDF格式的manual文件。

3. GapFiller的使用

直接运行主程序，会给出软件的参数说明，如下：

-l library文件
-s scaffold序列的fasta文件
-m default:29 和gap边缘重叠的最小碱基数,该数值最好设置比reads的长度小一点点的数。比如36bp长度的reads，设置该值为30~35.
-o default:2 在补洞时，延伸一个碱基最小需要的reads数.
-r default:0.7 在补洞时，至少有该比例reads的碱基一致，才能对该碱基位点进行延伸。
-d default:50 gap部分序列的允许的最大差异。填补gap后，若值“填补上的序列长度 - gap长度”大于该阈值，则停止补洞；若小于该阈值，则不进行融合。
-n default:10 在一个scaffold中对邻近的两个contigs进行融合所需要最小重叠的碱基数。
-t default:10 由于gap边缘的碱基大部分是低质量碱基，补洞时需要先将gap边缘该数目的碱基trim掉，作为N处理。
-i default:10 迭代的最大次数。
-g default:1 使用bowtie进行比对的时候允许的最大的gap数，和bowtie中的-v参数一致
-T default:1 运行时使用的线程数
-S 跳过重新读取输入文件
-b 输出文件的basename。

-l 参数所指向的library文件需要先行编辑好。该文件包含7列，每一列之间以空格(space)隔开.其例子和格式如下：

Lib1 bwa file1.1.fasta file1.2.fasta 400 0.25 FR 
Lib1 bowtie file2.1.fasta file2.2.fasta 400 0.25 FR 
Lib2 bowtie file3.1.fastq file3.2.fastq 4000 0.5 RF

第1列：library名称
第2列：使用的序列比对方法，如果reads长度<50,则使用bowtie；若长度>50并<150,则使用bwa；若长度很大，比如454的reads，则使用bwa。BWA和BWA-sw运行在默认模式下。 第3,4列：双末端测序的fastq文件或fasta文件。 第5,6列：插入片段的长度，以及承认的长度。比如上例子中插入片段长度为400bp，成对的reads的片段长度只有在[400-400*0.25,400+400*0.25]范围内才被承认。 第7列：双端测序reads的方向，有FF，FR，RF和RR几种。

4. 例子

编辑一个libraries.txt文件，内容如下：

Illumina_160bp bwa fragment.reads1.fastq fragment.reads2.fastq 156 0.25 FR
Illumina_6000bp bwa jumping.reads1.fastq jumping.reads2.fastq 6170 0.25 FR

运行GapFiller程序，如下：

$GapFillerHome/GapFiller.pl -l libraries.txt\
  -s genome.fasta -m 90 -T 8 -b species

USB安装64位的CentOS6.4

发表于由
回复

一、方法1

1. 官网下载系统安装包

2. 将DVD1包烧录至U盘

在windows系统下可以使用Ultraiso

3. 替换vesamenu.c32文件

此步骤为关键的步骤，若不进行，则会在U盘引导计算机时，启动后卡在下面的界面：Press the key to begin the installation process。

使用CentOS6.0版本的vesamenu.c32文件替换U盘中syslinux/目录下的vesamenu.c32文件，即可。

4. 镜像文件

64位的镜像文件较大，需要将DVD1文件复制到一个可以容纳该文件的文件系统中。使用Ultraiso烧录的U盘为fat文件系统，不能装入镜像文件；常用的windows下NTFS文件系统不能被centos识别。推荐另外找快容量较大的U盘或移动硬盘，格式换成linux系统识别的ext4文件系统。然后将DVD1文件复制到ext4文件系统的磁盘根目录下；同时提取iso文件中的image文件夹到该磁盘根目录下。

此步骤也关键，否则，“定位在 Hard drive 上, directory holding image”此关就过不去了。

参考自：http://blog.chinaunix.net/uid-22566367-id-3346333.html，可以在该站点下载到 vesamenu.c32.zip

此方法本人使用来一次，发现安装完毕后没能引导进入系统。推荐第二种方法，一次性成功。

二、方法2

方法来自：如何利用Usb盘安装CentOS操作系统20111208。

1. 准备好CentOS镜像和容量足够的U盘，8G以上为宜

2. 在Linux上对U盘进行分区和文件系统创建

以下方法将U盘分成两个区：第一个区大小150M，用于fat32文件系统；第二个分区大小为剩下的所有部分，用于ext2文件系统。

# fdisk /dev/sdb
m    帮助命令
d    删除磁盘分区，如果有已经存在多个分区，多次该命令
n    新建分区
p    选择分区为主分区
1    第一个分区
  回车选择分区起始为默认的磁柱1开始
+150M    选择第一个分区大小为150M
t    分区的文件系统类型
l    列出所有的文件系统类型
b    选择文件系统为fat32

n    新建分区
p    选择分区为主分区
2    第二个分区
  回车选择分区起始为默认的磁柱开始
  回车选择分区结束区为最后的磁柱末

a    选择启动分区
1    以第一分区为启动分区
w    写入分区格式

分区完后，对分区进行格式化：

# mkfs.vfat -n BOOT /dev/sdb1
# mkfs.ext2 -m 0 -b 4096 -L DATA /dev/sdb2

以上第2个命令可能需要等待很长时间。

3. 挂在U盘分区并拷入文件

3.1 在默认情况下，重新插上U盘，则会自动将U盘分区挂载到了/media目录下。若默认情况下不挂载，则手动挂载如下：

# mkdir /media/BOOT /media/DATA
# mount /dev/sdb1 /media/BOOT
# mount /dev/sdb2 /media/DATA
# mount -ro loop CentOS-6.4-x86_64-bin-DVD1.iso /mnt

以上命令中也将CentOS6.4的镜像文件挂载到了/mnt中。

接下来需要将镜像文件中的isolinux文件夹拷贝到U盘的第一分区中，并稍作修改，并让第一分区作为启动分区，用来启动CentOS的安装界面；将镜像文件拷贝到U盘第二分区，支持整个系统的安装，同时将镜像中的images/install.img拷贝到U盘第二分区，用于系统的安装。

# cp CentOS-6.4-x86_64-bin-DVD1.iso /media/DATA/
# mkdir /media/DATA/images
# cp /mnt/images/install.img /media/DATA/images/

# yum install syslinux
# syslinux -s /dev/sdb1    配置syslinux在U盘第一分区
# dd if=/usr/share/syslinux/mbr.bin of=/dev/sdb    写入MBR

# cp /mnt/isolinux/ /media/BOOT/syslinux
# mv /media/BOOT/syslinux/isolinux.cfg /media/BOOT/syslinux/syslinux.cfg
# rm -f /media/BOOT/syslinux/isolinux.bin

# umount /mnt
# umount /media/BOOT
# umount /media/DATA

至此，CentOS的U盘启动安装盘作好了，可以用于CentOS系统的安装。不需要像第一个方法中变更vesamenu.c32文件，这样可能是导致安装完毕后无法正常引导系统的原因。

需要注意的是，我是将第二分区内的镜像文件和images文件夹一起放在了centos文件夹下面；这样按如何利用Usb盘安装CentOS操作系统20111208中所述，则需要修改/media/BOOT/syslinux/syslinux.cfg文件,在文件中加入下述的斜体字部分。好处是让第二分区在进行正常文件存放时更整洁。

append initrd=initrd.img repo=hd:sdb2:/centos

同时，在安装过程中我卡在了“定位在 Hard drive 上, directory holding image”这一关卡上。首先硬盘驱动里面只有/dev/sda,没有/dev/sdb。于是我换了个USB插口，就能显示了(可能是USB口原因，或需要重新插拔U盘)。然后在下面的横线中出现的”/images/install.img”前面手动加上“/centos”，于是正确定位到了该文件（默认情况下不需要加“/centos”）。

NGS数据的Error correction方法

发表于由
回复

参考自文献：Yang X,Chockalingam SP,Aluru S. A survey of error-correction methods for next-generation sequencing. Brief. Bioinformatics 2013;14:56-66.。

现在进行error-correciton的算法有三种： k-spectrum-based、Suffix tree/array based 和 MSA based。本文对以下软件使用真实数据进行了评估：HSHREC、Reptile、Quake、SOAPec、HiTEC、ECHO、Coral（只有第一个和最后一个，共2个软件用于454和ion数据的评估；所有软件都用于了Illumina数据的评估）。

该文献的结果表明：适用于Illumina数据的软件(我个人认为的优先顺序)是 Reptile、HiTEC(在2*100的数据中结果最好) 和 ECHO。Reptile的参数选择比其它两个软件要麻烦很多；HiTEC不适合处理有‘N’的或不同长度reads。适用于454和ion数据的软件是Coral。

以上是对基因组的测序reads进行修正的方法，而新出了一个对RNA-Seq数据进行修正的软件：SEECER。其文献为：Probabilistic error correction for RNA sequencing。文中结果表示，以上对基因组reads进行修正的方法对RNA-Seq reads的修正效果不好。

以下是这几个软件的安装与使用方法：

1. HiTEC

HiTEC的下载网页：http://www.csd.uwo.ca/~ilie/HiTEC/。根据其安装源码包内的说明文件，HiTEC的安装需要先行安装gsl library,然后修改Makefile文件，最后进行make。

HiTEC的参数比较简单，运行方法如下：

$ ./hitec <inputReads> <correctedReads> \
  <GenomeLength> <per-base error rate>

$ ./hitec S.aureusReads.fasta \
S.aureusCorrectedReads.fasta 2820462 1

输入fasta或fastq文件，给定一个大致的基因组大小和碱基错误率，则会运行，并将结果输出到指定文件中。

HiTEC将忽略那些含有碱基 “N” 的reads，这些reds不会在输出文件中显示。输出结果为fasta文件，fasta的头被更换了，变成了从0开始的数字。

2. ECHO

ECHO的下载网页：http://uc-echo.sourceforge.net/。根据其安装源码包内的说明文件，ECHO的安装需要有GCC，Python，SciPy和NumPy，使用yum安装即可。ECHO的运行也比较简单：

$ python ErrorCorrection.py -b 2000000 --nh 256 \
  --ncpu 6 -o output/5Mreads.fastq 5Mreads.fastq

生成了结果文件和HiTEC一样，生成fasta文件，fasta的头被更换了，变成了从0开始的数字。

3. reptile

reptile的下载网页：http://aluru-sun.ece.iastate.edu/doku.php?id=reptile。

reptile的使用说明相比前两个软件，更加完整。其使用说明参考软件包中的readme文件和网页的turorial。以下为reptile使用实例：

3.1 将fastq文件转换成fa文件和q文件

对象为名为reads.fastq的文件，位于当前工作目录。reads.fastq为使用NGSQC-toolkit过滤掉低质量reads和含有N的reads后的文件。

$ $reptileHome/utils/fastq-converter-v2.0.pl ./ ./ 1

结果生成reads.fa和reads.q文件。这两个文件都为fast格式，fasta头为从0开始按顺序排列的数字，内容分别为序列和质量。

3.2 seq-analy配置文件参数调整

copy一个范本的seq-analy配置文件和Reptile配置文件。

$ cp $reptileHome/utils/seq-analy/config-example seq-analy-config
$ cp $reptileHome/src/config-example reptile-config

将该配置文件修改如下：

IFlag                   1
BatchSize               1000000
InFaFile                ./reads.fa
IQFile                  ./reads.q
KmerLen                 24
OKmerHistFile           ./kmerhist
QHistFile               ./qualhist
OKmerCntFile    
MaxErrRate              0.02
QThreshold              73
MaxBadQPerKmer          4
QFlag                   1

修改好输入和输出文件路径后，运行程序：

$ $reptileHome/utils/seq-analy/seq-analy seq-analy-config

以上命令执行后生成根据配置文件信息输出./kmerhist和./qualhist两个文件。

继续修改seq-analy-config配置文件信息：
1. reads长度为35bp的时候，MaxBadQPerKmer的值设为4-6；reads长度为100+bp的时候，MaxBadQPerKmer的值设为6-10.
2. KmerLen可以按自己的意愿设置，可以保留为默认的24.
3. 查看./qualhist文件信息，找出第三列累计频率 >=0.8 的那一行对应的碱基质量对应的ASCII值，作为QThreshold的值；同时找出第三列累计频率 <= 0.2 的那一行对应的碱基质量对应的ASCII值，作为下一步 Reptile配置文件中 Qlb 参数的值。 修改完seq-analy-config配置文件后，继续运行“seq-analy seq-analy-config”，直到不用继续修改为止；如果当前不用修改，则进入下一步。

3.3 Reptile配置文件参数调整

对范本Reptile配置文件修改，修改读取文件和输出文件，如下：

InFaFile                ./reads.fa
QFlag                   1
IQFile                  ./reads.q
OErrFile                ./reads.errors
BatchSize               1000000 
KmerLen                 12
hd_max                  1
Step                    12
ExpectSearch            16 
T_ratio                 0.5

######## The following parameters need to be tuned to specific dataset ########

QThreshold              73
MaxBadQPerKmer          4
Qlb                     62
T_expGoodCnt            16
T_card                  6

通时，对上数的一些参数进行修改：
1. T_expGoodCnt，设置为./kmerhist第三列累计频率最接近0.95的一行的第一列的值的2倍；
2. T_card，设置为./kmerhist第二列中第二大的值的一行所对应的第一列的值；
3. KmerLen，设置为 >= (seq-analy最终的配置文件中该参数的值的一半);
4. QThreshold, 和seq-analy最终的配置文件中该值的设置一致；
5. Qlb， 设置为./qualhist文件中，第三列累计频率 <= 0.2 的那一行对应的碱基质量对应的ASCII值； 6. MaxBadQPerKmer，和seq-analy最终的配置文件中该值的设置一致； 7. step 设置为【（seq-analy最终的配置文件中KmerLen值） - (本reptile参数文件中KmerLen的值)】的值

3.4 运行Reptile

$ $reptileHome/src/reptile-v1.1 reptile-config

生成了配置文件中设定的结果文件./reads.errors。该文件格式为：

ReadID ErrNum [pos to from qual] [pos to from qual] …

第一列是read ID; 第二列为该read错误的碱基数; 对每一个错误的碱基，有4个值给出来：该错误碱基的位置(pos)、新的修正后的碱基(to)、原来的错误的碱基(from) 和 该位点碱基质量对应的ASCII值(qual)。

3.5 生成修正后的fasta文件

$ $reptileHome/utils/reptile_merger/reptile_merger \
  ./reads.fa ./reads.errors out.fa

Velvet的安装和使用

发表于由
2

1. Velvet的安装

Velvet用于基因组的de novo组装，支持各种原始数据，包括Illumina的short reads和454的long reads。

首先下载velvet的安装包，直接使用make命令来编译，即可获得可执行主程序velveth和velvetg。安装如下：

$ wget http://www.ebi.ac.uk/~zerbino/velvet/velvet_1.2.10.tgz
$ tar zxf velvet_1.2.10.tgz
$ cd velvet_1.2.10
$ make 'CATEGORIES=10' 'MAXKMERLENGTH=57' 'LONGSEQUENCES=1' \
  'OPENMP=1' 'BUNDLEDZLIB=1'

值得注意的是，make后有多种参数，需要对velvet软件进行需要的设置：

CATEGORIES=10： 默认情况下velvet只能输入 2 groups of short reads，此处设置为10. 如果有多个文库或多种类型的原始数据，需要相应增加该值的大小。该值越大，耗内存越大。

MAXKMERLENGTH=57： 最大的Kmer长度，默认情况下该值为 31， 一般de novo组装基因组，31 是不够的，故需要增大该值。组装过程中kmer越大，越耗内存。

BIGASSEMBLY=1： 当超过 2.2G 的reads用于组装基因组的时候，需要设置该值。实际上很少会有如此之多的reads用于基因组组装，可以不用设置该值。设置该值后，会消耗更多的内存。

LONGSEQUENCES=1： 当组装出的contigs长度超过 32kb 长的时候，需要设置该值。会消耗跟多内存。

OPENMP=1：打开多线程运行。需要设置环境变量 OMP_NUM_THREADS 和 OMP_THREAD_LIMIT。 Velvet最多使用 OMP_NUM_THREADS+1 或 OMP_THREAD_LIMIT 个线程。也值哟部分的velvet算法支持多线程运行，不会造成运行时间的线性减少。

BUNDLEDZLIB=1： 默认velvet使用系统自带的zlib，如果系统没有zlib，则需要加入该参数来使用velvet源码包中自带的zlib。

在上述 make 编译完 Velvet 后，继续velvet的使用

2. 使用velveth来准备数据

velveth接受输入的文件，产生一个hash表。生成两个文件：Sequences和Roadmaps。velveth用法为：

$ velveth
$ velveth output_directory hash_length \
  -file_format -read_type filename1 filename2 ... \
  [-file_format -read_type filename1 filename2 ...]

在不带参数情况下，直接运行velveth会给出其帮助文件。而其参数如下：

output_directory：输出文件夹的路径

hash_length： 设置Kmer的大小。该值3点要求：1.必须为奇数；2.必须小于或等于编译velvet时设置的MAXKMERLENGTH值；3.必须小于reads的长度。

file_format： 支持的格式有fasta（默认）、fastq、bam等。

reads_type： short(默认)、shortpaired、short2、shortpaired2 … short10、shortpaired10、long、longPaired。 默认情况下short reads只保留了2个通道。在之前设置中将CATEGORIES值设置为10，则velvet最多支持10种不同类型的short reads数据。long支持长的reads，比如sanger测序数据和454测序数据。 如果reads_type一致，则同一个reads_type下可以有多个文件。

filename1： reads的文件名。如果是成对的reads，则需要将两个reads文件合并成一个文件。Velvet安装文件夹中提供了这样的一个文件。

一个velveth的例子。对5个Illumina测序的结果和一个454测序的结果使用velvet进行组装：

$ velveth output/ 31 \
  -shortPaired -fastq fragment1.reads.fastq \
  -shortPaired2 -fastq fragment2.reads.fastq \
  -shortPaired3 -fastq jumping1.reads.fastq \
  -shortPaired4 -fastq jumping2.reads.fastq \
  -long -fastq 454.fasta

3. 使用velvetg来进行基因组组装

velvetg是vlevet软件的进行de Bruijin图构建和操作的核心。在命令行下直接键入命令velvetg，这样描述：velvetg – de Bruijn graph construction, error removal and repeat resolution。其使用方法如下：

$ velvetg
$ velvetg directory [options]

$ velvetg output -exp_cov auto -cov_cutoff auto \
  -shortMatePaired3 yes -shortMatePaired4 yes \
  -clean yes -scaffolding yes -amos_file yes

velvetg的具体参数如下：

directory 工作的目录名，即为上一步骤velveth中的输出文件夹。

-cov_cutoff <float|auto> default： no removal
设置最低kmer覆盖度的值。默认下会生成很多nodes，而有些nodes很短，覆盖度较低，需要去除掉。auto则自动设置该值，是该值为所有nodes的kmer覆盖度值的median值的1/2。

-exp_cov <float|auto> default： no long or paired-end read resolution。
期望的kmer覆盖度。如果设置了auto，则该值为所有nodes的kmer覆盖度值的median值; 该值设置为auto，则同时自动设置-cov_cutoff为auto。如果对杂合基因组进行组装时，设置auto，却很难进行预测，组装结果肯定不好。auto适用于标准的基因组测序。

-long_cov_cutoff <float> default： no removal
移除低long-read覆盖度的nodes。

-ins_length <interger>
成对的short reads间的distance的期望值，即插入片段的平均长度。

-ins_length* <interger>
*代表第*组shortPaired reads, 和 velveth 中的参数相对应。

-ins_length_long <interger>
和前两个参数一样，代表longPaired reads的插入片段长度。

-ins_length*_sd <interger>
此时’*’代表’nothing、2、3…、long’等，和上面的3个参数相匹配；该值设置插入片段长度的标准差。有关插入片段长度和sd的如果不设置，velvet则会自动计算。velvet是将成对的reads比对到组装出来的nodes上，最后计算出一个insert size 和sd。这样做的话，可能估算的不准确，需要看velvet的运行输出信息，检测其预算是否准确。

-scaffolding <yes|no> defautl: yes
是否要使用paired end信息进行scaffolds组装

-max_branch_length <integer> default: 100
处理气泡(bubble)的参数。默认下，如果两条序列超过了100bp的位点处的kmer不一致，则将这两条序列分开成单独的contigs。

-max_divergence <float> default: 0.2
在一个bubble内两条branches最大的差异率（分母为较长的branch的长度）。

-max_gap_count <integer> default: 3
在一个bubble内两条branches比对结果中，运行的最大gap数。该参数和上两个参数为重要的参数，能很大程度上影响组装结果。如果设置得松点，可分别将值设为200,0.33,5。

-min_pair_count <integer> default: 5
默认，将两个contigs连成scaffold最少需要5对paired reads的验证。如果测序深度较大，则可以提高该值；测序深度低，则降低该值。

-long_mult_cutoff <integer> default: 2
默认下，融合两个contigs需要最少有2个long reads验证。

-max_coverage <float> default: no removal
最高的覆盖度，高于此覆盖度的nodes将被删除。

-coverage_mask <integer> default: 1
contigs最小的置信覆盖度，低于此覆盖度的contigs被maksed。

-shortMatePaired* <yes|no> default: no
如果哪一个shortPaired为mate-pair library测序的结果，则需要指定该参数为yes。

-conserveLong <yes|no> default: no
是否保留含有long reads的序列

-unused_reads <yes|no> default: no
是否输出unused reads, 保存到 UnusedReads.fa 中。

-alignments <yes|no> default: no
是否输出contig比对到参考序列的summary.

-exportFiltered <yes|no> default: no
是否输出由于覆盖度的原因被过滤掉的long nodes。

-clean <yes|no> default: no
是否删除所有的不能用于重新计算的中间文件

-very_clean <yes|no> default: no
是否删除所有的中间文件（删除后不能重新计算）

-min_contig_lgth <integer> default: hash length * 2
输出结果中最小的contigs长度

-amos_file <yes|no> default: no
是否输出AMOS文件

velvetg的默认输出结果

directory/contigs.fa 长度2倍长于kmer的contigs。参数-scaffolding决定生成的该fasta文件是否包含scaffold序列。
directory/stats.txt 用于决定覆盖度cutoff的统计表
directory/PreGraph 初始的de vruijin图
directory/Graph2 最终的de bruijin图 关于该文件中内容的解释，请见velvet PDF manual。
directory/velvet_asm.afg AMOS兼容的组装文件，能用于AMOS基因组组装软件包
directory/Log velvet的运行记录

4. Velvet提供了额外的一些scripts

这些scripts非常有用，位于$velvetHome/contrib/文件夹下。

4.1 estimate-exp_cov

在使用velvetg组装出一个初步结果后，利用结果文件stats.txt来估算出kmer的期望覆盖度。

$ $velvetHome/contrib/estimate-exp_cov/velvet-estimate-exp_cov.pl \
  [options] stats.txt

4.2 observed-insert-length.pl

在使用velvetg组装出一个初步结果后，以结果文件夹为输入，计算insert size和insert sd。

$ $velvetHome/contrib/observed-insert-length.pl/observed-insert-length.pl \
  [options] Velvet_directory

4.3 shuffleSequences

将对称的reads1和reads2文件合并成一个文件，用于velvet的输入。velvet不能将两个文件用于输入。

$ $velvetHome/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  reads1.fastq reads2.fastq reads.fastq

4.4 VelvetOptimiser

输入该脚本的命令，不加参数，给出详细的帮助信息。该程序用于选取不同的Kmer值，来对原始数据进行组装，同时计算出最优化的组装参数。最后给出最优化方案的组装结果。非常适合于简单基因组的自动化组装操作。其参数设置较为简单。当然需要能在$PATH中有velveth和velvetg两个主要的velvet命令。以下就不详细介绍该命令参数，仅给出一个例子:

$ $velvetHome/contrib/VelvetOptimiser-2.2.4/VelvetOptimiser.pl \
  --s 17 --e 97 --x 2 --a --t 4 --k 'n50' --c 'Lbp' \
  --f '-shortPaired -fastq fragment1.reads.fastq \
  -shortPaired2 -fastq fragment2.reads.fastq \
  -shortPaired3 -fastq jumping1.reads.fastq \
  -shortPaired4 -fastq jumping2.reads.fastq \
  -long -fastq 454.fasta' \
  -o 'shortMatePaired3 yes shortMatePaired4 yes' \
  [-g 40M]

以上命令意义为：使用kmer长度从17到97，每次运行kmer长度+2； –a表示要生成amosfile；使用N50来决定kmer的选取；以’Lbp’,即大于1kb的contigs的总碱基数来决定coverage cutoff的值的选取；–f后为velveth的参数；–o后接velvetg的参数；-g后为预计的基因组大小，用于计算该命令下内存的消耗，然后退出运行，由于该程序会同时运行多个velvet和velveg，消耗的内存巨大。

4.5 AssemblyAssembler

该程序能使用不同的kmer组装出很多Assemblies，然后将这些Assemblies组装融合出一个最终的结果。其使用方法位于该程序的通目录下的README文本文件中。给出一个例子如下：

$ $velvetHome/contrib/AssemblyAssembler1.3/AssemblyAssembler1.3.py \
  -s 17 -e 99 -v $velvetHome -c 5.1 -t 35 -i 300 -m 51 -r y \
  -f "-fastq -shortPaired reads.fastq"

以上命令意义为：使用从17到99的kmer，每种kmer都组装出基因组并最后合并成一个结果。-c为coverage cutoff值；-t设置为只有大于35的kmer时，才进行coverage cutoff；-i表示插入片段长度为300；-m表示期望的kmer覆盖度为51；-r表示结果给出short reads要包含到最后的组装结果中。

可以看出该程序适合于只有一种类型的short reads时，才适合。因为其参数-i只能设置一种插入片段长度大小。

4.6 其它程序

extractContigReads用于提取和指定conting相对应的reads；

afg_handing 用于检查 .afg 文件；

fasta2agp 用于将fasta格式的基因组组装结果(gaps用N表示)转换成AGP文件，用于提交到EMBL或NCBI；

show_repeats 指出Assembly中larger repeated contigs的长度；

read_prepare 和 select_paired 用于NGS数据的过滤

5. 有参考基因组下velvet的组装

转录组De novo组装软件

发表于由
1

1. Trinity

Trinity，本博客之前介绍了：Trinity的安装与使用

Citation: Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011 May 15;29(7):644-52. doi: 10.1038/nbt.1883. PubMed PMID: 21572440.

2. Oases

Oases是EMBL-EBI开发出来的，是一种使用short reads（such as Illumina, SOLiD, or 454）进行de novo转录组组装的软件。Oases是在velvet的基础上开发的，其运行需要安装Velvet。

Citation：M.H. Schulz, D.R. Zerbino, M. Vingron and Ewan Birney. Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics, 2012. DOI: 10.1093/bioinformatics/bts094.

上文献中写道：It was tested on human and mouse RNA-seq data and is shown to improve significantly on the transABySS and Trinity de novo transcriptome assemblers.

3. EBARDenovo

EBARDenovo, which stands for Extension, Bridging And Repeat-sensing Denovo. This software package (with patent pending) is free of charge for academic use only.

Citation：Hsueh-Ting Chu,William W. L. Hsiao,Jen-Chih Chen,Tze-Jung Yeh,Mong-Hsun Tsai,Han Lin,Yen-Wenn Liu,Sheng-An Lee,Chaur-Chin Chen,Theresa T. H. Tsao,and Cheng-Yan Kao。EBARDenovo: highly accurate de novo assembly of RNA-Seq with efficient chimera-detection Bioinformatics (2013) 29 (8): 1004-1010 first published online March 1, 2013 doi:10.1093/bioinformatics/btt092

上文献中写道：our algorithm is the most accurate among the examined programs, including de Bruijn graph assemblers, Trinity and Oases.

4. Trans-ABySS

Trans-ABySS: Analyze ABySS multi-k-assembled shotgun transcriptome data.

Trans-ABySS 1.4.4, Released Oct 09, 2012, Supports both transcriptome and genome assemblies. A more robust pipeline and improvement to the 3 analysis suites (fusions, indels, and splicing).
Citation:Robertson G,et al. De novo assembly and analysis of RNA-seq data. Nat. Methods 2010;7:909-912.

Variant calling 软件

发表于由
回复

1. SNP和INDEL caller

1.1 GATK

1.2 SAMtools

1.3 DISCOVAR

DISCOVAR的官方博客：http://www.broadinstitute.org/software/discovar/blog/。该软件和ALLPATHS-LG是同一个组织开发的，是在ALLPATHS-LG的基础上开发的，用于基因组的组装和Variants的发掘。2013.3 公布出了该软件。

1.4 PyroHMMvar

PyroHMMvar: a sensitive and accurate method to call short INDELs and SNPs for Ion Torrent and 454 data. 2013.8.28发表在bioinformatics上。

Illumina数据模拟软件

发表于由
1

现有的Illumina数据模拟软件有以下几种：

1. pIRS

piRS: Profile-based Illumina pair-end reads simulator. 是BGI开发出来的软件。使用C++和Perl编写，下载地址：ftp://ftp.genomics.org.cn/pub/pIRS/。pIRS同时提供了一个引入variants的工具，从而能模拟杂合基因组数据。

文献：Xuesong Hu,Jianying Yuan,Yujian Shi,Jianliang Lu,Binghang Liu,Zhenyu Li,Yanxiang Chen,Desheng Mu,Hao Zhang,Nan Li,Zhen Yue,Fan Bai,Heng Li,and Wei Fan。 pIRS: Profile-based Illumina pair-end reads simulator Bioinformatics (2012) 28 (11): 1533-1535
first published online April 15, 2012 doi:10.1093/bioinformatics/bts187

2. ART

ART能模拟Illumina, 454和SOLiD的这3中NGS数据，包括 single-end, paired-end和mate-pair三种类型的数据。ART是NIH开发出的软件。其官网网址为：http://www.niehs.nih.gov/research/resources/software/biostatistics/art/。现在ART支持Illumina MiSeq的数据模拟了。该软件更新较多较快。

Citation: Weichun Huang, Leping Li, Jason R Myers, and Gabor T Marth. ART: a next-generation sequencing read simulator, Bioinformatics (2012) 28 (4): 593-594

3. 其它软件

Wgsim from the Samtools package (Li et al., 2009)

MetaSim (Richter et al., 2008) for simulating metagenomic data

Mason (http://seqan.de/projects/mason.html) for both Illumina and 454 reads

SimSeq (https://github.com/jstjohn/SimSeq) for Illumina reads

FlowSim (Balzer et al., 2010) for 454 reads.

PBSIM: PacBio reads simulator–toward accurate genome assembly

计算paired-end测序的insert size

发表于由
回复

给出的条件为：
参考序列文件： ref.fasta
Illumnina paired-end文件： reads1.fq, reads2.fq

1. 结合 Bowtie2，SAMtools 和 Picard的使用来计算insert size

1.1 将paired-end数据比对到参考序列上

$ bowtie2-build ref.fasta ref
$ bowtie2 --rg-id librarylength_insert --rg "PL:ILLUMINA" \
  --rg "SM:Strain" -x ref -p 8 \
  -1 reads1.fq -2 reads2.fq -S align.sam

1.2 将sam格式比对结果转换为bam格式并排序

使用samtools将sam文件转换成bam文件；使用picard将bam文件进行排序.

$ samtools view -bS align.sam > align.bam
$ jar Xmx2g -jar $picardHome/SortSam.jar \
  I=align.bam O=align.sort.bam SO=coordinate

1.3 计算insertsize

输入排序后的bam文件，使用picard生成insertsize的txt的结果文件和PDF格式的图形结果。

$ jar Xmx2g -jar $picardHome/CollectInsertSizeMetrics.jar \
  I=align.sort.bam R=ref.fasta  \
  O=insertsize.txt H=insertsize.pdf

2. 使用Qualimap来计算insertsize

2.1 Qualimap的介绍

Qualimap 用于NGS比对数据的质量控制。使用JAVA和R编写的程序，有GUI和command-line两种运行方式。

文献：Fernando García-Alcalde, Konstantin Okonechnikov, José Carbonell, Luis M. Cruz, Stefan Götz, Sonia Tarazona, Joaquín Dopazo, Thomas F. Meyer, and Ana Conesa “Qualimap: evaluating next-generation sequencing alignment data.” Bioinformatics 28, no. 20 (2012): 2678-2679.
doi: 10.1093/bioinformatics/bts503

2.2 Qualimap的安装

Qualimap的正常使用需要 JAVA runtime vesion 6 或以上版本； R 2.14 或以上版本；需要一些 R 包。该软件默认使用最高线程数运行，速度不错。

$ wget http://qualimap.bioinfo.cipf.es/release/qualimap_v0.7.1.zip
$ unzip qualimap_v0.7.1.zip
$ cd qualimap_v0.7.1
$ sudo su
$ /usr/local/bin/Rscript scripts/installDependencies.r
$ R
> source("http://bioconductor.org/biocLite.R")
> biocLite("gplots","Repitools")

3. 注意事项

对于基因组的paired-end测序，将reads比对到参考基因组上；但是对于转录组测序的paired-end测序，则最好将reads比对到转录组序列上。而比对到基因组序列上reads数会偏少。

结合GATK和SAMtools从头挖掘SNPs和INDELs

发表于由
16

提出问题

给出以下的已有结果：
1. 物种species的基因组fasta文件: genome.fasta
2. 对该物种的一个名为sample的基因组DNA样，进行了Illumina hiseq2000 paired-end测序(300bp插入片段文库，~20x基因组碱基覆盖度)，结果文件： reads1.fq, reads2.fq
问怎么进行SNPs和INDELs的calling？

文章有很多地方引自Nowind的博文：GATK使用简介 Part2/2
流程参照：http://www.broadinstitute.org/gatk/guide/topic?name=best-practices

以下给出具体的答案和步骤, 以上述已有结果的3个文件所在的目录为当前工作目录，各种软件的主目录以$software表示。多线程的程序以8个线程运行。

1. Raw reads的处理

使用NGSQC toolkit使用默认设置对raw reads进行过滤。

$ perl $NGSQCToolkitHome/QC/IlluQC_PRLL.pl \
  -pe reads1.fq reads2.fq 2 5 -c 8 \
  -o QC

2. 将过滤后的reads比对到基因组上

$ bowtie2-build genome.fasta genome
$ bowtie2 --rg-id sample --rg "PL:ILLUMINA" --rg "SM:sample" \
  -x geome -1 ./QC/reads1.fq_filtered -2 ./QC/reads2.fq.filtered \
  -p 8 -S sample.sam

3. 将sam文件转换为Bam文件并标记出PCR duplicates

$ samtools view -bS sample.sam > sample.bam
$ java -Xmx10g -jar $picardHome/SortSam.jar \
  INPUT=sample.bam OUTPUT=sample.sort.bam \
  SORT_ORDER=coordinate
$ java -Xmx10g -jar $picardHome/MarDuplicates.jar \
  INPUT=sample.sort.bam OUTPUT=sample.dd.bam \
  METRICS_FILE=sample.dd.metrics \
  MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

4. 构建索引文件

$ samtools faidx genome.fasta
$ java -Xmx10g -jar $picardHome/CreateSequenceDictionary.jar \
  R=genome.fasta O=genome.dict
$ samtools index sample.dd.bam

5. 对INDEL周围进行重新比对(realignment)

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T RealignerTargetCreator \
  -I sample.dd.bam -o sample.realn.intervals
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T IndelRealigner \
  -targetIntervals sample.realn.intervals \
  -I sample.dd.bam -o sample.realn.bam

6. 第1遍Variant calling

6.1 使用GATK和samtools进行SNP和INDEL calling

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T UnifiedGenotyper \
  -I sample.realn.bam -o sample.gatk.raw1.vcf \
  --read_filter BadCigar -glm BOTH \
  -stand_call_conf 30.0 -stand_emit_conf 0
$ samtools mpileup -DSugf genome.fasta sample.realn.bam | \
  bcftools view -Ncvg - > sample.samtools.raw1.vcf

6.2 对Variant结果进行综合和过滤

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T SelectVariants \
  --variant sample.gatk.raw1.vcf \
  --concordance sample.samtools.raw1.vcf \
  -o sample.concordance.raw1.vcf
$ MEANQUAL=`awk '{ if ($1 !~ /#/) { total += $6; count++ } } \
  END { print total/count }' sample.concordance.raw1.vcf`
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T VariantFiltration \
  --filterExpression "QD < 20.0 || ReadPosRankSum < -8.0 || \
  FS > 10.0 || QUAL < $MEANQUAL" \
  --filterName LowQualFilter --variant sample.concordance.raw1.vcf \
  --missingValuesInExpressionsShouldEvaluateAsFailing
  --logging_level ERROR -o sample.concordance.flt1.vcf
$ grep -v "Filter" sample.concordance.flt1.vcf > \
  sample.concordance.filter1.vcf

7. 对第5步获得的realign的bam文件的进行校正

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T BaseRecalibrator \
  -I sample.realn.bam -o sample.recal_data.grp \
  -knownSites sample.concordance.filter1.vcf
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T PrintReads \
  -I sample.realn.bam -o sample.recal.bam \
  -BQSR sample.recal_data.grp

8. 第2遍Variant calling

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T UnifiedGenotyper \
  -I sample.recal.bam -o sample.gatk.raw2.vcf \
  --read_filter BadCigar -glm BOTH \
  -stand_call_conf 30.0 -stand_emit_conf 0
$ samtools mpileup -DSugf genome.fasta sample.recal.bam | \
  bcftools view -Ncvg - > sample.samtools.raw2.vcf
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T SelectVariants \
  --variant sample.gatk.raw2.vcf \
  --concordance sample.samtools.raw2.vcf \
  -o sample.concordance.raw2.vcf
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T VariantFiltration \
  --filterExpression "QD < 10.0 || ReadPosRankSum < -8.0 || \
  FS > 10.0 || QUAL < 30" \
  --filterName LowQualFilter --variant sample.concordance.raw2.vcf \
  --missingValuesInExpressionsShouldEvaluateAsFailing
  --logging_level ERROR -o sample.concordance.flt2.vcf
$ grep -v "Filter" sample.concordance.flt2.vcf > \
  sample.concordance.filter2.vcf
$ cp sample.concordance.filter2.vcf sample.final.vcf

最后的结果文件为sample.final.vcf。关于具体的vcf的格式详解见：http://www.hzaumycology.com/chenlianfu_blog/?p=1514