# How do I rewrite a read group using pysam?

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)


## 1 Answer

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.

• 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. Oct 12 '18 at 12:11