1
$\begingroup$

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)```
$\endgroup$
0
1
$\begingroup$

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

read_file = "/file/here.fq"
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)
    
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.