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.
sequenz
for? You are storing every second line of the file in a list; if the file has 1.5GB, that means using hundreds of MBs of memory even if Python stored them very efficiently. It's also not very CPU friendly, since the list is constantly having to be expanded. \$\endgroup\$