# 10bp sliding window quality score

I am having trouble creating a 10 basepair sliding window that goes through seq and trim them all 10 basepairs that have an average quality score (qual) less than 15. I am trying to approach this by reading in the file and seperating the sequences and quality scores, then read in the 10 basepair window in the sequence.

#!/usr/bin/envpython
with open('fastaQ.fq') as f:
lines = [line.rstrip() for line in f]
seq= []
qual= []
tenFrame= []

#Seperates the sequences and the quality score of the file
for i in range(len(lines)):
if lines[i].startswith('@'):
seq.append(lines[i+1])
qual.append(lines[i+3])
i=i+4

#
for i in range(0,len(seq)):
tenFrame.append(seq[i][i:i+10])
if len(tenFrame)%10 == 0:
print(tenFrame)$$$$


I'm not sure what you mean by "trim". The code below will discard any reads that have a 10bp window with a mean phred score below your cutoff. I'll ammend the post to also include a solution that "trims" basepairs on both ends until the window achieves an average score of above the cutoff. It would be ideal if you could provide a minimal reproducible example that shows what your code does and what you'd like it to do.

As far as the code goes, I'd first suggest using a library like Biopython to parse the Fastq file. It will make your life a lot easier. For your code, it looks like seq actually contains all the sequences, so len(seq) is the number of reads in the file. However, when you iterate over i in range(len(seq)), it looks like you're using i to both index which read to iterate over as well as get the 10 base pair window. You need an additional index variable to obtain the window:

for window_start in range(len(seq[i]) - window_len) + 1)
window = seq[i][window_start:window_start + window_len]


Additionally, some pythonic tips for the future:

• If a range starts from 0, don't specify it as an arugment i.e. use range(x) as opposed to range(0, x)
• Try to avoid anything along the lines of range(len(x)). In this case, you could just iterate over the sequences directly: for s in seq: or, if you wanted the index, you could enumerate the sequences via for seq_idx, s in enumerate(seq)
• Use good variable names. seq implies a single sequence, yet is a list of sequences. Same with qual

Here is some code that will retain only sequences that have all windows with a passing mean phred score

from Bio import SeqIO

w = 10
cutoff = 15
valid_seqs = []
with open(fq_file) as fq_fd:
seqs = SeqIO.parse(fq_fd, "fastq")
for seq in seqs:
scores = seq.letter_annotations["phred_quality"]
mean_scores = (np.mean(scores[i:i+w]) for i in range(len(scores) - w))
if any(mean_score < cutoff for mean_score in mean_scores):
print(f"{seq.id} does not pass the filter!")
continue
valid_seqs.append(seq)

`