Answers.Body: <p>If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: <a href="https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups" rel="nofollow noreferrer">https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups</a>.</p>
  
  <p>Modifying headers using PySAM might be tricky – the <a href="https://github.com/pysam-developers/pysam/issues/618" rel="nofollow noreferrer">internal representation of headers</a> has recently been changed and seems not to be stable across PySAM versions.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-12T12:11:10.383Z
Answers.Text: I want to assign all read groups to the same sample, where read groups are from different flowcells, but not rewrite all alignments to the same read group, so the Picard AddOrReplaceReadGroups is not the tool I am looking for. I have actually now written something using htsjdk in java to perform the operation I want, but I will leave the question for pysam.
Answers.UserDisplayName: mattm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-11T22:45:54.400Z
Answers.LastActivityDate: 2018-10-11T22:54:02.530Z
Answers.LastEditDate: 2018-10-11T22:54:02.530Z
Answers.UserDisplayName: Karel Břinda
Body: <p>I am trying to rewrite a SAM/BAM file with altered read group entries using pysam. In this simplified version, I want to take a BAM, and rewrite the SM tags in all read groups to the same string and copy the alignments with new header into a new BAM. The contents of the read group dictionaries I am changing appear to have no effect on the header that prints or is written to file output. How can I alter the read groups? </p>
  
  <pre class="lang-python prettyprint-override"><code>import argparse
  import pysam
  
  def rewrite_read_groups(input_filename, output_filename, sample):
      try:
          bam_input = pysam.AlignmentFile(input_filename, 'rb')
          # replace header fields
          output_header = bam_input.header
          read_groups = output_header['RG']
          print(read_groups)
          for read_group in read_groups:
              read_group['SM'] = sample
              print(read_group)
          output_header.to_dict()
          print(output_header)
      finally:
          bam_input.close()
  
      print(output_header)
  
      try:
          bam_input = pysam.AlignmentFile(input_filename, 'rb')
          bam_output = pysam.AlignmentFile(output_filename, 'wb', header=output_header)
  
          # copy alignments
          for read in bam_input.fetch(until_eof=True):
              bam_output.write(read)
  
      finally:
          bam_input.close()
          bam_output.close()
  
  if __name__ == "__main__":
      parser = argparse.ArgumentParser(description="Rewrite bam read groups")
      parser.add_argument('-s', "--sample", help="Set sample for all read groups")
  
      parser.add_argument('-i', "--input_file", help="Input bam", required=True)
      parser.add_argument('-o', "--output_file", help="Output bam", required=True)
  
      args = parser.parse_args() 
  
      rewrite_read_groups(args.input_file, args.output_file, args.sample)
  </code></pre>
  
Tags: <sam><pysam>
Title: How do I rewrite a read group using pysam?
ThreadText: # How do I rewrite a read group using pysam?
  
  (Asked by: username_0 on 2018-10-09)
  
  I am trying to rewrite a SAM/BAM file with altered read group entries using pysam. In this simplified version, I want to take a BAM, and rewrite the SM tags in all read groups to the same string and copy the alignments with new header into a new BAM. The contents of the read group dictionaries I am changing appear to have no effect on the header that prints or is written to file output. How can I alter the read groups?
  
  ```python
  import argparse
  import pysam
  
  def rewrite_read_groups(input_filename, output_filename, sample):
      try:
          bam_input = pysam.AlignmentFile(input_filename, 'rb')
          # replace header fields
          output_header = bam_input.header
          read_groups = output_header['RG']
          print(read_groups)
          for read_group in read_groups:
              read_group['SM'] = sample
              print(read_group)
          output_header.to_dict()
          print(output_header)
      finally:
          bam_input.close()
  
      print(output_header)
  
      try:
          bam_input = pysam.AlignmentFile(input_filename, 'rb')
          bam_output = pysam.AlignmentFile(output_filename, 'wb', header=output_header)
  
          # copy alignments
          for read in bam_input.fetch(until_eof=True):
              bam_output.write(read)
  
      finally:
          bam_input.close()
          bam_output.close()
  
  if __name__ == "__main__":
      parser = argparse.ArgumentParser(description="Rewrite bam read groups")
      parser.add_argument('-s', "--sample", help="Set sample for all read groups")
  
      parser.add_argument('-i', "--input_file", help="Input bam", required=True)
      parser.add_argument('-o', "--output_file", help="Output bam", required=True)
  
      args = parser.parse_args() 
  
      rewrite_read_groups(args.input_file, args.output_file, args.sample)
  
  ```
  
  **Tags:** sam, pysam
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2018-10-11. Score: 2
  
  If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: <https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups>.
  
  Modifying headers using PySAM might be tricky – the [internal representation of headers](https://github.com/pysam-developers/pysam/issues/618) has recently been changed and seems not to be stable across PySAM versions.
  
  > **Comments:**
  > **username_0**: I want to assign all read groups to the same sample, where read groups are from different flowcells, but not rewrite all alignments to the same read group, so the Picard AddOrReplaceReadGroups is not the tool I am looking for. I have actually now written something using htsjdk in java to perform the operation I want, but I will leave the question for pysam.  
  
  ---
  
  