Assessing Quality of Error Correction

Question

The main idea:

Given fastq/fasta files containing reads(simulated from ref. genome, with errors introduced) and another file with the ground truth(simulated but without errors) and then the resulting fastq/fasta files from error correction program ABC, is there an existing program that will give some kind of score and analysis on how well some error correction program ABC did.

In depth:

I have 2 tools I am using to synthetically generate reads from Illumina machines. They both take in a reference genome and then produce 2 fastq files that are paired end reads with errors. One of these programs produces a second pair of fastq files that are the same reads but without errors(ground truth). The other program produces a bam/sam file for ground truth.

The fastq files with errors are passed to error correction programs which then produce a new set of fasta files. I am looking to compare the quality of the error correction between these programs, as well as my own method. I am not bias towards one method or the other and could generate alignments or sequences from alignments. I am just looking for a program that exists that is scalable and fast at assesing the quality of the ouputs of error correction programs. If there are multiple programs I have to chain together, that would work too, I just don't want to write my own if one already exists.

Answer 1

Partial answer, but I want to discuss something that is only rarely mentioned in any read simulation / ground truth studies.

Given that you're working with Illumina reads, you need to be very careful with evaluating that ground truth and identifying repetitive sequences. In particular, don't assume that a particular high-complexity sequence will have a unique match to the ground-truth genome; you need to either ignore repetitive sequences (which would lead to biases), or have some way of incorporating repetition into your model.

As a demonstration of the issue, I did a 100bp repeat survey of the Telomere-to-Telomere CHM13 assembly a few years ago:

https://zenodo.org/doi/10.5281/zenodo.4762023

One of these plots (from chr1) is shown here:

Without getting too deep into the details, each point on this plot represents a location where there is a sequence of at least 100bp that is found identically in at least two different locations in the chromosome. Because most of these locations are sparsely-scattered high-complexity sequences, many repeat-finder algorithms (e.g. see here) will ignore them. Any program that maps sequences to these locations and tries to evaluate ground truth as "sequence X must be in location Y" will produce false results.

Bear in mind that this is only looking for identical sequences; slight variations (e.g. 1 SNP differences) are not represented on this plot.