I am currently writing an R package that includes a module to run featureCounts (gene quantification tool) from Rsubread. I wanted to be able to specify the correct strandedness option to featureCounts without the user needing to specify whether their protocol is unstranded, forwardly stranded, or reversely stranded. Typically when I don't know the protocol for a given BAM file, I use RSeQC's infer_experiment module, which uses the BAM plus an annotation (GTF) file to give an output like so:
Unstranded:
This is PairEnd Data
Fraction of reads failed to determine: 0.0172
Fraction of reads explained by "1++,1--,2+-,2-+": 0.4903
Fraction of reads explained by "1+-,1-+,2++,2--": 0.4925
Reversely stranded (R2 same strand as parent gene):
This is PairEnd Data
Fraction of reads failed to determine: 0.0072
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0487
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9441
Does anyone know of an R package that does something similar?
If I really needed to I could always require that users have RSeQC installed and just run infer_experiment as a system call. But it would really be better to run it from an existing R package if available.
I know I could also write something up myself (use rsamtools to downsample the BAM, check ratio of SAM flags 99/147 vs. 83/163 in regions matching positive and negative strand genes). But would rather not reinvent the wheel here if there is something already available to do what I need.
