This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Pull Out Qual From Multi Line Fastq?

From http://biostar.stackexchange.com/questions/1389/how-to-generate-a-consensus-fasta-sequence-from-sam-tools-pileup We use awk and perl to get fasta by

awk '/^@chr1$/,/^+$/' consensus.fastq | perl -pe "s/@/>/ ; s/\\+//" > chr1.fasta

The sample data

@chr1
nnnnnnnagatagaaataCACGATGCGAGCAATCAAATTTCATAACATCACCATGAGTTT
GGTCCGAAGCATGAGTGTTTACAATGTTTGAAtaCCTTATACAGTTCTTATACATACTTT
ATAAATTATTTCCCaagctgttttgatacactcactaacagaTATTCTATAGAAGGAAAA  
+
!!!!!!!????????BBBEEEHKKKKKKKKKKKKNNNNNNQQNNNNNNQWWWW7WZWWWZ
ZZZ]]]`````````>]]]]]]]ZTQQQTQNNKKKKKKHKKHHHEEEEEEEHHHHHHHHH
HKKHHHHHHEEEEEBBBBBBBBBBBB?=B@BBBBBB??BBBBEEEEEEEEEEHKNNNNNQ

Now I want to extract the qual for each chromosome. But the qual score for each chromosome all beginning with "+", how can I extract them?

Thanks.

fastq quality
• 6,886 views • link
updated by Haluk Haluk • written by Zhshqzyc • Zhshqzyc

It quality doesn't wrap lines can't you just grep for "+"?

• 0 views • link
written by Zev.Kronenberg Zev.Kronenberg

2 answers

Here is quick and extremely dirty python script to pull out fastq quality not making assumptions about line wrapping:

import sys

inFile = open(sys.argv[1],'r')

header = ''
qual = ''

seqs = False
quals = False
for line in inFile:
    if line[0] == "@":
        if header != '':
            print ">" + header
            print qual
        header = line[1:].strip()
        quals = False
        qual = ''
    elif line[0] == "+":
        quals = True
    else:
        if quals:
            qual += line.strip()

print ">" + header
print qual

Save as yourName.py and run it by: python yourName.py yourFile.fastq > output.quals

edit * There was a mistake in the code where the qual variable wasn't being reset.

edit **

Here is the script to get just one entry:

import sys

inFile = open(sys.argv[1],'r')
queryHeader = sys.argv[2]

header = ''
qual = ''

seqs = False
quals = False
for line in inFile:
    if line[0] == "@":
        if header != '':
            if header == queryHeader:
                print ">" + header
                print qual
                sys.exit()
        header = line[1:].strip()
        quals = False
        qual = ''
    elif line[0] == "+":
        quals = True
    else:
        if quals:
            qual += line.strip()

print ">" + header
print qual

Save it as yourName.py. Use it by: python yourName.py yourFile.fastq "entry name"

• 0 views • link
written by Damian Kao Damian Kao

My fastq file is 6GB, I got the quals file is more than 12GB. Is it right?

• 0 views • link
written by Zhshqzyc • Zhshqzyc

No, it is 37.8GB. Could you please double check your code?

• 0 views • link
written by Zhshqzyc • Zhshqzyc

Oops, there was a mistake. I've fixed it.

• 0 views • link
written by Damian Kao Damian Kao

Now I got:There are many returns, is it correct?

chr2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ]]]`G]]`````cfTcffcfcOilifffEciiiififfcccccccfccffil Ullllliorrrrooiiooorueuuuuruxuuuxxrrrroorrrrrolllliilloooororoololoroollllllllllliilllrr QQQQQQQHQQQQQQQQQQQQQQQQTTTTQQQQQHQNNNNNNNNNNNNNKKKKKKKKHHK

NNKKKNTTTTWWWQWZZZZT`]ccccifcllllilllooorrooooorlflllooooil

]]]WWWWWTTTTTTTQWZZ]Z]]```]]Z]Z]]]]]]]]]cccfcliioouuuuuuu

{xxuuuuuruxxxx{{xx{xuuurrorlllloooorrrrooooooooouorrloouuuu

• 0 views • link
written by Love • Love

updated: If I only want to get qual for chromosome 1, how to change code?

• 0 views • link
written by Love • Love

I've added the code to just get one entry.

• 0 views • link
written by Damian Kao Damian Kao

I've added the scrip for getting just one entry. Make sure you typed the name of the entry correctly with quotes if there are spacing. If you didn't type it correctly, it will just give you the last entry.

• 0 views • link
written by Damian Kao Damian Kao

Many many many thanks.

• 0 views • link
written by Love • Love

line[0] can be '+' for one of the qual lines if the phred33 score is 10. So maybe this code may break in some cases.

0 1 2 3 4 01234567890123456789012345678901234567890 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI

• 0 views • link
written by Chris Penkett • Chris Penkett

line[0] can be '+' for one of the qual lines if the phred33 score is 10. So maybe this code may break in some cases.

• 0 views • link
written by Chris Penkett • Chris Penkett

Also @ is ASCII 64 which is 64 - 33 = PHRED Q 31. The quality line can contain and begin with both + and @.

• 1 views • link
written by Shaun Jackman Shaun Jackman

fq2fa tool extracts qual scores from fastq file.

Here is the link : http://www.phyloware.com/Phyloware/XSTK.html

• 1 views • link
written by Haluk Haluk

Log in to answer this question.

More posts like this