Hello,
I have RNAseq data where I want to generate bigwigs with just the 3' end position of the data, which corresponds with the first base of my read 1 in a reverse-stranded library.
My initial bamCoverage command to generate these bigwigs was:
bamCoverage -b ${sample} -o $bigwig_directory/${base}_CPM_bs10_forw2.bw -of bigwig --filterRNAstrand forward --normalizeUsing CPM --Offset -1 --binSize 10
bamCoverage -b ${sample} -o $bigwig_directory/${base}_CPM_bs10_rev2.bw -of bigwig --filterRNAstrand reverse --normalizeUsing CPM --Offset -1 --binSize 10
Two questions have come out of this:
1) I believe I am mistaken and should have actually set offset to 1 to get the first basepair of read 1, is that correct?
2) More confusingly, it appears that this command is resulting in peaks in my bigwig file corresponding to both the end of read 1 as well as the end of read 2.
For example, here I have read 1s ending in one bin, and read 2s ending in a separate, adjacent bin. I'd expect only one of the two bins to have signal when I convert to bigwig using the command above.
However, the corresponding bigwig file has signal both from 700-710 and from 710-720, whereas I'd expect only 1 of those 2 bins of 10 to have signal if my binsize=10 and my reads were only being mapped to the end of one strand? I have seen other examples of this for many genes, always the same pattern, and not only for cases where read ends are adjacent to each other.
Perhaps I am missing something obvious, but could someone please explain where the discrepancy is coming from? Devon Ryan
Thank you for your reply!
Does filterRNAstrand not remove read 2?
From the documentation for filterRNAstrand: "Selects RNA-seq reads (single-end or paired-end) originating from genes on the given strand. This option assumes a standard dUTP-based library preparation (that is, –filterRNAstrand=forward keeps minus-strand reads, which originally came from genes on the forward strand using a dUTP-based method)."
My understanding was that "–filterRNAstrand=forward keeps minus-strand reads" means that only read 1 is kept for a + strand gene (where read 1 is the minus-strand read for my reverse stranded library). Perhaps I misunderstood?
Yes, the mix up is maybe that your RNA can come from plus or minus strand.
Depending on which strand the RNA originated from, Read1 will correspond the 3'end of that fragment. So read 1 can map to the plus or minus strand, it just depends on which strand the RNA came from.
Accordingly, read2 is just the other side of your fragment, so it can also map to the plus or minus strand, depending on the original RNA fragment.
The ability to identify strand comes from a known fact on whether Read1 is sense to antisense to your RNA (forward or reverse stranded library). This is determined by specific steps in the library prep. So if your read1 is sense with the original RNA, then if the original RNA comes from the plus strand, read1 will map to the plus strand and read 2 will map to the minus strand.
If your RNA came from the minus strand, then Read1 maps to the minus strand and Read2 maps to the plus strand.
If your library is reverse stranded, then you expect the opposite, where Read1 maps to the opposite strand from the original RNA, and so forth.
So filterRNAstrand operates in the assumption of a reverse stranded library. And wants to give you information about the original RNA fragment.
So if the original RNA fragment came from the plus (forward) strand, then Read1 will map to the minus strand, BUT read2 will map to the plus strand. It will count both of these reads as a fragment on the plus strand.
I know it's easy to get turned around with these, so let me know if I'm not clear.
In short, if you are positive that Read1 is what you want, you can use samtools to filter on the sam flags for first in pair reads, then make bigwigs.
Oh I see! Yes, this is very clear, thank you so much for the in-depth answer :)
I think what I was missing was exactly what you said - that read2 is also still going to be mapped, and that when the documentation says "–filterRNAstrand=forward keeps minus-strand reads," it means those reads are also kept, not that they are exclusively kept.
I will play around with the samflag exclude/include parameters in bamCoverage in this case.
Thank you very much again!