Is there a way to filter fastq reads based on Qscore before paired end merging

[simmimourya]

I have forward and reverse FASTQ files from Illumina paired-end sequencing, and I'm merging these reads using vsearch, which works great. However, I need to filter the reads based on a quality score (Qscore) threshold of 30 before merging. Specifically, within each forward and reverse FASTQ file, I want to discard any reads that have a quality score below this threshold. Once the reads are filtered, I want to perform a paired-end merge on the filtered reads.

Did I miss something? Also, is this the right way to approach the problem?

Here’s the code I’m using and the error I encountered during filtering:

def filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold=30, merge_override=False):

    filtered_r1_path = os.path.join(output_dir, "filtered_r1.fastq")
    filtered_r2_path = os.path.join(output_dir, "filtered_r2.fastq")
    
    # Filter R1 and R2 reads
    subprocess.run(f"vsearch --fastq_filter {r1_path} --fastqout {filtered_r1_path} --fastq_qmin {qscore_threshold}", shell=True, check=True)
    subprocess.run(f"vsearch --fastq_filter {r2_path} --fastqout {filtered_r2_path} --fastq_qmin {qscore_threshold}", shell=True, check=True)
    
    # Merge filtered reads
    merged_output_prefix = os.path.join(output_dir, "merged")
    merged_output_file = f"{merged_output_prefix}.fastq"
    
    if not os.path.exists(merged_output_file) or merge_override:
        subprocess.run(f"vsearch --fastq_mergepairs {filtered_r1_path} --reverse {filtered_r2_path} --fastqout {merged_output_file}", shell=True, check=True)
    else:
        print(f"Using existing merged file: {merged_output_file}")
    
    return merged_output_file

# Example usage
r1_path = "forward.fastq.gz"
r2_path = "reverse.fastq.gz"
output_dir = '../results/'
filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold=30, merge_override=False)

This is the error I'm seeing upon running the script:

vsearch v2.28.1_linux_x86_64, 124.5GB RAM, 16 cores
https://github.com/torognes/vsearch

Reading input file

Fatal error: FASTQ quality value (16) below qmin (30)
---------------------------------------------------------------------------
CalledProcessError                        Traceback (most recent call last)
Cell In[35], line 36
     34 r2_path = "reverse.fastq.gz"
     35 output_dir = '../results/'
---> 36 filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold=30, merge_override=False)

Cell In[35], line 19, in filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold, merge_override)
     16 filtered_r2_path = os.path.join(output_dir, "filtered_r2.fastq")
     18 # Filter R1 and R2 reads
---> 19 subprocess.run(f"vsearch --fastq_filter {r1_path} --fastqout {filtered_r1_path} --fastq_qmin {qscore_threshold}", shell=True, check=True)
     20 subprocess.run(f"vsearch --fastq_filter {r2_path} --fastqout {filtered_r2_path} --fastq_qmin {qscore_threshold}", shell=True, check=True)
     22 # Merge filtered reads

File /opt/conda/lib/python3.9/subprocess.py:528, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
    526     retcode = process.poll()
    527     if check and retcode:
--> 528         raise CalledProcessError(retcode, process.args,
    529                                  output=stdout, stderr=stderr)
    530 return CompletedProcess(process.args, retcode, stdout, stderr)

CalledProcessError: Command 'vsearch --fastq_filter forward.fastq.gz --fastqout ../results/filtered_r1.fastq --fastq_qmin 30' returned non-zero exit status 1.

[frederic-mahe]

Hi @simmimourya

in your python script, the option --fastq_qmin 30 tells vsearch to stop with an error if a fastq file contains positions with Q<30. Hence the error you get when running it on your data.

With vsearch, you can filter reads based on length, number of unknown positions (Ns), expected error (EE), or the number of copies:

The sequences may be filtered using the options --fastq_maxee,
--fastq_maxee_rate, --fastq_maxlen, --fastq_maxns, --fastq_minlen
(default 1), --fastq_trunclen, --maxsize, and --min‐size. Sequences
not satisfying the requirements are discarded. For pairs of sequences,
both sequences in a pair must satisfy the requirements, otherwise both
are discarded.

As quality tends to decrease towards the 3' end of reads, the usual approach is to trim reads when the cumulated expected error (EE) drops below a certain threshold, or after the first position with a Q-value equal or lesser than a certain threshold. For example, you can trim reads (and read pairs) after the first low Q-value with --fastq_truncqual 30. Here is a toy-example:

vsearch \
    --fastx_filter <(printf "@fwd\nAA\n+\nI?\n") \
    --reverse <(printf "@rev\nTT\n+\nII\n") \
    --quiet \
    --fastq_truncqual 30 \
    --fastqout - \
    --fastqout_rev -

@fwd
A
+
I
@rev
TT
+
II

A Q-value of ? (30) is equal to our threshold, so fwd is truncated.

If you decide to go that way and to use --fastq_truncqual, you can use it directly when merging. It will truncate the reads before trying to merge them. Here is a toy-example:

vsearch \
    --fastq_mergepairs <(printf "@fwd\nAAAAAAAATA\n+\nIIIIIIIIIH\n") \
    --reverse <(printf "@rev\nTATTTTTTTT\n+\nIIIIIIIIII\n") \
    --fastq_truncqual 39 \
    --fastqout -

Pairs that failed merging due to various reasons:
         1  overlap too short

H (Q=39), so fwd is truncated, and is not long enough anymore to be merged (merging requires a minimal overlap of 10 nucleotides). With --fastq_truncqual 38, there is no trimming and reads are merged.

[simmimourya]

Thank you for your response. I was able to come up with something similar yesterday, glad to see that our approaches align. I saw documentation about truncqual here and came up with this.

def filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold=30, merge_override=False):
    filtered_r1_path = os.path.join(output_dir, "filtered_r1.fastq")
    filtered_r2_path = os.path.join(output_dir, "filtered_r2.fastq")
    # ---------------------- Code change -------------------------------------------
    # Filter R1 and R2 reads
    subprocess.run(f"vsearch --fastq_filter {r1_path} --fastqout {filtered_r1_path} --fastq_truncqual {qscore_threshold}", shell=True, check=True)
    subprocess.run(f"vsearch --fastq_filter {r2_path} --fastqout {filtered_r2_path} --fastq_truncqual {qscore_threshold}", shell=True, check=True)
    # ---------------------- Code change -------------------------------------------
    # Merge filtered reads
    merged_output_prefix = os.path.join(output_dir, "merged")
    merged_output_file = f"{merged_output_prefix}.fastq"
    if not os.path.exists(merged_output_file) or merge_override:
        subprocess.run(f"vsearch --fastq_mergepairs {filtered_r1_path} --reverse {filtered_r2_path} --fastqout {merged_output_file}", shell=True, check=True)
    else:
        print(f"Using existing merged file: {merged_output_file}")
    return merged_output_file
# Example usage
r1_path = "forward.fastq.gz"
r2_path = "reverse.fastq.gz"
output_dir = '../results/'
filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold=30, merge_override=False)

Which helped me filter the reads but I'm seeing this issue now:

375836 sequences kept (of which 138630 truncated), 10387 sequences discarded.
vsearch v2.28.1_linux_x86_64, 124.5GB RAM, 16 cores
https://github.com/torognes/vsearch
Reading input file 100%
360321 sequences kept (of which 176234 truncated), 25902 sequences discarded.
vsearch v2.28.1_linux_x86_64, 124.5GB RAM, 16 cores
https://github.com/torognes/vsearch
Merging reads
Fatal error: More forward reads than reverse reads
---------------------------------------------------------------------------
CalledProcessError                        Traceback (most recent call last)
Cell In[36], line 38
     36 r2_path = "reverse.fastq.gz"
     37 output_dir = '../results/'
---> 38 filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold=30, merge_override=False)
Cell In[36], line 29, in filter_and_merge_reads(r1_path, r2_path, output_dir, qscore_threshold, merge_override)
     26 merged_output_file = f"{merged_output_prefix}.fastq"
     28 if not os.path.exists(merged_output_file) or merge_override:
---> 29     subprocess.run(f"vsearch --fastq_mergepairs {filtered_r1_path} --reverse {filtered_r2_path} --fastqout {merged_output_file}", shell=True, check=True)
     30 else:
     31     print(f"Using existing merged file: {merged_output_file}")
File /opt/conda/lib/python3.9/subprocess.py:528, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
    526     retcode = process.poll()
    527     if check and retcode:
--> 528         raise CalledProcessError(retcode, process.args,
    529                                  output=stdout, stderr=stderr)
    530 return CompletedProcess(process.args, retcode, stdout, stderr)
CalledProcessError: Command 'vsearch --fastq_mergepairs ../results/filtered_r1.fastq --reverse ../results/filtered_r2.fastq --fastqout ../results/merged.fastq' returned non-zero exit status 1.

Filtering before merging can lead to issues due to unequal read lengths, which can prevent successful merging. Instead, I switched to merging first and then filtering, and this approach works really well! Now, I want to do a sanity check on some of the filtered sequences to assess metrics like Qscore distribution, cumulative percentage (AccPct), mean quality, etc. I’ve already generated the stats using --fastq-stats. What key metrics would you recommend examining to verify the quality of the FASTQ data after filtering and merging?

# generate stats
def fastq_stats(fastq_file, output_log, fastq_ascii='sanger', fastq_qmin=0, fastq_qmax=40):
    command = (
        f"vsearch --fastq_stats {fastq_file} "
        f"--log {output_log}"
    )

    subprocess.run(command, shell=True, check=True)
    print(f"Statistics generated and saved to {output_log}")

# Example usage
fastq_file = '../results/merged_filtered.fastq'
output_log = '../results/merged_filtered_output_stats.log'
fastq_stats(fastq_file, output_log)

[frederic-mahe]

Filtering before merging can lead to issues due to unequal read lengths, which can prevent successful merging.

This is not exactly what caused the second issue Fatal error: More forward reads than reverse reads. By filtering your forward and reverse datasets independently, you desync'ed them. Trimming probably resulted in empty fastq entries (zero-length reads) in your reverse dataset, and empty reads are discarded by vsearch. The trimmed reverse dataset now contains less entries than the trimmed forward dataset, which triggers an error when trying to merge the two.

When working with pairs of fastq files, vsearch can process them simultaneously, to preserve synchronization (note the --reverse):

vsearch \
    --fastx_filter <(printf "@fwd\nAA\n+\nI?\n") \
    --reverse <(printf "@rev\nTT\n+\nII\n") \
    --quiet \
    --fastq_truncqual 30 \
    --fastqout - \
    --fastqout_rev -

What key metrics would you recommend examining to verify the quality of the FASTQ data after filtering and merging?

There are several available quality metrics (see FastQC). I work mostly with metabarcoding data, and I tend to focus on the expected error (EE). It is a measure of the risk of errors in a given sequence. I use it to filter out sequences that are both rare and probably erroneous (high EE, or high (EE / length), as EE is length-dependent).

[simmimourya]

This is what I had come up with for my usecase:

# Quality filter paired end reads
# Using qscore_threshold = 30
def filter_reads(r1_path, r2_path, output_dir, params, filter_override=False):
    """
    Apply quality filtering to paired-end reads.

    Generates forward and reverse filtered reads in the specified output directory.
    """
    forward_filtered = os.path.join(output_dir, "forward_filtered_read.fastq")
    reverse_filtered = os.path.join(output_dir, "reverse_filtered_read.fastq")

    qscore_threshold = params['fastq_qscore_threshold']
    command = (
          f"vsearch --fastq_filter {r1_path} --reverse {r2_path} "
          f"--fastqout {forward_filtered} --fastqout_rev {reverse_filtered} "
          f"--fastq_truncqual {qscore_threshold}"
    )
    subprocess.run(command, shell=True, check=True)
    
    return forward_filtered, reverse_filtered

Follow-up: Looks like the --fastq_truncqual {qscore_threshold} is truncating the reads at the first base that falls below the specified quality threshold. This results in the read being shortened instead of entirely discarding it if it contains low-quality segments.
We instead want to discard the reads that are below a certain quality threshold (30), basically the q30 filtering should filter out reads, not bases. I am approaching this problem using --fastq_maxee. Basically discard the reads exceeding the specified maximum expected errors or max_expected_errors based on quality scores. But I am not quite sure what should be the value of max_expected_errors since vsearch technically will not use qscore_threshold for --fastq_maxee.
Am I approaching this correctly?

[simmimourya]

Please let me know if you need any specific error logs! :)

[colinbrislawn]

Hi Simmi. In case you are still working on this:

We instead want to discard the reads that are below a certain quality threshold (30), basically the q30 filtering should filter out reads, not bases.

Congratulations! You are well on your way to replicated the finding of Robert Edgar, who made the usearch tool vsearch is based on! Filtering full reads, not bases, is correct! Next, because you can't average q-scores, Expected Error filtering is recommended instead.

But I am not quite sure what should be the value of max_expected_errors since vsearch technically will not use qscore_threshold for --fastq_maxee.

Q30 * 250 bp = correct base pairs
0.999 * 250
249.75 = correct bases
expected error of 0.25

And this is exactly why average q-scores is a bad idea. So ignore that and use EE instead!

So how many errors are okay to keep in a sequence? Well, here is how many basepairs differentiate OTUs at two identity thresholds:

250 bp OTUs clustered at 99% = 250*0.01 = 2.5 bases
250 bp OTUs clustered at 97% = 250*0.03 = 7.5 bases

So using an Expected Error threshold under that of your clustering threshold makes sense to me.

The default value is also defensible.

[torognes]

I can add a fastq_minqual option to fastq_filter/fastx_filter that will discard reads with any base having a quality score below the given threshold.

[torognes]

I have now added the --fastq_minqual option in vsearch version 2.30.0. It should do exactly what you want.

[frederic-mahe]

--fastq_minqual tests added to our test-suite (frederic-mahe/vsearch-tests@9ecdd15)