This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Number of bases with a certain quality in FASTQ file

Hi, I'm wondering is there a convenient way to count quantity of bases in FASTQ file with Phred quality score equal or higher then 30? Thanks!

next-gen sequencing
• 4,509 views • link
updated by h.mon h.mon • written by Denis Denis

FastQC provides a graphical quality per base position summary, and if you dig into its output files, you will find mean, median, lower quartile, upper quartile, 10th percentile and 90th percentile quality per base position.

Why do you need total count?

• 1 views • link
written by h.mon h.mon

We often use the number of bases with quality score >=30 as a quick assay of data quality. Besides, our sequencing facilities commonly use that statistic.

• 1 views • link
written by Denis Denis

2 answers 

 gunzip -c  input.fastq.gz | paste - - - - | cut -f 4 | fold -w 1 | awk '($1>="?")' | wc -l
  • gunzip -c input.fastq.gz decompress input
  • paste - - - - linearize 4 lines
  • cut -f 4 get the quality column
  • fold -w 1 fold to one column
  • awk '($1>="?")' select strings greater than "?" (ASCII code=63 = 33(base fastq)+30.
  • wc -l count
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks! But from where awk may know that for example "=" is less than "?" ? How can it compare just two characters?

• 1 views • link
written by Denis Denis

Based on ASCII code decimal values. It is comparing value contained in $1 (which would be individual Q score) for being more than or equal to to ?.

• 1 views • link
written by GenoMax GenoMax

Hi genomax! Thanks a lot! Now it's clear.

• 1 views • link
written by Denis Denis

Should i add something like this tr -d '\n' somewhere in the pipe to remove new line characters? echo \n | awk '($0>="?")' results in n

• 1 views • link
written by Denis Denis

reformat.sh from BBTools / BBMap package can also do this:

reformat.sh in=file.fq.gz qchist=file.qchist.txt

For paired end files

reformat.sh in1=file1.fq.gz in2=file2.fq.gz qchist=file.qchist.txt
• 1 views • link
written by h.mon h.mon

Log in to answer this question.

More posts like this