How to edit a SAM file using pysam
0
0
Entering edit mode
2.1 years ago
Apex92 ▴ 280

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 • 608 views
ADD COMMENT

Login before adding your answer.

Traffic: 2289 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6