Read counts IGV vs SamTools for intronic region
1
0
Entering edit mode
5.6 years ago
mmagnus ▴ 150

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.3k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
5.6 years ago
mmagnus ▴ 150

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

ADD COMMENT

Login before adding your answer.

Traffic: 1756 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6