Delete all 4 lines of a fastq read from a fastq file using read ID

Question

I have the following error when running bowtie2:

Error: Read HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG has more read characters than quality values.

I now want to remove all 4 lines of this specific read from the fastq file.

How can I use awk or sed to do this?

Answer 1

Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ. For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.

If the FASTQ is malformed, then you should really ask yourself how this happened in the first place and fix the source of the problem. If the record is valid FASTQ, then I suggest parsing the read with for example FastqGeneralIterator and dumping the parsed result back to FASTQ in a 4-line-per-record form.

Answer 2

If you are 100% sure the read only has 4 lines (they can have more), you can use this sed command:

sed -i.bak '/^@HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG/,+3d'

The -i.bak makes sed modify the original file and create a backup copy with the same name and the extension .bak. The command just means "delete the line matching the pattern and the next three lines".

Answer 3

When I can I like to do file processing line-by-line in the spirit of UNIX tools. You can read 4 lines from a Fastq file into 4 tab-separated values on a single line using paste, and then use grep to filter out the record in question. (Then you just have to turn the tabs back into line breaks.)

paste - - - - < reads.fastq \
    | grep -v 'HWI-D00466:116:CC62WANXX:3:1102:7363:63646' \
    | tr '\t' '\n' \
    > reads-fixed.fastq

However, it's likely that other reads in your Fastq file are corrupted as well, in which case it's probably better to write a script in Python or some other language that discards all reads whose read length doesn't match the quality string length.

But of course the most important question is: how did the data become corrupted in the first place? Where did this file come from, and what pre-processing steps were performed on it? Is getting rid of this read (or all reads with mismatched sequence/quality length) going to solve the problem for you? These aren't questions that we can answer for you, but probably deserve good answers. :-)

Answer 4

Here's an awk one-liner that works as intended, unlike the bioawk answer.

This removes all instances of fastq records in which the length of the seq and qual fields do not match. Again, append | gzip > filtered.fastq.gz to gzip the output.

zcat bad_file.fq.gz | awk '{ pos = NR%4; if ( pos == 1) {h=$0} else if (pos == 2) \
  {s=$0} else if (pos==3){c=$0} else if ( pos==0 ) \
  {q=$0; if ( length(q) == length(s) ){printf("%s\n%s\n%s\n%s\n", h, s, c, q)}}}'

Delete the \ characters at the end of the first two lines to make it easy to paste into your terminal.

Here is a shorter version of the same thing thanks to @terdon

zcat file.fastq.gz | awk '{ record[++k]=$0; if(NR%4==2){slen=length($0)} else if(NR%4==0){if(length($0) == slen){ for(i in record){print record[i]}} k=0; }}'

Answer 5

To generalize this to remove all reads that do not have the same length of read and sequence quality, you can use this bioawk one-liner:

bioawk -cfastx '{ if (length($seq) == length($qual) ){printf("@%s %s\n%s\n+\n%s\n", $name, $comment, $seq, $qual)} }' my_fastq.gz

This strips anything in the "fastq comment field", +, but that almost never has anything there anyway.

Append | gzip > filtered.fastq.gz to the above command to gzip the output.