Comparing two RNAseq experiments using summarizeOverlaps
3
0
Entering edit mode
@walter-f-baumann-12439
Last seen 3.8 years ago

Hi,

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 • 957 views
1
Entering edit mode
@martin-morgan-1513
Last seen 9 weeks 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.

1
Entering edit mode
@darya-vanichkina-6050
Last seen 4.6 years ago
Australia/Centenary Institute Universit…

The rseqc tool has a very useful tool called infer_experiment.py 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.

1
Entering edit mode
@herve-pages-1542
Last seen 1 day 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.

H.