Answers.Body: <p>I'm not sure what you mean by &quot;trim&quot;. 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 &quot;trims&quot; 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.</p>
  <p>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 <code>seq</code> actually contains all the sequences, so <code>len(seq)</code> is the number of reads in the file. However, when you iterate over <code>i in range(len(seq))</code>, it looks like you're using <code>i</code> 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:</p>
  <pre><code>for window_start in range(len(seq[i]) - window_len) + 1)
      window = seq[i][window_start:window_start + window_len]
  </code></pre>
  <p>Additionally, some pythonic tips for the future:</p>
  <ul>
  <li>If a range starts from 0, don't specify it as an arugment i.e. use <code>range(x)</code> as opposed to <code>range(0, x)</code></li>
  <li>Try to avoid anything along the lines of <code>range(len(x))</code>. In this case, you could just iterate over the sequences directly: <code>for s in seq:</code> or, if you wanted the index, you could enumerate the sequences via <code>for seq_idx, s in enumerate(seq)</code></li>
  <li>Use good variable names. <code>seq</code> implies a single sequence, yet is a list of sequences. Same with <code>qual</code></li>
  </ul>
  <p>Here is some code that will retain only sequences that have all windows with a passing mean phred score</p>
  <pre><code>from Bio import SeqIO
  
  read_file = &quot;/file/here.fq&quot;
  w = 10
  cutoff = 15
  valid_seqs = []
  with open(fq_file) as fq_fd:
      seqs = SeqIO.parse(fq_fd, &quot;fastq&quot;)
      for seq in seqs:
          scores = seq.letter_annotations[&quot;phred_quality&quot;]
          mean_scores = (np.mean(scores[i:i+w]) for i in range(len(scores) - w))
          if any(mean_score &lt; cutoff for mean_score in mean_scores):
              print(f&quot;{seq.id} does not pass the filter!&quot;)
              continue
          valid_seqs.append(seq)
      
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-06T00:46:15.280Z
Answers.LastActivityDate: 2021-05-06T00:46:15.280Z
Answers.UserDisplayName: Throckmorton
Body: <p>I am having trouble creating a 10 basepair sliding window that goes through <code>seq</code> and trim them all 10 basepairs that have an average quality score (<code>qual</code>) 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.</p>
  <pre><code>#!/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)<span class="math-container">```</span>
  </code></pre>
  
Tags: <python>
Title: 10bp sliding window quality score
ThreadText: # 10bp sliding window quality score
  
  (Asked by: username_0 on 2021-05-05)
  
  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)```
  
  ```
  
  **Tags:** python
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-05-06. Score: 1
  
  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)
      
  
  ```
  
  
  ---
  
  