Question: Sequence Number Count In Fastq.Gz File
gravatar for Bioscientist
6.3 years ago by
Bioscientist1.6k wrote:

Do we have any easy, fast way to know how many sequences contained in paired-end fastq.gz file? One simple way I think to calculate is to count the # of lines in fastq file divided by 4. (Four line for one sequence)

But I'm afraid it'll be quite slow for the calculation due to the large size of fastq.gz

What about the sequences in .sai/.sam file? (Because I want to guarantee the quality by comparing sequences in fastq file, and how many sequences been processed in .sai/.sam file)

The only way I can think of, is that we have the error report of .sai/.sam file generation, in which we get the information like "10000000 sequences have been processed". I can print this line by python, but how can I print out this number out of the line?


To Ale: thx, but what's the rationale for your first command? I mean about ^@ Just look at the fastq.gz file below, where is ^@ ? I see many @, but not ^


Edit again: Also I try both of the commands, and find the result is 2-fold difference. I compared with .sai file error report, which said 853813 sequences have been processed, then should say 'zgrep' command double counts.

[x@hpcc01 bwa2]$ zcat ERR002356_1.recal.fastq.gz | echo $((`wc -l`/4))
[x@hpcc01 bwa2]$ zgrep -c '^@' ERR002356_1.recal.fastq.gz

EDIT: THX. But I'm not dealing with one file. Plz look at below:


import os,sys,re
import gzip
import commands

path = "/home/xug/nearline"

for file in os.listdir(path):
  if re.match('.*\.recal.fastq.gz', file):
    fullpath = os.path.join(path, file)
    result = commands.getoutput('zcat fullpath |wc -l')
    numseqs = int(result)/4.0
  print numseqs

Problem is, I define 'fullpath' here for all fastq files, but after being put under ' ', seems this fullpath doesn't work...How can I solve this problem? thx

Problem solves....fullpath is variable. 'zcat '+fullpath+' |wc -l'

fastq sequence counts • 31k views
ADD COMMENTlink modified 2.0 years ago by Biostar ♦♦ 20 • written 6.3 years ago by Bioscientist1.6k

@bioscientist, this is a very good point - [z]grep cannot be used (like for fasta) because @ and + can also occur in the quality line.

ADD REPLYlink written 6.3 years ago by Aleksandr Levchuk3.0k

^ means --> begins with ...

ADD REPLYlink written 6.3 years ago by Aleksandr Levchuk3.0k

The following should do the trick:

zgrep -c '^@.*/[0-9]$' fastq.gz

ADD REPLYlink written 6.3 years ago by Andreas2.3k

In general, checking for a terminal digit will not work because Sanger-encoded Fastq qualities include characters in the range ASCII 33 to 126, that is to say the qualities may contain digits. Your regex will still potentially match quality lines.

ADD REPLYlink written 6.2 years ago by iw9oel_ad6.0k
gravatar for Aleksandr Levchuk
6.3 years ago by
United States
Aleksandr Levchuk3.0k wrote:

Your method should also work and it should be pretty fast:

zcat my.fastq.gz | echo $((`wc -l`/4))

For SAM it's:

cat my.sam | grep -v '^ *@' | wc -l

I'm not sure how to interact with .sai files.

NOTE1: You can call the above Bash commands from Python:

import commands
result = commands.getoutput('zcat my.fastq.gz | wc -l')
numseqs = int(result) / 4.0

If you are doing Bioinfromatics you want to have the GNU toolset which includes the Bash shell (it's default on Linux and Mac, on Win you can install Cygwin)

NOTE2: I'm assuming that all input files have a newline at the end. It's easy to check: run

zcat my.fastq.gz | tail -n 1 | cat -E

if you see a dollar sign ($) then everything is OK. Conversely, you may have generated the files on Windows or from a newbie/ignorant code that does not place a newline at the end of the file, then you will not see a the dollar sign. In that case sanitize the files like this:

zcat my.fastq.gz  | ruby -pne '$_.chomp!; $_ << "\n"' | gzip > good.fastq.gz

NOTE3: I'm assuming that your FASTQ files don't wrap lines. One line per read. You can test like this:

zcat my.fastq.gz  | ruby -e '
skip_line = false
last_char = "+"

while gets
  if skip_line
    skip_line = false

  if $_ =~ /^@/ and last_char == "+"
    last_char = "@"

  elsif $_ =~ /^\+/ and last_char == "@"
    last_char = "+"

    STDERR.puts "WARNING: fastq lines wrap at line #{$.}" 

  skip_line = true

STDERR.puts "fastq lines do not wrap"

NOTE4: Instead of all of the above you can use a filesystem with online compression (i.e. ZFS) and do all the analysis on fastq files instead of fastq.gz

ADD COMMENTlink modified 6.2 years ago • written 6.3 years ago by Aleksandr Levchuk3.0k

Uh, you are making the assumption that the Fastq file isn't wrapping the sequences. Although this is uncommon, it's allowed within the standard. Also, SAM files may have a separate header section, you need to make sure only to count the aligned reads.

ADD REPLYlink written 6.3 years ago by Ketil3.8k

thx Ale, but this is shell command, right? I want to do this within Python script.

ADD REPLYlink written 6.3 years ago by Bioscientist1.6k

yeah the line counting thing to easily flummoxed by missing end-of-file carriage returns (bad but common) and the SAM thing is just outright wrong in headered SAM files.

ADD REPLYlink written 6.3 years ago by Jeremy Leipzig17k

@Jeremy, @Ketil - I took into account your comments.

ADD REPLYlink written 6.2 years ago by Aleksandr Levchuk3.0k
gravatar for summerela
3.1 years ago by
United States
summerela70 wrote:

I'm a bit late to the party, but here's an answer using gnu parallel that should do the trick.. and quickly: 

parallel "echo {} && gunzip -c {} | wc -l | awk '{d=\$1; print d/4;}'" ::: *.gz

This assumes your fastq files are gzipped and in standard fastq format. 

ADD COMMENTlink modified 3.1 years ago • written 3.1 years ago by summerela70

tks summerela \o/

ADD REPLYlink written 20 months ago by renatopuga10
gravatar for Matt Shirley
3.1 years ago by
Matt Shirley7.8k
Cambridge, MA
Matt Shirley7.8k wrote:

Here's a method to count the number of reads from the BAM index file:

The speed of this method is independent of the number of reads in your BAM file.

ADD COMMENTlink written 3.1 years ago by Matt Shirley7.8k
gravatar for Alex Reynolds
3.1 years ago by
Alex Reynolds20k
Seattle, WA USA
Alex Reynolds20k wrote:

I posted a solution here that uses GNU Parallel: C: Gnu Parallel - Parallelize Serial Command Line Programs Without Changing Them

ADD COMMENTlink modified 3.1 years ago • written 3.1 years ago by Alex Reynolds20k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 942 users visited in the last hour