NGS数据的Error correction方法

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的参数比较简单，运行方法如下：

1. $ ./hitec <inputReads> <correctedReads> \
2. <GenomeLength> <per-base error rate>
3.  
4. $ ./hitec S.aureusReads.fasta \
5. 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的运行也比较简单：

1. $ python ErrorCorrection.py -b 2000000 --nh 256 \
2. --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后的文件。

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

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

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

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

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

将该配置文件修改如下：

1. IFlag 1
2. BatchSize 1000000
3. InFaFile ./reads.fa
4. IQFile ./reads.q
5. KmerLen 24
6. OKmerHistFile ./kmerhist
7. QHistFile ./qualhist
8. OKmerCntFile
9. MaxErrRate 0.02
10. QThreshold 73
11. MaxBadQPerKmer 4
12. QFlag 1

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

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配置文件修改，修改读取文件和输出文件，如下：

1. InFaFile ./reads.fa
2. QFlag 1
3. IQFile ./reads.q
4. OErrFile ./reads.errors
5. BatchSize 1000000
6. KmerLen 12
7. hd_max 1
8. Step 12
9. ExpectSearch 16
10. T_ratio 0.5
11.  
12. ######## The following parameters need to be tuned to specific dataset ########
13.  
14. QThreshold 73
15. MaxBadQPerKmer 4
16. Qlb 62
17. T_expGoodCnt 16
18. 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

1. $ $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文件

1. $ $reptileHome/utils/reptile_merger/reptile_merger \
2. ./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