Count the number of mapped and annotated reads in bam files

Question

[Katon]

Hi All,

Recently I am working on a fruit fly RNAseq project. After indexing and mapping with STAR, now I need to calculate the number of uniquely mapped and annotated reads as well as the number of uniquely mapped but unannotated reads in the bam file.

Does anyone know whether bedtools or other tools could be used to solve this problem?

Thanks in advance!

Answer 1

[h.mon]

STAR can output that for you, if you use the parameter --quantMode GeneCounts. You can also use featureCounts or HTSeq.

[GenoMax]

uniquely mapped but unannotated reads

Those are going to require bedtools/bedops.

[h.mon]

From STAR counts file, I think the N_noFeature is what he wants (head -n 4 sample1.ReadsPerGene.out.tab):

N_unmapped           764798         764798        764798
N_multimapping       814209         814209        814209
N_noFeature         8073899       24202648       8165278
N_ambiguous           91993            141         29605
MRPL17                 1123              1          1122
...

P.S:: the three counts columns correspond to unstranded, 1st read strand aligned with RNA (htseq-count option -s yes), and 2nd read strand aligned with RNA (htseq-count option -s reverse).

[Katon]

Thank you soooo much for your help!! I have reprocessed the STAR with quantMode and got the tab file. Just curious, could we also get the number of genes from coding/non coding regions respectively by dealing with these output data?

Thanks again!

[Katon]

Thank you for your suggestion! It seems bedtools need two files to do comparison. Since I have already have an A file. Do you have any idea how to build an annotated gene B file based on those fasta and GTF files?

Thanks again!