Get per-nucleotide mapping position for each mapped reads from bam file

Question

[Changhao Li]

Hello!

I'm trying to get per-nucleotide mapping position for each mapped read from a bam file.

Specifically, for each read from a bam file, I want to know the mapping position for the 13th nt in the read. Because its RNAseq bam, there should be lots of splicing reads. So we cannot just count 13 from the first mapping position.

So does anyone have any scripts to solve this problem? Many thanks in advance.

Answer 1

[Pierre Lindenbaum]

I wrote sam2tsv : http://lindenb.github.io/jvarkit/Sam2Tsv.html

$  java -jar dist/sam2tsv.jar -R src/test/resources/toy.fa src/test/resources/toy.bam 


#Read-Name  Flag  MAPQ  CHROM  READ-POS0  READ-BASE  READ-QUAL  REF-POS1  REF-BASE  CIGAR-OP
r001        163   30    ref    0          T          .          7         T         M
r001        163   30    ref    1          T          .          8         T         M
r001        163   30    ref    2          A          .          9         A         M
r001        163   30    ref    3          G          .          10        G         M
r001        163   30    ref    4          A          .          11        A         M
r001        163   30    ref    5          T          .          12        T         M
r001        163   30    ref    6          A          .          13        A         M
r001        163   30    ref    7          A          .          14        A         M
r001        163   30    ref    8          A          .          .         .         I
r001        163   30    ref    9          G          .          .         .         I
r001        163   30    ref    10         A          .          .         .         I
r001        163   30    ref    11         G          .          .         .         I

[Changhao Li]

Thank you very much Pierre! I will try your method.

[Changhao Li]

Hi Pierre,

I have one more question is that do you consider strandness in your scripts. For example, if one read mapped to the negative strand, the 5' of the read should map to the 3' of the reference.

[Pierre Lindenbaum]

have one more question is that do you consider strandness in your scripts. For example, if one read mapped to the negative strand, the 5' of the read should map to the 3' of the reference.

no, it's always from 5' to 3' of the REF.

[Changhao Li]

So for the negative mapped reads, if the READ-POS0 is 0, actually it should be the last nucleotide, right?

[Pierre Lindenbaum]

So for the negative mapped reads, if the READ-POS0 is 0, actually it should be the last nucleotide, right?

yes !

[Changhao Li]

OK, thank you very much!