3
$\begingroup$

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.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.