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?
•
written
by
manekineko
With the BBMap package:
pileup.sh in=mapped.sam covstats=stats.txt rpkm=rpkm.txt
It's fast.
•
written
by
Brian Bushnell
How about using a gff or gft file of your reference genome and a dedicated program such as:
summarizeOverlaps from the R package GenomicAlignments
•
written
by
Antonio R. Franco
Another one http://rseqc.sourceforge.net/
•
written
by
Ashutosh Pandey
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.
Check swbarnes2's comment above.
Log in to answer this question.
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.
Why you don't use standard tools such as HTSeq, summarizeOverlaps in R or RSeQC?
What about samtools idxstats?
I think this should work for him.
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 :)