6
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ Run tail on your FASTQ file. If a job aborted while writing to disk, this read will be the last line. This error message likely means the file is truncated and missing reads. That read is not as much of a concern as the rest of the data. It’s best to run your data processing again to make sure the file is complete. $\endgroup$ – Tom Kelly Nov 23 '18 at 7:34
9
$\begingroup$

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.

$\endgroup$
6
$\begingroup$

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".

$\endgroup$
  • $\begingroup$ Maybe the "@" should be added to the pattern. What would happen if the read name was also on the "+" line? $\endgroup$ – bli Oct 31 '18 at 12:44
3
$\begingroup$

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. :-)

$\endgroup$
1
$\begingroup$

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; }}'
$\endgroup$
  • 1
    $\begingroup$ Heh, I was just writing one too :) Here's what I came up with, it's a bit shorter than yours: 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; }}' feel free to add it to your answer if you want to. $\endgroup$ – terdon Oct 31 '18 at 15:51
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This doesn't seem to work for me (Arch Linux, bioawk version 20110810). Even after correcting the syntax error (corrected in your answer), it seems to exit after the first record it finds where the lengths are not the same. Can you confirm it works on your machine? $\endgroup$ – terdon Oct 31 '18 at 14:27
  • $\begingroup$ Thanks for the edit! Yep, I just confirmed that it works for me on all records. I'm using Ubuntu 16.04. >>bioawk --version returns awk version 20110810. It looks like the last time that I compiled bioawk was from the github repo in 2016. $\endgroup$ – conchoecia Oct 31 '18 at 15:15
  • $\begingroup$ That's really weird. I have the same version. Can you try a case where the record with different lengths is the second one in the file? Does it not just print the first record there? $\endgroup$ – terdon Oct 31 '18 at 15:18
  • $\begingroup$ I got the same error as you - I made the second entry in the fastq file have a different length of seq and qual fields and bioawk just exited before printing that record. I wonder if this is intended behavior? I have to add that in my above comment, I just tried the bioawk one-liner on a normal fastq file without any length mismatches. Oops! $\endgroup$ – conchoecia Oct 31 '18 at 15:23
  • $\begingroup$ I don't see how it can be intended. Sounds like a bug in bioawk. $\endgroup$ – terdon Oct 31 '18 at 15:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.