[BioC] Rsubread featureCounts unexpected strand-specifc counts with Ensembl GTF
Kaur Alasoo [guest]
guest at bioconductor.org
Wed Aug 27 16:53:58 CEST 2014
Hi,
I am trying to use featureCounts to perform strand-specifc RNA-Seq read counting with an Ensembl GTF file. However, it seems to count reads from each strand separately rather than for each gene from the strand where the gene is located. I have created a small BAM that contains reads from only two genes (one from each strand) to illustrate the problem.
First, if I specify strandSpefic = 0 I get 97.7% of the reads aligned to two genes:
a = featureCounts("merged.bam", annot.ext = "Homo_sapiens.GRCh38.76.gtf",
strandSpecific = 0, isGTFAnnotationFile = TRUE, isPairedEnd = TRUE)
> a$stat
Status counts.merged.bam
1 Assigned 16179
2 Unassigned_Ambiguity 0
3 Unassigned_MultiMapping 56
4 Unassigned_NoFeatures 335
And the reads counts for the two genes are:
> a$counts[c("ENSG00000066336","ENSG00000136869"),]
ENSG00000066336 ENSG00000136869
11759 4419
However, when I try to specifiy strandSpecific = 1 or 2 I get reads from only one strand counted in either case
b = featureCounts("counts/merged.bam", annot.ext = "../../annotations/GRCh38/genes/Homo_sapiens.GRCh38.76.gtf", strandSpecific = 1, isGTFAnnotationFile = TRUE, isPairedEnd = TRUE)
> b$stat
Status counts.merged.bam
1 Assigned 5295
2 Unassigned_Ambiguity 0
3 Unassigned_MultiMapping 56
4 Unassigned_NoFeatures 11219
> b$counts[c("ENSG00000066336","ENSG00000136869"),]
ENSG00000066336 ENSG00000136869
952 4343
ENSG00000136869 is located on the forward strand, gets high counts with strandSpecific = 1.
#####
c = featureCounts("counts/merged.bam", annot.ext = "../../annotations/GRCh38/genes/Homo_sapiens.GRCh38.76.gtf", strandSpecific = 2, isGTFAnnotationFile = TRUE, isPairedEnd = TRUE)
> c$stat
Status counts.merged.bam
1 Assigned 10884
2 Unassigned_Ambiguity 0
3 Unassigned_MultiMapping 56
4 Unassigned_NoFeatures 5630
> c$counts[c("ENSG00000066336","ENSG00000136869"),]
ENSG00000066336 ENSG00000136869
10807 76
ENSG00000066336 is located on the reverse strand, gets high counts with strandSpecific = 2.
Furthermore,
a$counts[c("ENSG00000066336","ENSG00000136869"),] = c$counts[c("ENSG00000066336","ENSG00000136869"),] + b$counts[c("ENSG00000066336","ENSG00000136869"),]
I downloaded the Ensembl 76 GTF file from here:
ftp://ftp.ensembl.org/pub/release-76/gtf/homo_sapiens
The sample BAM files is available from here:
https://dl.dropboxusercontent.com/u/6006832/merged.bam
Is this a bug of featureCounts or am I doing something obviously wrong?
Kind regards,
Kaur Alasoo
PhD student at Sanger Institute
-- output of sessionInfo():
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] Rsubread_1.12.6 BiocInstaller_1.12.1
loaded via a namespace (and not attached):
[1] tools_3.0.2
--
Sent via the guest posting facility at bioconductor.org.
More information about the Bioconductor
mailing list