This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to edit a SAM file using pysam

Dear all - I have a template sam file and I want to change one of the columns (template_length) and replace it with a new value. The new value is a quick mathematical operation.

template sam file:

@HD VN:1.0  SO:unsorted
@SQ SN:test_ref LN:17637
SRR1524970.144283   16  test_ref    1706    255 25M *   0   0   TGCTGATGAAGCAGAACAACTTTAA   ]YG[^baaaa^W`ab]]````aaba   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:25 YT:Z:UU
SRR1524970.316478   16  test_ref    1706    255 24M *   0   0   TGCTGATGAAGCAGAACAACTTTA    `\X_`aaaaaY]``b_aa_aaaaa    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:24 YT:Z:UU

I use pysam to make this change but based on my code it outputs exact columns of the input sam file with =out changing the wanted column - am I missing something?

my code:

import pysam

genome_length = 400

samfile = pysam.AlignmentFile("test.sam", "r")

with pysam.AlignmentFile("test_new.sam", "w", template=samfile) as outputsam:
    for i in samfile.fetch():
        a = pysam.AlignedSegment(outputsam.header)
        a.query_name = i.query_name
        a.query_sequence = i.query_sequence
        a.reference_name = i.reference_name
        a.flag = i.flag
        a.reference_start = i.reference_start
        a.mapping_quality = i.mapping_quality
        a.cigar = i.cigar
        a.next_reference_id = i.next_reference_id
        a.next_reference_start = i.next_reference_start
        a.template_length = genome_length/pow(4,len(i.query_sequence))
        a.query_qualities = i.query_qualities
        a.tags = i.tags
        outputsam.write(a)

samfile.close()
outputsam.close()

Thank you in advance.

sam pysam mapping algnment bam

0 answers

No answers yet.

Log in to answer this question.