Hi. I have been working with mouse RNA-Seq data and I aimed to use Rsubread package to align, annotate and quantify the reads. My issue is that whenever I use the inbuilt annotation for the mouse genome practically none of the reads are mappped:

align (index = "GRCm38.p4", readfile1 = all.fastq1, readfile2 = all.fastq2, input_format = "FASTQ", output_format="BAM", nthreads=8, output_file = all.BAM)
featuresC <- featureCounts(all.BAM, annot.inbuilt="mm10")

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file mm10_RefSeq_exon.txt ...                              ||
||    Features : 222996                                                       ||
||    Meta-features : 27179                                                   ||
||    Chromosomes/contigs : 43                                                ||
||                                                                            ||
|| Process BAM file 1A.bam...                                                 ||
||    Paired-end reads are included.                                          ||
||    Assign alignments to features...                                        ||
||    Total alignments : 93836764                                             ||
||    Successfully assigned alignments : 29766 (0,0%)                         ||
||    Running time : 3,96 minutes                                             ||
||                                                                            ||

I tried with the latest patch too (p6) and I obtained the same result. It only works by explicitly stating an external annotation file (.gtf). But in that case genes are being annotated by name and not by entrezid which I require for gene ontology analysis and so on.

Could somebody help me with this? Thank you very much in advance.

The problem is that NCBI has renamed all the chromosomes in its recent genome releases. A few years ago, chromosome 1 was called "chr1" in NCBI's fasta files, but it is now called "NC_000067.6". The Rsubread in-built annotation uses "chr1", "chr2" etc, which is compatible with UCSC and Gencode but no longer with NCBI.

Note that it's not a matter of which patch (p1, p4, p6) was used, as all the reference chromosomes are unchanged between patches. It's just an arbitrary change of notation by NCBI.

The easiest solution is to use a NCBI fasta file from before they renamed the chromosomes and gene ids. That's what I do, and it would seem that James MacDonald is doing the same.

Your cutting'n'pasting seems off - you are running featureCounts on 'all.BAM' but featureCounts appears to be running on 1A.bam.

Anyway, do you have a mismatch between the genomes you are using? The inbuilt annotation is based on the NCBI genome, which has chr1, chr2, etc. If you used an Ensembl genome to align, then all the chromosomes will be 1, 2, 3, etc. You can check using Rsamtools, doing something like

head(scanBamHeader("all.BAM")$targets)

Where the names of the resulting vector are the chromosome names:

> bf <- BamFile("../data/aligned_nodups/6-E-4Aligned.sortedByCoord.out.bam")
> head(scanBamHeader(bf)$targets)
     chr1     chr10     chr11     chr12     chr13     chr14 
195471971 130694993 122082543 120129022 120421639 124902244 

And my BAM file is aligned against the mm10 genome from NCBI.