Read counts IGV vs SamTools for intronic region

Question

[mmagnus]

I'm trying to get read counts of a region that I see in IGV has around 2k read counts. When I use samtools I get 10x more. The number that it's very similar to the level of read counts for exons in this region.

How to get a number of reads for this intronic region from the Terminal?

 ~ samtools view -c MEX_WT_7Aligned.sortedByCoord.out.bam VII:439190-439190
24029

[Pierre Lindenbaum]

IGV using default parameters removes supplementary reads, duplicates, failing QC, etc... check this using the -f/-F option of samtools view

Answer 1

[mmagnus]

OK, I know what's going on. This bam file was produced by Star. And it seems that when we ask for, even 1bp, localization inside the intron (marked with ! ), samtools gives read counts for intron itself (A) AND also to splice junctions = reads that are spanned from one exon to another (B). Now this seems to be kind of obvious.

So we can use samtools to get only reads mapped to intron:

$ samtools view MEX_WT_7Aligned.sortedByCoord.out.bam VII:439190-439190  | cut -f 6 | grep -v 'N' | wc -l
2721