This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Read counts IGV vs SamTools for intronic region

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

rna-seq
• 3,808 views • link
written by mmagnus mmagnus

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

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

1 answer

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

enter image description here

• 0 views • link
written by mmagnus mmagnus

Log in to answer this question.

More posts like this