Hello everyone!
I am trying to rewrite an existing bam file via pysam. I have reference.fasta file and file_to_be_changed.bam. I would like to change the nucleotide on a read given the position of the nucleotide and the reference position.
For example:
pos 0123456
ref ATCGGTC
bam ATCG
CGT
GGTC
I think that the problem could be solved with:
if pos_on_read == 0:
if nucl_on_reference == 'A':
nucl_on_read = 'C'
However, I have some problems with realising the algorithm via pysam.. Give me a hand please!
•
written
by
Anst •
1 answer
I would say that the problem you pose is be a lot more difficult to do than you think.
You could iterate for each alignment and you will get the POS field of each, then as you iterate on the read (on the aligned region) you could modify the stored sequence.
Then you could write out the alignment into a new BAM file.
But frankly, that alignment would be incorrect, the score would be wrong, the mapping quality would be incorrect, the CIGAR string would be incorrect, some of the TAGS would be incorrect, and even the alignment itself may be invalid, the modified read may then align somewhere else better.
I can't help by wonder what is the purpose of this fix? Your goal sounds a bit like going after a problem with a premature fix. Perhaps ask a different question on the problem you are after then trying to modify a BAM file.
But then some of the alignments would not be correct. How come that is not a concern? To have incorrect alignments?
An alignment operates in the context of the entire information in the genome, you can't expect that locally changing both the reference and read would still produce valid alignments across the genome. Sure, some, or even the majority might stay correct, but that ratio depends on many factors and the error rate would be highly uneven, might be the most error-prone in the locations where the information between ancient and non-ancient genomes is more meaningful.
that would work
you could extract your reads from the BAM file then realign them to the new genome
Log in to answer this question.
