This is a test version of Biostars. For the public version, visit https://www.biostars.org.
featureCounts() in Rsubread freezes on gene_biotype as the meta-feature, command-line version doesn't

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

0 answers

No answers yet.

Log in to answer this question.