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)