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.