Biostar Beta. Not for public use.
Question: Number of bases with a certain quality in FASTQ file
0
Entering edit mode

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!

ADD COMMENTlink 19 months ago Denis • 70 • updated 19 months ago h.mon 25k
Entering edit mode
0

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?

ADD REPLYlink 19 months ago
h.mon
25k
Entering edit mode
0

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.

ADD REPLYlink 19 months ago
Denis
• 70
5
Entering edit mode
 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
ADD COMMENTlink 19 months ago Pierre Lindenbaum 120k
Entering edit mode
0

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

ADD REPLYlink 19 months ago
Denis
• 70
Entering edit mode
1

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

ADD REPLYlink 19 months ago
genomax
68k
Entering edit mode
0

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

ADD REPLYlink 19 months ago
Denis
• 70
Entering edit mode
0

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

ADD REPLYlink 19 months ago
Denis
• 70
1
Entering edit mode

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
ADD COMMENTlink 19 months ago h.mon 25k

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.0