pysam rewrite help

Question

[Anst]

Hello!

I am trying to rewrite the existing bam file by replacing some nucleotides in given reads. For example

ATCCG -> TTCCG

I though about a prototype of if condition like

str1 = 'ATCCG'

if str1[0] == 'A' : str2[0] = 'T'

Although, did not find a good explaination in documentation how to do it. Maybe there is more accurate way of rewriting file via pysam.

Thanks!

[Pierre Lindenbaum]

see How to introduce artificial mutation in bam ; Modifying fastq base at specific reference location on different length reads ;

Answer 1

[dariober]

Since you have pysam in the tile, maybe something along these lines?

import pysam

samfile = pysam.AlignmentFile("in.bam", "rb")
outsam = pysam.AlignmentFile('out.bam', 'wb', template= samfile)

for aln in samfile:
    if aln.query_sequence == 'ATCCG':
        q = aln.query_qualities
        aln.query_sequence = 'TTCCG'
        aln.query_qualities = q
    outsam.write(aln)

samfile.close()
outsam.close()

[Anst]

Thanks! But how I can interate though the nucleotides? For example, I'd like to change each third nucleotide to the custom one?

I thought about:

for pos in read_positions:
    if reference_nucl == 'A':
       read_nucl = 'C'

but I don't know how to organise it in pysam