Cross-posted on Biostars (with no answers currently). I hope that's OK.
I thought that iff a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case.
This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.
What's wrong in my understanding?
Updated code (and updated output) moved to Gist.
Edit: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.
E00489:44:HNNVYCCXX:1:2120:26524:19135 163 5 345248 0 48M102S = 345446 350 TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)<)<--)<)--<-FAJFJFFFF-< XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M MD:Z:48 RG:Z:HNNVYCCXX.1 NM:i:0 MQ:i:60 AS:i:48 XS:i:48 E00489:44:HNNVYCCXX:1:2120:26524:19135 83 5 345446 60 6M2D144M = 345248 -350 TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG JJAF7-FJJJFJFF-AFJF7--<A--<F<7-7--FA--F-J--A-7--7-FJFJ7JFJJF<JA77-7----7----<F7-<JJF-<<FAJJJJA-7-<JAJJ<7FFJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA MC:Z:48M102S MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1 NM:i:7 MQ:i:0 AS:i:119 XS:i:69
I guess this explains why the read appears in more than
150 pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in
198 = (
144) pileup columns, at each position where one of them maps?
Update: I re-ran this script, using @Bioathlete's suggestion of using the flags to distinguish the read and its mate. In case
pileup() was filtering reads with low base quality, I also calculated how many positions the read has
baseQ >= 13. It's
150 positions, but this read still appears in pileup columns for only