NGS数据的Error correction方法

发表评论
2,371

A+

所属分类：

Genomics

现在进行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为使用

NGS

QC-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

参考自文献：

Yang X,Chockalingam SP,Aluru S. A survey of error-correction methods for next-generation sequencing. Brief. Bioinformatics 2013;14:56-66.

原文来自：http://www.hzaumycology.com/chenlianfu_blog/?p=1702

原文链接：https://blog.csdn.net/u010608296/article/details/102398125