Count number of mapped raw reads for each transcript in SAM file

Question

[manekineko]

Is there a better war to count how many reads are mapped to each sequence/transcript and get a count table with names and counts than this:

perl -ne ' if (/^\@SQ/) { @F = split(/\t|:/, $_); print $F[2]."\n" } ' SAMFILE > ID.txt
perl -ne ' chomp($_); print $_."\t".`grep -c "\t$_" SAMFILE ` ' ID.txt > COUNTTABLE

The sam file is 8GB nevertheless the transcripts bowtie index was containing only 70 transcripts I need, and it takes me forever (couple of hours) to get this count-table?

[swbarnes2]

What about samtools idxstats?

[Ashutosh Pandey]

I think this should work for him.

[manekineko]

AFAIK the numbers reported by samtools idxstats represent the number of alignments of reads that are mapped to the sequence, not the (non-redundant) number of reads what I need?

EDIT:that worked fine :)

[darxsys]

I wrote an implementation in C++ for my purposes using seqan library to do just that. Works fairly fast. Didn't try it on such big files though. You can try doing something similar. With seqan, it should be relatively painless.

[Ashutosh Pandey]

Why you don't use standard tools such as HTSeq, summarizeOverlaps in R or RSeQC?

Answer 1

[Brian Bushnell]

With the BBMap package:

pileup.sh in=mapped.sam covstats=stats.txt rpkm=rpkm.txt

It's fast.

Answer 2

[Antonio R. Franco]

How about using a gff or gft file of your reference genome and a dedicated program such as:

• summarizeOverlaps from the R package GenomicAlignments
• htseq-count which is written in phyton HTSeq (Python)
• featureCounts from the R package Rsubread
• simpleRNASeq from the R package easyRNASeq

[Ashutosh Pandey]

Another one http://rseqc.sourceforge.net/

[manekineko]

The problem is that I do not want to make it complicated involving whole genome or transcriptome analysys, I just have 70 sequences which I index for bowtie and map my reads to a SAM file....and want to count the reads in those 70 sequences.

[Ashutosh Pandey]

Check swbarnes2's comment above.