Program to process fastq files

Question

My program evaluates fastq files. The fastq file is just a huge text file (about 1,5gb in my case), which looks like this:

[email protected] 1 length=490
NTCCGGATGATGTCGCTGTTGCTGACAATGGTAATACGTTGACGGGGCAATATGCAGTTCGCTGCATACCGGTCCGACCCCGTACTGCTCACGCAGCTTATCCAGCAGTGGCATCATTTTTTCCAGAGGCGGTCGAACTCCGCCTTCGCAAAAAAAAAGGGAGCCCGGCGGAGGAGAACGTTACTGCGGCGGAGGTTACGATTTTTCCGGTTCCGCTCCTTTAGAAGCCGGACGTCTACCCGGCTCTTTTTGCTGAACGTCAGCGTCTGAAAGAGCTGGAACGTGAAAATCGTGAACTGCGCCGCAGTAACGATATCCTTCGCCAGGCTTCCGCTTATTTTGCGAAGGCGGAGTTCGACCGCCTCTGGAAAAAATGATGCCACTGCTGGATAAGCTGCGTGAGCAGTACGGGGTCGGACCGGTATGCAGCGAACTGCATATTGCCCCGTCAACGTATTACCATTGTCAGCAACAGCGACATCATCCGGAT
+SRR587217.1 1 length=490
BT[^^^^^caccccacffffffhhhgghhhhefeggghhghghgfgggghfgggeghhhhhgggghhgggggggggggggggdgggggfffffffffffffddfffffbfffffffffffffffffffffffMFFFHOIIPHMWHOIIHHMM^IFTMFHOHIHMHWMWFMFFFTHHIIIHIOIHIIIHMFFFFFFOIIIHHIOH^IIIHMSFOIHMHMFPYIIIIIIIIFFMFMOWIIIPIHFFMOIYIII^^``^aaaccccccccffffffghghhggghghgghghhhggggefgfededbe`cgeghggegaefgghegggggggggccdggggceeggfT``dddfffdfffffac]Zbffdfdbdbdfffddd^^dfdfdfddffffdbdbdfffffffbdddf``c`b[]ccQW`MSHYbdffff]cbbdbbbdYI^bd^[WdbbbfMXbbbf^b   ^^b^bdbdddIX^d]W]QW^dd^^dcfW
@SRR587217.4 4 length=502
NTCCATGACTTCTGCACGCGTCGGCATCGGGTTAGTAATCATTGACTCCATCATCTGGGTCGCCGTGATTACCGCTCGGTTTAGCTGACGCGCACGACGGATCAACGCTTTCTGAATGCCGACCAGTTCCGGGTCGCCAATTTCCACACCGAGGTCGCCACGTGCAACCATTACCACGTCAGAGGCGAGGATGATGTCATCCATTGCATCCTGGCTGCAAACGGCTTCCGCACGTTCAACCTTGGCAACAAGTGATGCGAAAATTGTTGCCAAGGTTGAACGTGCGGAAGCCGTTTGCAGCCAGGATGCAATGGATGACATCATCCTCGCCTCTGACGTGGTAATGGTTGCACGTGGCGACCTCGGTGTGGAAATTGGCGACCCGGAACTGGTCGGCATTCAGAAAGCGTTGATCCGTCGTGCGCGTCAGCTAAACCGAGCGGTAATCACGGCGACCCAGATGATGGAGTCAATGATTACTAACCCGATGCCGACGCGTGCA
+SRR587217.4 4 length=502
BTT^^^aaccccccccffffeeggggghgggegghffghggggghhhhhhhhhhgghhhegegggdgghhgggggggggcgggffffffffffffffffffffffffffffffffffbffffffffdffffffb`fffffffffffffffffffbfffffffffffffffffffffffdfffffffffffWYbddfdddffffffffffffdfffbffffffffcW]ffff]cfddddffffffffffffd^^^^^aa`ccccccccffffffhhhghghghggggggggghgbegggfgghgeegggggghhhhhhhhhhgfghghfggbcdggggg^cceeecegeeffffeeedfffaaffdfbdbddfdfffddffffffcff[`^bfdfffffWbffddbdfffdb]dddbbfdfcfffffSWbdd^bP^db]QWcccOSWbbbddb`ffW`fff^bbI^YYbbYYYOYY^^b^^bbddbbWW]bdbQMMW]cS]dW

Each sequence is divided into 4 lines. First line shows us the sequence name. Second line shows us the sequence itself. Third line shows us some comments (irrelevant for this programm) and the fourth line shows as the quality of the bases in ascii character (Phred score)

First of all i subdivide the three necessary lines (line1, line2 and line4)

Line one (sequence names) are saved in one list. Line two are saved in another list as well. In line four, I convert the ASCII character into numbers and save these numbers for each line in a different list. Example:

file=
"abcdefghijkl \n
"abcdefghijkl"
list of numbers=[[1,2,3,4,5,6,7,8][1,2,3,4,5,6,7,8]]

with those lists i calculate means, medians, variance for some plotting.

---

The code works fine with smaller "test.files" but if I use the big one it overloads my PC. Is there a way to compress the code? Should I use less lists/dicts?

def main():
    parser = ArgumentParser(prog='fastq', description=desc)
    parser.add_argument('dateien', metavar='Files', nargs='*',
                    default=['sys.stdin'], help='Input of the file')
    parser.add_argument('-p', '--phred',dest='phred',nargs=  '?',action='store',
                    choices= [33 ,64 ], required=True, type=int ) 
    parser.add_argument('-c', '--cutoff',dest = "cutoff",nargs='?',
                    action='store', default=1, type=int,
                    help='deletes whole Reads with bad quality')
    parser.add_argument('-t', '--trim',dest = 'trim',nargs='?',
                    action='store', default=1, type=int,
                    help='Trimming Readends with bad quality')
    #parser.add_argument('-v','--version',action='version', version="{prog}: {version}".format(prog=args.prog, version='0.1'))

    args = parser.parse_args()

    names = [] ; x=0 ;names1 = ""
    seq_dict={} ; mean_dict={} ; qual_dict={}
    qual_liste=[] ; i=0 ; sequenz=[] ; seq_name=[]
    mean1=[] ; summe=0
    quality_code=[]
    gc_dict={}
    qually=[]
    quality=[]
    with open('test.fastq', 'r') as seq:
      for line in seq:                
        quality_code=[]
        line=line.rstrip('\n')   
        rest=i%4                   #modulo to focus on specific lines
        i+=1                        
        #print(rest)

        if rest==0:     #line including the sequenz name
            name=line[:-13]    #some embellishment
            name1=name.strip("@")
            name2=name1.strip("length=" or "length" or "lengl" or "leng")
            name3=name2.rstrip()
            seq_name.append(name3) #List which includes lists of the seq names
        elif rest == 1:  #Line including sequenz nucleotides
            seq=line[:-args.trim] #trimming the sequenz
            sequenz.append(seq)
        elif rest ==3:    #line including quality score for sequenz
            qual=line[:-args.trim]           
            if args.phred == 33:
                quality_code=[[(ord(ii)-33) for ii in i] for i in qual.split('\n')]
            elif args.phred == 64:
                quality_code=[[(ord(ii)-64) for ii in i] for i in qual.split('\n')]
            for list in quality_code:
                mean=sum(list)/len(qual)              
            mean1.append(int(mean))    

            if mean >= int(args.cutoff):
                seq_dict[name3]=seq
                mean_dict[name3]=int(mean)   #filter/cutoff of reads with worse quality
                qual_dict[name3]=quality_code

The whole code includes some more lines about plotting to show the results.

Answer 1

This won't solve your problem, but here are some suggestions regarding the organization of the code:

• You can use a function like this to fetch four lines at once, which I think makes your code easier to understand than using rest:

  with open('test.fastq', 'r') as seq:
      for line1, line2, line3, line4 in group_it(4, seq):
          ...
  

• Instead of doing if args.phred == 33: you can just use the variable in the calculation, avoiding the repetition:

  quality_code=[[(ord(ii)-args.phred) for ii in i] for i in qual.split('\n')]