Krdytkyu

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

file=@SRR587217.1 1 length=490

NTCCGGATGATGTCGCTGTTGCTGACAATGGTAATACGTTGACGGGGCAATATGCAGTTCGCTGCATACCGGTCCGACCCCGTACTGCTCACGCAGCTTATCCAGCAGTGGCATCATTTTTTCCAGAGGCGGTCGAACTCCGCCTTCGCAAAAAAAAAGGGAGCCCGGCGGAGGAGAACGTTACTGCGGCGGAGGTTACGATTTTTCCGGTTCCGCTCCTTTAGAAGCCGGACGTCTACCCGGCTCTTTTTGCTGAACGTCAGCGTCTGAAAGAGCTGGAACGTGAAAATCGTGAACTGCGCCGCAGTAACGATATCCTTCGCCAGGCTTCCGCTTATTTTGCGAAGGCGGAGTTCGACCGCCTCTGGAAAAAATGATGCCACTGCTGGATAAGCTGCGTGAGCAGTACGGGGTCGGACCGGTATGCAGCGAACTGCATATTGCCCCGTCAACGTATTACCATTGTCAGCAACAGCGACATCATCCGGAT

+SRR587217.1 1 length=490

BT[^^^^^caccccacffffffhhhgghhhhefeggghhghghgfgggghfgggeghhhhhgggghhgggggggggggggggdgggggfffffffffffffddfffffbfffffffffffffffffffffffMFFFHOIIPHMWHOIIHHMM^IFTMFHOHIHMHWMWFMFFFTHHIIIHIOIHIIIHMFFFFFFOIIIHHIOH^IIIHMSFOIHMHMFPYIIIIIIIIFFMFMOWIIIPIHFFMOIYIII^^``^aaaccccccccffffffghghhggghghgghghhhggggefgfededbe`cgeghggegaefgghegggggggggccdggggceeggfT``dddfffdfffffac]Zbffdfdbdbdfffddd^^dfdfdfddffffdbdbdfffffffbdddf``c`bccQW`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.

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')]