Hi, Question about counting multimapped reads, things like pirna and tRNA. For instance, I ran star aligning small rna mouse reads to mm10. Performed counting in htseq-count as well as in featurecounts using -M command, but still don't seem to get counts for these reads..just indication that I have a ton of non-unique alignments (expected given interest in tRNA). Any suggestions for counting trna/non unique reads? And quick simple question...how do I get ncrna features annotated in my counts.txt (e.g. miRNA, tRNA)? Only info in my output is ensembl number and chromosome.
Thanks a bunch, Greg
Second
1 answer
I made a pipeline adapted to small RNA-seq with multimapped reads. It is based on bowtie2, allowing reads to map at up to 10 locations (option -k 10), but keeping only best-scoring ex-aequo. If I remember well, I divided the contribution of a given read to the total counts by the number of such ex-aequo locations. The counting is done in a custom python script based on pysam.
The pipeline is here: https://bitbucket.org/blaiseli/pirna-pipeline
It is for Drosophila and we were not interested in tRNAs, so the counting of tRNAs is only global, not at the gene level. We were more interested in piRNAs. We counted them in a second step by re-mapping the reads identified as piRNA candidates to piRNA clusters and transposable elements. Besides, I was then unaware of lots of tools, formats and ways of doing common in bioinformatics, so my approach may sound complicated. It would probably be simpler to write your own counting script, but feel free to fork, adapt, use ideas...
Log in to answer this question.