<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>

CC BY-SA 4.0
2018-10-12T12:11:10.383Z
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.
mattm
CC BY-SA 4.0
2018-10-11T22:45:54.400Z
2018-10-11T22:54:02.530Z
2018-10-11T22:54:02.530Z
Karel Břinda
<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>

<sam><pysam>
How do I rewrite a read group using pysam?
# 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? 
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.  

---

