Python Programming for Biology: Bioinformatics and Beyond

Reading sequences from FASTQ files

A FASTQ-format sequence file is read by using the imported reader class FastqReader  to

make an object that represents the open file. It is then a simple matter to extract each of

the sequence records that is represented in the file by looping through this object. Each of

---

the  seqRes  objects,  as  it  appears  in  the  loop,  is  a  SequenceWithQualities  class  from  the

HTSeq library, and as such it comes with lots of inbuilt functionality and some of this is

demonstrated below: printing the name, the sequence itself and the reverse complement of

the  read  (here  the  slice  notation[::-1]  gives  the  sequence  in  reverse  relative  to  the  main

sequence).

fileObj = FastqReader(fastqFile)

for seqRead in fileObj:

print(seqRead.name)

print(seqRead.seq)

print(seqRead.get_reverse_complement()[::-1])

This gives a result like:

r999

AGGATAATGAGGCGAGCCGGGGGAACTGAAANTGG

TCCTATTACTCCGCTCGGCCCCCTTGACTTTNACC

Given  that  these  sequence  records  come  from  reading  a  file  format  that  incorporates

quality  scores  we  can  naturally  interrogate  those  scores.  In  this  example  we  generate  a

graph of the mean score along the alignment positions. The meanQual initially starts out

as None and for the first sequence is set to a NumPy array of the scores seqRead.qual. For

subsequent  records  the  scores  are  then  added  (element  by  element  as  is  the  standard

NumPy way) to this array, so that at the end the whole array can be divided by numReads

to give the average value along the sequence, which is plotted with the pyplot library.

numReads = 0.0

meanQual = None

for seqRead in fileObj:

print(seqRead.qual)

if meanQual is None:

meanQual = array(seqRead.qual)

else:

meanQual += seqRead.qual

numReads += 1.0

if numReads:

pyplot.plot(meanQual/numReads)

pyplot.show()