How do I rewrite a read group using pysam?

Question

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)

Answer 1

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 has recently been changed and seems not to be stable across PySAM versions.