featureCounts() in Rsubread freezes on gene_biotype as the meta-feature, command-line version doesn't
0
1
Entering edit mode
4.6 years ago
jirkanov ▴ 10

Hello,

I want to use gene_biotype as the meta-feature. However, featureCounts() from Rsubread (version 1.34.7) freezes when I do this:

res <- featureCounts(
  files = "gsnap_out/SRR1039508.bam",
  annot.ext = "../Homo_sapiens.GRCh37.75.gtf",
  isGTFAnnotationFile = TRUE,
  GTF.featureType = "exon",
  GTF.attrType = "gene_biotype",
  isPairedEnd = TRUE,
  countMultiMappingReads = FALSE,
  allowMultiOverlap = FALSE,
  countChimericFragments = FALSE,
  nthreads = 4
)

I kept it running for several hours, then I had to kill it.

When I try to do the same in command-line featureCounts (version 2.0.0), it works:

$ featureCounts -a ../Homo_sapiens.GRCh37.75.gtf -t exon -g gene_biotype -o counts.tsv -p -T 4 gsnap_out/SRR1039508.bam
$ cut -f 1,7 counts.tsv
# Program:featureCounts v2.0.0; Command:"featureCounts" "-a" "../Homo_sapiens.GRCh37.75.gtf" "-g" "gene_biotype" "-o" "counts.tsv" "-p" "-T" "4" "gsnap_out/SRR1039508.bam" 
Geneid  gsnap_out/SRR1039508.bam
pseudogene  54632
lincRNA 6359
protein_coding  109356
antisense   4546
processed_transcript    818
snRNA   110
sense_intronic  107
miRNA   157
misc_RNA    355
snoRNA  81
rRNA    24
3prime_overlapping_ncrna    0
polymorphic_pseudogene  0
sense_overlapping   24
IG_V_gene   0
IG_C_gene   0
IG_J_gene   0
IG_V_pseudogene 0
TR_C_gene   0
TR_J_gene   0
TR_V_gene   0
TR_V_pseudogene 0
IG_C_pseudogene 0
TR_D_gene   0
TR_J_pseudogene 0
IG_J_pseudogene 0
IG_D_gene   0
processed_pseudogene    0
Mt_tRNA 0
Mt_rRNA 0

Note that using gene_id as the meta-feature works normally:

res <- featureCounts(
  files = "gsnap_out/SRR1039508.bam",
  annot.ext = "../Homo_sapiens.GRCh37.75.gtf",
  isGTFAnnotationFile = TRUE,
  GTF.featureType = "exon",
  GTF.attrType = "gene_id",
  isPairedEnd = TRUE,
  countMultiMappingReads = FALSE,
  allowMultiOverlap = FALSE,
  countChimericFragments = FALSE,
  nthreads = 4
)

res$stat

                          Status SRR1039508.bam
1                       Assigned         124787
2            Unassigned_Unmapped          17518
3           Unassigned_Read_Type              0
4           Unassigned_Singleton              0
5      Unassigned_MappingQuality              0
6             Unassigned_Chimera         631415
7      Unassigned_FragmentLength              0
8           Unassigned_Duplicate              0
9        Unassigned_MultiMapping         384231
10          Unassigned_Secondary              0
11           Unassigned_NonSplit              0
12         Unassigned_NoFeatures         658303
13 Unassigned_Overlapping_Length              0
14          Unassigned_Ambiguity           6665

Am I doing something wrong or is it some bug? Anyway thanks in advance for help!

subread featurecounts rsubread RNA-Seq • 1.3k views
ADD COMMENT

Login before adding your answer.

Traffic: 3298 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6