Answers.Body: <p>For can do this by accessing the basecall qualities from <a href="https://pysam.readthedocs.io/en/latest/api.html#pysam.PileupRead.alignment" rel="nofollow noreferrer"><code>PileupRead.alignment</code></a>. For example:</p>
  
  <pre><code>_BASE_QUAL_CUTOFF = 30
  with pysam.Samfile(in_file, 'rb') as bam_in:
      for column in bam_in.pileup():
          for read in column:
              qual = read.alignment.query_qualities
              if all([ord(c)-33 &gt;= _BASE_QUAL_CUTOFF for c in qual]):
                  # do some stuff
  </code></pre>
  
  <p>This involves checking every read at every position. You could memoise this task:</p>
  
  <pre><code>def check_quality(read):
      try:
          return passed_quality[read.query_name]
      except KeyError:
          qual = read.alignment.query_qualities 
          passed = all([ord(c)-33 &gt;= _BASE_QUAL_CUTOFF for c in qual])
          passed_quality[read.query_name] = passed
          return passed_quality[read.query_name]
  
  _BASE_QUAL_CUTOFF = 30
  passed_quality = dict()
  with pysam.Samfile(in_file, 'rb') as bam_in:
      for column in bam_in.pileup():
          for read in column:
              if check_quality(read):
                  # do some stuff
  </code></pre>
  
  <p>Or, using <code>functools.lru_cache</code> (thanks @winni2k):</p>
  
  <pre><code>from functools import lru_cache
  
  @lru_cache(maxsize=1000)
  def check_quality(read):
      qual = read.alignment.query_qualities 
      return all([ord(c)-33 &gt;= _BASE_QUAL_CUTOFF for c in qual])
  
  _BASE_QUAL_CUTOFF = 30
  passed_quality = dict()
  with pysam.Samfile(in_file, 'rb') as bam_in:
      for column in bam_in.pileup():
          for read in column:
              if check_quality(read):
                  # do some stuff
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-23T09:30:24.017Z
Answers.Text: FYI: There's a decorator in the library that does this and more! https://docs.python.org/3/library/functools.html#functools.lru_cache
Answers.UserDisplayName: winni2k
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-21T18:12:16.310Z
Answers.LastActivityDate: 2018-11-24T18:25:05.413Z
Answers.LastEditDate: 2018-11-24T18:25:05.413Z
Answers.UserDisplayName: Scott Gigante
Body: <p>I have a code below:</p>
  
  <pre><code>def FilterReads(in_file, out_file):
  
  def read_ok(read):
      """
      read_ok - reject reads with a low quality (&lt;5) base call
      read - a PySam AlignedRead object
      returns: True if the read is ok
      """
      if any([ord(c)-33 &lt; _BASE_QUAL_CUTOFF for c in list(read.qual)]):
          return False
      else:
          return True
  
  _BASE_QUAL_CUTOFF = 30
  
  bam_in = pysam.Samfile(in_file, 'rb')
  bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)
  
  
  for read in bam_in.fetch():
      if read_ok(read):
          bam_out.write(read)
  </code></pre>
  
  <p>This code works fine by </p>
  
  <ul>
  <li>Rejecting reads having a base with phred quality score below five</li>
  <li>But it first takes a BAM file </li>
  <li>And then creates a filtered BAM file for further analysis but this takes a lot of time. </li>
  </ul>
  
  <p>So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it? 
  Or can I modify <a href="https://bioinformatics.stackexchange.com/questions/2131/filtering-bases-based-on-phred-qualities-with-pysam">this code</a> to perform the same function</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-31T21:11:55.713Z
Comments.Text: Why are you sorting the file prior to iterating over it?  Your function doesn't require it to be sorted and since you are not using the random indexing pysam doesn't require it either.
Comments.UserDisplayName: Bioathlete
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-24T12:20:38.803Z
Comments.Text: `if …: return False else: return True` is an anti-pattern. Just write `return not …`. Or, better yet, logically invert the condition; in your case, you could write `return all([ord(c) - 33 >= _BASE_QUAL_CUTOFF for c in list(read.qual)])` (this follows from straightforward application of [de Morgan’s laws](https://en.wikipedia.org/wiki/De_Morgan%27s_laws)).
Comments.UserDisplayName: Konrad Rudolph
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-23T09:26:48.170Z
Comments.Text: This is probably not the cause of the slowdown, but you can avoid two list creations per read by refactoring line 9 as `any(ord(c)-33 < _BASE_QUAL_CUTOFF for c in read.qual):`
Comments.UserDisplayName: winni2k
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-23T09:32:47.557Z
Comments.Text: I don't think pileup is going to save you any time here, as it has to load the same data that you are loading with the `fetch` method.
Comments.UserDisplayName: winni2k
Tags: <sam><ngs><phylogenetics><python><pysam>
Title: Reject reads with low quality bases from a Bam file through pysam
ThreadText: # Reject reads with low quality bases from a Bam file through pysam
  
  (Asked by: username_0 on 2017-07-26)
  
  I have a code below:
  
  ```
  def FilterReads(in_file, out_file):
  
  def read_ok(read):
      """
      read_ok - reject reads with a low quality (<5) base call
      read - a PySam AlignedRead object
      returns: True if the read is ok
      """
      if any([ord(c)-33 < _BASE_QUAL_CUTOFF for c in list(read.qual)]):
          return False
      else:
          return True
  
  _BASE_QUAL_CUTOFF = 30
  
  bam_in = pysam.Samfile(in_file, 'rb')
  bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)
  
  for read in bam_in.fetch():
      if read_ok(read):
          bam_out.write(read)
  
  ```
  
  This code works fine by
  
  * Rejecting reads having a base with phred quality score below five
  * But it first takes a BAM file
  * And then creates a filtered BAM file for further analysis but this takes a lot of time.
  
  So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it?
  Or can I modify [this code](https://bioinformatics.stackexchange.com/questions/2131/filtering-bases-based-on-phred-qualities-with-pysam) to perform the same function
  
  **Tags:** sam, ngs, phylogenetics, python, pysam
  
  > **Comments:**
  > **username_1**: Why are you sorting the file prior to iterating over it?  Your function doesn't require it to be sorted and since you are not using the random indexing pysam doesn't require it either.  
  > **username_2**: `if …: return False else: return True` is an anti-pattern. Just write `return not …`. Or, better yet, logically invert the condition; in your case, you could write `return all([ord(c) - 33 >= _BASE_QUAL_CUTOFF for c in list(read.qual)])` (this follows from straightforward application of [de Morgan’s laws](https://en.wikipedia.org/wiki/De_Morgan%27s_laws)).  
  > **username_3**: This is probably not the cause of the slowdown, but you can avoid two list creations per read by refactoring line 9 as `any(ord(c)-33 < _BASE_QUAL_CUTOFF for c in read.qual):`  
  > **username_3**: I don't think pileup is going to save you any time here, as it has to load the same data that you are loading with the `fetch` method.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_4 on 2018-10-21. Score: 1
  
  For can do this by accessing the basecall qualities from [`PileupRead.alignment`](https://pysam.readthedocs.io/en/latest/api.html#pysam.PileupRead.alignment). For example:
  
  ```
  _BASE_QUAL_CUTOFF = 30
  with pysam.Samfile(in_file, 'rb') as bam_in:
      for column in bam_in.pileup():
          for read in column:
              qual = read.alignment.query_qualities
              if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                  # do some stuff
  
  ```
  
  This involves checking every read at every position. You could memoise this task:
  
  ```
  def check_quality(read):
      try:
          return passed_quality[read.query_name]
      except KeyError:
          qual = read.alignment.query_qualities 
          passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
          passed_quality[read.query_name] = passed
          return passed_quality[read.query_name]
  
  _BASE_QUAL_CUTOFF = 30
  passed_quality = dict()
  with pysam.Samfile(in_file, 'rb') as bam_in:
      for column in bam_in.pileup():
          for read in column:
              if check_quality(read):
                  # do some stuff
  
  ```
  
  Or, using `functools.lru_cache` (thanks @username_3):
  
  ```
  from functools import lru_cache
  
  @lru_cache(maxsize=1000)
  def check_quality(read):
      qual = read.alignment.query_qualities 
      return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
  
  _BASE_QUAL_CUTOFF = 30
  passed_quality = dict()
  with pysam.Samfile(in_file, 'rb') as bam_in:
      for column in bam_in.pileup():
          for read in column:
              if check_quality(read):
                  # do some stuff
  
  ```
  
  > **Comments:**
  > **username_3**: FYI: There's a decorator in the library that does this and more! https://docs.python.org/3/library/functools.html#functools.lru_cache  
  
  ---
  
  