This is actually a follow-up from a previous problem I had described in this post.
In this experiment we've taken cultured human cells, then done the NEBNext ultra directional library prep using their ribosomal RNA depletion module, before sequencing 2x43bp PE on a NextSeq and mapping with TopHat.
In that earlier post I struggled to get SeqMonk to show what I thought were reasonable values in the RNA QC report plot, which I thought I had fixed by using the default Ensembl genome annotation (instead of a custom genome based on hg19, which I had used for my mapping). It gave low values of rRNA reads (which is what I would hope for, given they were ostensibly depleted), and other seemingly reasonable (and crucially, non > 100%) values. I then later re-mapped the data using the Ensembl distribution from iGenomes, re-ran SeqMonk and got the same thing (attached pic).
However, parallel to that I had also been quantitating reads with featureCounts (using the appropriate whole-genome gtf annotation file), and seemed to be getting fewer feature counts than expected, given how many reads were mapped. To double-check the rRNA content, I then downloaded an rRNA/tRNA specific gtf from UCSC and re-ran featurecounts with that. Summing the total number of features counted using the rRNA gtf output about five times the number of features detected when using the whole genome gtf (about 50% and 10% of the total fragments detected respectively).
So it seems to me like these two strategies are given rather contradictory accounts: SeqMonk is directly saying very little of our mapped data is rRNA, while featureCounts is indirectly saying that most of it is (which is entirely possible, given that I do not have a positive control to show whether the rRNA depletion worked or not). I would like to use the non-overlapping features of both programs, but I can't really justify that if I can't get them to agree in the areas where their output does overlap.
I assume that I am making a silly technical or conceptual mistake, but I've been over it a few times now and I can't see what that would be. I assume I have done something wrong with the data import or quantitation in SeqMonk, as it's the area I'm less familiar with, but I can't think what. As always, I'd be very grateful for any input!
PS on a related note to this I would also love to hear from people using the same kit what proportions of reads you might expect to typically a) map at all, and b) map to rRNA/fail to map to coding sequences - this is my first exploration of rRNA-depleted data, and I don't yet have a feel for how well this protocol should perform in the real world.
In this experiment we've taken cultured human cells, then done the NEBNext ultra directional library prep using their ribosomal RNA depletion module, before sequencing 2x43bp PE on a NextSeq and mapping with TopHat.
In that earlier post I struggled to get SeqMonk to show what I thought were reasonable values in the RNA QC report plot, which I thought I had fixed by using the default Ensembl genome annotation (instead of a custom genome based on hg19, which I had used for my mapping). It gave low values of rRNA reads (which is what I would hope for, given they were ostensibly depleted), and other seemingly reasonable (and crucially, non > 100%) values. I then later re-mapped the data using the Ensembl distribution from iGenomes, re-ran SeqMonk and got the same thing (attached pic).
However, parallel to that I had also been quantitating reads with featureCounts (using the appropriate whole-genome gtf annotation file), and seemed to be getting fewer feature counts than expected, given how many reads were mapped. To double-check the rRNA content, I then downloaded an rRNA/tRNA specific gtf from UCSC and re-ran featurecounts with that. Summing the total number of features counted using the rRNA gtf output about five times the number of features detected when using the whole genome gtf (about 50% and 10% of the total fragments detected respectively).
So it seems to me like these two strategies are given rather contradictory accounts: SeqMonk is directly saying very little of our mapped data is rRNA, while featureCounts is indirectly saying that most of it is (which is entirely possible, given that I do not have a positive control to show whether the rRNA depletion worked or not). I would like to use the non-overlapping features of both programs, but I can't really justify that if I can't get them to agree in the areas where their output does overlap.
I assume that I am making a silly technical or conceptual mistake, but I've been over it a few times now and I can't see what that would be. I assume I have done something wrong with the data import or quantitation in SeqMonk, as it's the area I'm less familiar with, but I can't think what. As always, I'd be very grateful for any input!
PS on a related note to this I would also love to hear from people using the same kit what proportions of reads you might expect to typically a) map at all, and b) map to rRNA/fail to map to coding sequences - this is my first exploration of rRNA-depleted data, and I don't yet have a feel for how well this protocol should perform in the real world.