I am trying to identify enriched regions(peaks) between m6A RNA methylation RIP-seq and input data. For that, the human genome need to be split into non-overlapping 20 nt windows, then I need to calculate the coverage for each window. Finally, I will use the Fisher test to test the window enrichment between IP and input.
The available files in GEO ( https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE83561 ) are bigwig files ( signal values represent normalized reads per million of mapped reads) for each IP and input replicates.
Can I convert the bigwig to bedGraph file and use these coordinates for calculating the coverage, or I need the original Sam to get the coordinates of all the mapped reads?
Bioinformatics.
In theory, it should be ok. But I wouldn't recommend it. There are big differences in how the signal could have been normalized - is the sum of all the reads calculated from all the mappings (=each read could be counted more than once), all the multi-mapped reads or unique mappings only? Also, even though it should be the same in different tools (I mean the CPM normalization) I have seen huge differences between STAR and for example DeepTools tracks.
Sometimes, the conversion bedGraph->bigWig "rounds" up the signals so you don't have the exact counts as you would have in the original bedGraph. This is because the main purpose of bigWig is visualization and for that the "rounding" is fine. I would go for the original SAM/BAM if I were you but if you want to have a quick check the bigWig should work as well.
Hello Jan Oppelt. Thanks for your answer. Indeed, I converted the bigwig into bedgraph and use it (the authors claim they selected only uniq reads), but the results were not why I expected. I followed the authors protocols, so, when I performed the fisher test for each genome window (IP vs. input), I used a matrix as follows:
IP Input
bin IP read counts bin Input read counts
Total IP mapped read Total Input mapped reads
Do you think the matrix could have generated some problem? Should I change its orientation? As you see, the matrix contain read counts, that's why I think I should use the original bam file with all mapped reads.
Something I noted is that the fisher odds ratio was several time "Inf", because there are several windows without mapped reads.
An enriched windows must have a p-adjusted
Well, if the table needs to have Read count or Total mapped reads you cannot use the normalized bigWig because it doesn't preserve the information about the actual number of reads (as you've mentioned, it's normalized as cpm).
Have you tried to use some peak caller (usually used for ChIP-seq analysis) on your bedgraph files? Some accept bedgraph as input, such as MACS / MACS2 and I'm pretty sure you can specify a window or a "fragment size" to look at, if you want to keep you 20 bp window. Good luck!
I found the Macs function "bdgpeakcall " which calls peaks from bedGraph output, but it only works with one input file. Instead I need to identify enriched peak in IP sample compared with input. Another detail is this command works with p-value or q-value scored bedgraphs files; the files I have contain only RPM coverage signal.
My mistake, I was sure you could perform a peak calling (macs2 callpeak) could manage bedgraphs input but it' not the case. For other macs2 function using bedgraphs (diffpeak, bdgpeakcall), they require bedgraphs from macs output, which have p- or q-value scores as you said... One solution might be to use bigwigCompare from DeepTools which outputs one bedgraph file of a bw file subtracted from another and then run macs2 bdgpeakcall. It may be more easy to ask for the bam files, or even to download the fastq files in SRA and restart from alignment ! Good luck!
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics