How to Read FASTQ file using Python for Quality controlling

Q score and Nucleotide graph

Hi, I’m Aniket and in this story, we will try to read FASTQ file using only python method. Basically, in bioinformatics, FASTQ file is the text-based file format which is used to store sequence and their corresponding quality score (each encoded as an ASCII codes). For example:

1. @_SEQ_READNO_TIME| THIS IS THE LABEL length=52
2. ACTAGCGAGCGTACGACTGCTATGCTCTAGTAGCCCGGTATTTTGCATGCTCA
4. AAAAAAAA:9&@11!!;;<:.>>>AA’’””’’;??//.??./#$$556&&*89(**6A!2??<281AAS79(

This is the single read of a FASTQ file. Where line 1 is called header or label which is denoted as a symbol @ and just after line 2 will be sequence of a single read, in line 3 plus (+) symbol is use as a separator to separates the sequence and quality scores. Line 4 shows quality score or Phred score or Q score, this Q score was generated by the sequencing machine in the form of ASCII character through the probabilistic outcome of each base. In a real case of *.fastq file can have more than one reads, just like up to 1,000 or more.

With the help of python programming, we’ll analyze these FASTQ file format to check the quality of the sequence and find out the length of that sequence, dividing a FASTQ file into multiple sub files, trimming and so on.

let’s create a first function that’s return a dictionary of read of sequence: and corresponding ASCII code:

def readFASTQfile(fastq_file):
        # seqData is the empty list use for storing data of sequence and ASCII only,
        # seqData variable drop all the sep (+) and header line starting from '@'
        seqData = []
        # seq and ascii list variable stores sequence and thier ASCII code seperately...
        (seq, ascii) = ([], [])
        with open(fastq_file, 'r') as FASTQ:
            # open FASTQ file and read as 'FASTQ', reading start from 2nd lines and store in 'line' variable
            # ie. try to drop first header.
            line = FASTQ.readlines()[1:]
            for i in range(len(line)):
                if i%2 == 0:         # seqData stores those line which contain ASCII and sequence both (even lines)
                    seqData.append(line[i].replace('\n', ''))    # replace all new lines (\n) of the each lines

        # this is the seperate part of this program, makes a 'seq' and 'ascii' list from 'seqData', and than 
        # a dict (seqdict) from 'seq' and 'ascii'...
        for i in range(len(seqData)):
            if i%2 == 0:       # this part is same as above even line part, where 'seq' store sequence data from 'seqData'
                seq.append(seqData[i])
            elif i%2 != 0:      # except all (ASCII CODE) stores 'ascii'
                ascii.append(seqData[i])
        seqdict = {key: value for key, value in zip(seq, ascii)}  # use dict. comprehension to make dict of 'seq' and 'ascii'
        # return dict of final sequence and thier respective ASCII code for finding the phred score of that sequence
        return seqdict

f = 'out_1000_ERR101899.1.fastq'
print(readFASTQfile(f))

output:
{'AACGTAATGATTGGCAAGACTAATTTCATTGGATGTTCTGGCGCAAATGTATGCACAAAAGCATCACCATAAGGTGGATTAGAAACAAACATCACATATGCCACAACTAAAATCATCACAATACCAAGAATCATTGAAACAACGTCCATA': '@CCFFFFFGGHGHJIGIJIIIJIJIJIGIIJJIJJJJHGGIFIIHIJIIGGGG@GGIIJIGIJGIC@EHEFHHFCDFFEDEEDEDDDDDB?CC@@CDCADCDCCBD@CCDC@>CCCCCDDDD>C:>:<B?A@CDCCCCCCD?858<<BDC', 'CTTTATTACAGAGTGAATCGGATTGGTGAAAATCGAAATTTTGAGATTTTTACCAATTCGATTTTTTTCATAGAAATTAAAAAAGCCAACAAGGCTCTTGAAACCTTGTTGGCGTAAAACTTAGTCATCACTAATTAGTGAATGAAGTTA': '@@CFFFFFHHFFHGIHIIJIJEDHHIHHGIIIGHIDHIGIJJFBFBGHIJJGHIIGCGDFCFFGIJHFBCDD@>CECEDDDDDD=B<CD<B@?BB:8CC@C@ACBDDDCD?B4959BDD>ACCCCCDEDCD>@CDD@C@C@CDCCC@###', 'ATCCACCATCGCCAGTAATCGGCAAACGTCACATGTGAATTCGGCATTATAAATAATGATGCAAGATAAAATAACGGAATTGCAATTGCTCCAACAAACAATAGCGTCAATAAGTGACGTTTATCATGATTCACTTGTACTTTGTTGCGA': ';??ADBDDDDFFFFGFFGG??@EG@D?@FFFBFIIBFGGGCDFHAEHIGIBHIIHIIGIEE?ACED@DB;C>CC;A@@BBCC:ACCDC:@>>@>C?<<88ACC@>0985<>>@>(+:4<<?CBDCC@>>A>:@@CC>:44::3::?<:<5'}

Logic behind this function, I’m created two list of seq and ascii for storing sequences of each reading and their ascii chars equally in seq and ascii lists, than using even and odd method for just after first line I mean except first line of this file, even line stored in seq list rather odd line is appended in ascii list.. (don’t forget to replace or remove new line (‘\n’) character from each line in the FASTQ file).

Now, How we can find the length of the FASTQ file:

Here, we can add some other features like, ‘The length of FASTQ file, also known as numbers of reads in a file’; Most of the important thing is that, python dictionary stored unique key, value pairs, So here we don’t need to remove duplicates. just simply use len() function of python for above output dictionary.

fastq_ = 'ERR101900_1.fastq'
fastq_dict = len(readFASTQfile(fastq_))
print(fastq_dict)

# output will be: 319739

In this ‘ERR101900_1.fastq’ file have 319739 numbers reads.

A another feature, to select some of the above reads by head() method:

Also, I’m created a head(n) function for selecting n number of the above reads.

def head(fastq_dict, top=5):
        # takes an argument 'top' for reading or returning limited sequence from top 
        # read and select all fastq and makes a list of total tuple items
        lis_limit = list(fastq_dict.items()) 
        return dict(lis_limit[:top])  # and create a dictionary of top limited list

fastq_dict = readFASTQfile(fastq_)
print(head(fastq_dict))
print(len(head(fastq_dict)))

output:
{'AACGTAATGATTGGCAAGACTAATTTCATTGGATGTTCTGGCGCAAATGTATGCACAAAAGCATCACCATAAGGTGGATTAGAAACAAACATCACATATGCCACAACTAAAATCATCACAATACCAAGAATCATTGAAACAACGTCCATA': 'CCCFFFFFGHFHHIHIGGGIJGIJJJJJJJIJJJJGIIJHCGGIGIJJIGIJIJIJJIJJIJJJJJIHHHHHHF@DFFEDEEEEDDDDDDDDDDDDDCEDDDDDDBDCCCDDEDDDDCDDDCCCDDDBDDDCDCACADCDDD@@8ABBCC', 'CTTTATTACAGAGTGAATCGGATTGGTGAAAATCGAAATTTTGAGATTTTTACCAATTCGATTTTTTTCATAGAAATTAAAAAAGCCAACAAGGCTCTTGAAACCTTGTTGGCGTAAAACTTAGTCATCACTAATTAGTGAATGAAGTTA': '@@CFFFFFHHFFHGIHIIJIJEDHHIHHGIIIGHIDHIGIJJFBFBGHIJJGHIIGCGDFCFFGIJHFBCDD@>CECEDDDDDD=B<CD<B@?BB:8CC@C@ACBDDDCD?B4959BDD>ACCCCCDEDCD>@CDD@C@C@CDCCC@###', 'ATCCACCATCGCCAGTAATCGGCAAACGTCACATGTGAATTCGGCATTATAAATAATGATGCAAGATAAAATAACGGAATTGCAATTGCTCCAACAAACAATAGCGTCAATAAGTGACGTTTATCATGATTCACTTGTACTTTGTTGCGA': ';??ADBDDDDFFFFGFFGG??@EG@D?@FFFBFIIBFGGGCDFHAEHIGIBHIIHIIGIEE?ACED@DB;C>CC;A@@BBCC:ACCDC:@>>@>C?<<88ACC@>0985<>>@>(+:4<<?CBDCC@>>A>:@@CC>:44::3::?<:<5', 'ATATAGAACGTAATCATATTATGATATGATAATAGAGCTGTGTAAAAAAATGAAAATAGACAGTGGTTCTAAGGTGAATCATGTTTTAAATAAGAAAGGAATGACTGTACGATGAGCTTTGCAGCAGAAATGAAAAATGAATTAACTAGA': '<@@DDDDDHHFDHIIGIHI>FHIGGIIEEIGCBCFEBFIDDHGHIIIIIIIEIIGGIEHGH;@=@@7DCEHEEE.;@DEDCC@.>C;A@CDDCC>;?C?(:ACCCCCCCCBB8?>>ACAACCD((889?>CC4>:>ABACCD(:CCCAA3', 'GTAATAAAATAATACGAATTACCAATACAAGGATAATAATAGCTAAACCATAATTGTCGTTTAATAAGTTATTTCCCAACCAATCCAATACATTTTTCATTGGATCTACGAATGTATTGTAGAAAAACCCAGTACGTTTTTCATGTTTAG': '@@@DDDDD?,AFBH9?FGH=B@DHI@FBHH<EFHIHCHHIE4?FGIIIGDAGI>G:DFAFFHI:@GIDGCEEGG>EHFFF<BD?C;@>CCACD;(;@;@>@::ACCCCC(528(4:::>:>:C3:?B9?00(+::(+<CBB#########'}
5

All above using object-oriented:

# created a class for reading FASTQ file format
# with takes an argument of file location from the local machine
class FASTQformat:

    def __init__(self, file_Location):
        # initialize argument of file location
        self.fileLocation = file_Location

    # function of reading FASTQ from FASTQformat class, returns dictionary of sequence and thier respective ASCII characters
    # ie. dict = {'seq': 'ASCII', ....}
    def readFASTQfile(self):
        # seqData is the empty list use for storing data of sequence and ASCII only,
        # seqData variable drop all the sep (+) and header line starting from '@'
        seqData = []
        # seq and ascii list variable stores sequence and thier ASCII code seperately...
        (seq, ascii) = ([], [])
        with open(self.fileLocation, 'r') as FASTQ:
            # open FASTQ file and read as 'FASTQ', reading start from 2nd lines and store in 'line' variable
            # ie. try to drop first header.
            line = FASTQ.readlines()[1:]
            for i in range(len(line)):
                if i%2 == 0:         # seqData stores those line which contain ASCII and sequence both (even lines)
                    seqData.append(line[i].replace('\n', ''))    # replace all new lines (\n) of the each lines

        # this is the seperate part of this program, makes a 'seq' and 'ascii' list from 'seqData', and than 
        # a dict (seqdict) from 'seq' and 'ascii'...
        for i in range(len(seqData)):
            if i%2 == 0:       # this part is same as above even line part, where 'seq' store sequence data from 'seqData'
                seq.append(seqData[i])
            elif i%2 != 0:      # except all (ASCII CODE) stores 'ascii'
                ascii.append(seqData[i])
        seqdict = {key: value for key, value in zip(seq, ascii)}  # use dict. comprehension to make dict of 'seq' and 'ascii'
        # return dict of final sequence and thier respective ASCII code for finding the phred score of that sequence
        return seqdict

    # function to return length of the FASTQ file 
    # read sequences in single file
    def FASTQlen(self):
        fastq = self.readFASTQfile()   
        # returns the length of the dict...
        return len(fastq)

    # function for returning top most limit of sequence where default value of top reads is 5
    def head(self, top=5):
        # takes an argument 'top' for reading or returning limited sequence from top 
        fastq_limited = self.readFASTQfile()
        # read and select all fastq and makes a list of total tuple items
        lis_limit = list(fastq_limited.items()) 
        return dict(lis_limit[:top])  # and create a dictionary of top limited list

fastq_ = 'ERR101900_1.fastq'
FASTQ = FASTQformat(fastq_)

# print(FASTQ.readFASTQfile())
# FASTQlen() returns 0 if file will empty..
print(FASTQ.FASTQlen())
# head(top=3) -> return three reads from top [return {} if file will empty]
print(FASTQ.head(top=3))

You can check my next story here to visualize FASTQ quality scoring graph by matplotlib and checking quality by Phred score (Q score) of based 33 and 64.