# Program to process fastq files

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^^^aaaccccccccffffffghghhggghghgghghhhggggefgfededbecgeghggegaefgghegggggggggccdggggceeggfTdddfffdfffffac]Zbffdfdbdbdfffddd^^dfdfdfddffffdbdbdfffffffbdddfcb[]ccQWMSHYbdffff]cbbdbbbdYI^bd^[WdbbbfMXbbbf^b   ^^b^bdbdddIX^d]W]QW^dd^^dcfW
@SRR587217.4 4 length=502
NTCCATGACTTCTGCACGCGTCGGCATCGGGTTAGTAATCATTGACTCCATCATCTGGGTCGCCGTGATTACCGCTCGGTTTAGCTGACGCGCACGACGGATCAACGCTTTCTGAATGCCGACCAGTTCCGGGTCGCCAATTTCCACACCGAGGTCGCCACGTGCAACCATTACCACGTCAGAGGCGAGGATGATGTCATCCATTGCATCCTGGCTGCAAACGGCTTCCGCACGTTCAACCTTGGCAACAAGTGATGCGAAAATTGTTGCCAAGGTTGAACGTGCGGAAGCCGTTTGCAGCCAGGATGCAATGGATGACATCATCCTCGCCTCTGACGTGGTAATGGTTGCACGTGGCGACCTCGGTGTGGAAATTGGCGACCCGGAACTGGTCGGCATTCAGAAAGCGTTGATCCGTCGTGCGCGTCAGCTAAACCGAGCGGTAATCACGGCGACCCAGATGATGGAGTCAATGATTACTAACCCGATGCCGACGCGTGCA
+SRR587217.4 4 length=502
BTT^^^aaccccccccffffeeggggghgggegghffghggggghhhhhhhhhhgghhhegegggdgghhgggggggggcgggffffffffffffffffffffffffffffffffffbffffffffdffffffbfffffffffffffffffffbfffffffffffffffffffffffdfffffffffffWYbddfdddffffffffffffdfffbffffffffcW]ffff]cfddddffffffffffffd^^^^^aaccccccccffffffhhhghghghggggggggghgbegggfgghgeegggggghhhhhhhhhhgfghghfggbcdggggg^cceeecegeeffffeeedfffaaffdfbdbddfdfffddffffffcff[^bfdfffffWbffddbdfffdb]dddbbfdfcfffffSWbdd^bP^db]QWcccOSWbbbddbffWfff^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)
default=['sys.stdin'], help='Input of the file')
choices= [33 ,64 ], required=True, type=int )
action='store', default=1, type=int,
action='store', default=1, type=int,

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.

• What do you use 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. – André Paramés Feb 15 '18 at 9:13
• You might want to look into this project: github.com/JohnLonginotto/SeQC – Austin Hastings Mar 17 '18 at 12:37

• 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:

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