I have some relatively recent strand specific RNA-seq data which comes as 'second-strand', that means all reads are reversed with respect to the genomic sequence of the transcript and annotation. With featureCounts I would have to use -s 2 option, for example, to get a 70+% rate of assignments to transcripts. But now I am very concerned that the counting went terribly wrong. Does easyRNAseq have an option to set the strand for read assignment, or does it determine those automatically?

I have counted those with the following settings:

bamParam <- BamParam(paired = FALSE, stranded = TRUE, yieldSize = 1000000L)
param <- RnaSeqParam(annotParam = annotParam, bamParam = bamParam, countBy = c("transcripts"), precision = "read")
ret <- simpleRNASeq(bamFiles = getBamFileList(args), param = param, nnodes = 120, verbose = TRUE, override = TRUE)

Hi Michael,

FWIW, you could double-confirm your results by counting again with summarizeOverlaps() from the GenomicAlignments package. If your data is "reversely stranded", then invert the read strand with invertStrand(reads) before passing them to summarizeOverlaps(). Alternatively, if the reads are paired-end, your can specify strandMode=2 when loading them with readGAlignmentPairs().

It might not be as efficient as featureCounts() so maybe you want to try this on a subset of your data only, just to confirm the "reversely stranded" thing.

Also there is a simple test in ?junctions (in GenomicAlignments) to find out whether RNAseq data is stranded or not. It doesn't use any annotation at all, it only looks at the read junctions.

Cheers,

H.

I have done a comparison using the featureCounts program, seemingly the worst-case scenario happened. The number of assigned reads in featureCounts using -s 2 is on average 10x-15x higher per sample. So it seems that all our counts on recent stranded Illumina data could be wrong! This is possibly a big bummer, everyone should check their counts if they were achieved using the following parameters:

• recent Illumina strand-specific protocol (reverse sense output)
• counted using easyRNAseq with option stranded = TRUE