Comparing two RNAseq experiments using summarizeOverlaps
Entering edit mode
Last seen 4.7 years ago


I have two RNA-seq experiments. I first checked them to see whether they are single or pairedEnd using:

bamtools view -c -F 1 myfile.bam, bamtools view -c -f 1 myfile2.bam

There I saw that all reads are singleEnd. In the next step I uploaded my data in R, and counted the reads per feature by summarizeOverlaps. 

x <- sumarizeOverlaps(features= myfeatures, reads= myrnaseq1, mode= IntersectionNotEmpty, singleEnd=T)

I repeated this with the second experiment. In the next step I performed a quality control to see how similar both experiments are: Basically by a dotplot in ggplot and adding a linear model (geom_smooth(method="lm")).. Surprisingly, the line was horizontal. However, after adding the argument ignore.strand=T, it works perfectly. 

Do you have an explanation for that? Might there be problems with the protocol used to generate the data? 



summarizeoverlaps R bioconductor rnaseq rna-seq • 1.2k views
Entering edit mode
Last seen 10 hours ago
United States

Many RNA seq protocols are not sensitive to strand, so setting ignore.strand=TRUE sounds reasonable. It's a little hard to be sure, because you haven't described what myfeatures and myrnaseq1 are.

Entering edit mode
Last seen 5.5 years ago
Australia/Centenary Institute Universit…

The rseqc tool has a very useful tool called which you can run on your bam file (you'll need a reference annotation) to see which stranding protocol your data looks like it came from. This is very useful post-hoc even when you know the protocol of how your library was generated to check that "stranding" worked, but in your case it can give you the confidence that setting ignore.strand=TRUE is an actually valid thing to do (i.e. justifiable) not just that it works so that's why you did it. 

Entering edit mode
Last seen 21 hours ago
Seattle, WA, United States

FWIW you could also try the strandedTest() function defined in the examples of ?junctions in the GenomicAlignments package to check whether the RNA-seq protocol was stranded or not.



Login before adding your answer.

Traffic: 264 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6