View Single Post
Old 09-24-2019, 02:51 AM   #2
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 622
Default

Dear Ali,

Thanks for your kind words, and for your thoughtful questions. As you will see below, I will probably not be able to give you a satisfactory answer to all questions you raised, but I will try to share my view on some of the issues nevertheless.

RRBS data has always looked quite ‘funky’ when it comes to M-bias plots. We have so far mostly chosen to simply accept this ‘as is’, especially given that we have more or less not used RRBS ourselves for more than six years…

To Question 1:
I don’t think the overall bisulfite conversion efficiency should necessarily be judged based on overall methylation percentage. The report you attached shows that the non-CG methylation is ~2.8% overall (which would mean a conversion efficiency of at least 97.2%), but one can see that the M-bias plots are not at all behaving uniformly. It rather looks like the overall non-CG methylation is well under 1% for most positions (just mouse-over in the plot, they are probably ~0.4-0.6% mostly), but there are some positions that show around 16-20% methylation. Such positions (see also Q2), will have a big impact on the average methylation percentage, and therefore get an unfair say in judging the conversion efficiency. We would argue that conversion efficiency should not discriminate by position (or context), so the lowest methylation average you see anywhere in the reads or in the genome can be used as a proxy for conversion efficiency.

This means that if you see a non-CG methylation of ~0.4% for most parts of all reads, that number has to be the combination of i) true non-CG methylation, ii) bisulfite conversion failure and iii) mismapping effects. If you now assume that there are hardly that many mismapping effects, and that there is hardly any non-CG methylation in the cell type you are looking at, then it would mean that virtually all of the 0.4% methylation are conversion errors (in reality it is probably a combination of all three effects though). So in the worst case, I would argue that the conversion efficiency must have been 99.6% efficient, or a bit more even. A value that I would find perfectly acceptable.

To Questions 2 and 3:
Spikes at individual positions: I would think that such positions come from very repetitive regions in the genome, and are possibly the result of mis-mapping or suffer from conversion failure because of some kind of higher order structure (something demonstrated very convincingly by your colleagues for methylation of the mitochondrium, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5671948/).

Spikes at the 3’ end of Read 2 are the result of stringent adapter trimming. The Illumina adapter starts with AGATC…, so reads will never end in A, AG, AGA, etc. Read 1 will participate in methylation calling only at C or T positions, however since Read 2 is the reverse complement of R1, methylation calling occurs at G and A positions. Since reads may never end in A, this also means that the very last position of a Read 2 may never be found in an unmethylated state. While it is true that this in theory introduces a bias for that very last position, one should take into consideration that:

a) it really only ever occurs for at the very last position of R2 which is not always the same for every read (arguably more likely for RRBS thought),

b) the total number of calls at the very last position is typically quite low (just mouse over for details)

c) the very last position is of R2 is also subject to overlap removal if the read overlaps with R1 (fairly likely).

In other words: Yes, this position is biased towards being called methylated, but it will almost certainly not have any impact on your results as a whole whatsoever.

Regarding the spiked positions again: You should be able to look at the genomic distribution of alignments. I would predict that there will be certain positions in the genome (e.g. close to the edges of chromosomes or centromeres, the MT etc). where there you will find thousands of reads aligned to the very same position (which could harbour the conversion artefacts). Depending on how you move on with downstream analysis, these positions might be completely irrelevant for your further results. While these positions can have quite some influence on the overall numbers and average stats (the ones you find in the Bismark report), but if you would call the average methylation over larger regions you could collapse the methylation values of tens of thousands of reads down to a single methylation percentage. In such an analysis, the artefactually high read coverage would have no higher say than any other region the genome.

To conclude, I would not hesitate to continue working with the data. And in any case, once you found potential regions of interest you should go back to the original data and convince yourself that you trust the underlying signal at that position.

I hope this helps a little.
All the best, Felix
fkrueger is offline   Reply With Quote