2
$\begingroup$

I want to remove reads from FASTQ file that contain homopolymers > 10bp and remove reads with <35 average quality score across the entire read.

To remove homopolymers > 10bp, I tried this on a Linux machine, but it only removes the sequence line:

zcat file.fastq.gz | 
  awk '!/A{10,}/&&!/C{10,}/&&!/G{10,}/&&!/T{10,}/ {print}' cleaned_file.fastq
Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Please edit your question and explain what you need in more detail. What do you mean by "less than the 35 quality score"? Do you want to remove all reads that have even a single base with <35 score? Do you want to calculate the average quality across the entire read and remove reads with <35 average quality score? Something else? Why do you want to remove these from the fastq file instead of handling it in the aligned bam or as part of your variant calling process? How are you defining homopolymers here? Do you just want to remove reads if they have >10 consecutive identical bases? $\endgroup$
    – terdon ♦
    Jan 25 at 0:12

1 Answer 1

Reset to default
4
$\begingroup$

Use an off the shelf tool for read preprocessing.

Here is one:

./fastq_qual_trimmer -i test.fq -m 35 -H 10

That one does exactly what you want but seems a little old/unmaintained, so here is another that does more or less the same thing:

fastq-mcf --qual-mean 35 --homopolymer-pct {X} adapters.fa reads.fq

where {X} is 10 / read length, adapters.fa is an adapter file (which I believe can be empty or filled with dummy sequences).

You could also use a library like biopython or dnaio to write a quick script to do this, but it hardly seems worth it.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for your help. Sorry I didn't write in detail but i need to calculate the average quality over the whole read and remove the reads with average quality score <35. As I understand it, these tools can't handle it. $\endgroup$
    – user16821
    Jan 25 at 19:45
  • $\begingroup$ @n_n sorry, used wrong options. Both of these tools can handle it across the read. Updated answer. $\endgroup$
    – Maximilian Press
    Jan 25 at 23:04
  • 1
    $\begingroup$ Thank you so much for taking the time to help me. Your assistance is greatly appreciated. $\endgroup$
    – user16821
    Jan 26 at 12:26
  • $\begingroup$ @user_01 if this answer has solved your problem could you very kindly please mark this as accepted? It would be a small way of showing gratitude and helps the site in several ways. $\endgroup$
    – M__ ♦
    Jan 27 at 0:53
  • $\begingroup$ Just one more question I tried fastq-mcf as you suggest for my fastq files but i got this Error: number of input files must match number of '-o' output files. May you help me to understand whats wrong with code. for file in *.fastq.gz; do read_length=$(head -n 2 ${file} | tail -n 1 | wc -m) homopolymer_pct=$((8/${read_length})) fastq-mcf ${file} adapters.fa --qual-mean 35 --homopolymer-pct ${homopolymer_pct} -o filtered/${file} done $\endgroup$
    – user16821
    Jan 27 at 21:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.