Bismark - A New Tool for Mapping and Analysis of Bisulfite-Seq Data

[pig_raffles]

Hi Felix,

I managed to get the script to work by having all the files (input, output) in the same directory as the bismark2Bedgraph script.

I could not find anything wrong with the file locations/file names that I had used previously but it works now!

Thanks again for your time and assistance

[daisyko]

Originally Posted by fkrueger

Hi seqfast,

When you are generating and sequencing PCR amplified regions it is quite common that you see only the top strand (with both the OT and CTOT strands) or the bottom strand (with both OB and CTOB as in your case), depending on which strand you targeted when designing the primers.

To help getting your head around this it really helps to draw the sequence of an amplicon of interest out on a sheet of paper. You should see that the bisulfite treatment will change one of the two strands so much that the oligos will only amplify one of the two strands, and this PCR product is then usually sequenced from both sides, e.g. the CTOB and OB strands, but this may be different depending on how the library preparation was done.

So in short: PCR amplicons are not normally directional libraries but you only sequence both versions of either the top or the bottom strand, depending on how the primers were designed. I hope this helps, let me know if I can be of any further assistance. Cheers, Felix

Hi Felix or any other experts,

I encountered similar problem which most of the reads from a paired-end library aligned to the bottom strand after running Bismark (with either directional option or with the non-directional option). I would like to know that how we could process the further analysis and how we would calculate the coverage for each site if the reads from bottom strands and from top strands are extremely unbalanced.
Sorry for the inconvenience caused and hope to get a response from any of you guys.

Best regards,
Daisy

[fkrueger]

Hi Daisy,

Were you expecting to see imbalances in your library, i.e. was it some kind of target enrichment or (PCR) amplification library? Depending on what you expect from the library design you should only analyse alignments from the strands you are expecting, so if you only expect OT or OB alignments you should not perform non-directional alignments.

If your library strategy was designed in a way that you only amplified or pulled down one the two strands prior to preparing the libraries you will indeed only see the methylation of either the top or the bottom strands, but not both. The coverage in that case would simply be the number of reads you see for a region (unless you used UMIs, in which case you could also take uniqueness of fragments into account...).

[daisyko]

Originally Posted by fkrueger

Hi Daisy,

Were you expecting to see imbalances in your library, i.e. was it some kind of target enrichment or (PCR) amplification library? Depending on what you expect from the library design you should only analyse alignments from the strands you are expecting, so if you only expect OT or OB alignments you should not perform non-directional alignments.

If your library strategy was designed in a way that you only amplified or pulled down one the two strands prior to preparing the libraries you will indeed only see the methylation of either the top or the bottom strands, but not both. The coverage in that case would simply be the number of reads you see for a region (unless you used UMIs, in which case you could also take uniqueness of fragments into account...).

Dear Felix,

Thanks a lot for answering.

So I think in that case my enrichment sequencing is targeted at one strand, I should only look into that strand. That really solves my problem

[shawpa]

I am not sure if this is the proper thread to post this on but since I am using bismark output, I thought I would ask here. I have very low coverage WGBS data. Ideally I would have higher coverage and/or more samples per group but right now, I only have 2 samples I am comparing with low coverage. I was trying to see if I could identify any DMR using a sliding window analysis. After I did that calculation, I wanted to see what the individual sites in those regions looked like in terms of per CpG methylation. What I think I am seeing in many cases is the same CpG site but the forward and reverse strand so it appears 1 base apart. However, they are showing pretty different methylation values at what should be the same CpG site. It also shows really different coverage values too. I assume this is some kind of mapping bias but I am not sure. Is what I am seeing "normal" for this type of data or should it be more consistent at the sample CpG site? I am having trouble trusting my window analysis when I see the methylation values jump around so wildly. Attached is a text file with all the sites in 1 region of my 2 samples.

[fkrueger]

Just generally, in order to get reads for both the forward and reverse strands you need to have a fairly good sequencing depth and have a library with a good complexity. While your sample 1 has several CpG dinucleotides covered on both strands, Sample 2 has hardly any covered on more than one strand.

By just glancing at over the values in your example the values don't seem to be differing all that wildly to be honest. If you for example take the example:

224012575	224012575	10	50	50
224012576	224012576	3	66.66666667	33.33333333

Here the 2 positions don't agree perfectly but with 3 reads in total there simply are only a limited number of percentages you can achieve, namely 0, 33, 66 or 100%. There will always be a problem to match numbers perfectly if you have a shallow read depth, but this is why you need to average the values over larger distances to even those limitations out.

[akramdi]

Hi Felix,

I am trying to make sens of my M-bias plot.
I have 150nt PE sequencing and the M-bias plot shows a gradual decrease in CHH calls for R2 starting from ~75nt.

R2 are trimmed and per base quality plot looks fine. What concerns me is that I could not find the same CHH calls when I looked into the methylation call strings (XM:Z). It works for R1; I am able to reproduce the exact same Mbias plot, but for R2, the methylation call strings do not seem to reflect the Mbias plot.

Moreover, when processing R2 as SE, the Mbias plot shows a fairly even distribution of CHH calls across the length of the reads (no decrease).

Could you help me understand this gradual decrease in PE mode? why does it disappear in SE mode?

Thanks for the help!
Cheers,

Amira

[fkrueger]

Hi Amira,

The reason for this is the overlap detection, and -removal, in paired-end files. What this simply does is to detect when R1 and R2 start overlapping, and as soon as this is the case the methylation extraction from R2 is halted. This mechanisms avoids introducing a coverage bias for parts of the reads that are overlapping. If you would specify the option --include_overlap you will not see this dramatic decrease anymore, but this is not to say that you should be using it. Does this help? Cheers, Felix

[akramdi]

Originally Posted by fkrueger

Hi Amira,

The reason for this is the overlap detection, and -removal, in paired-end files. What this simply does is to detect when R1 and R2 start overlapping, and as soon as this is the case the methylation extraction from R2 is halted. This mechanisms avoids introducing a coverage bias for parts of the reads that are overlapping. If you would specify the option --include_overlap you will not see this dramatic decrease anymore, but this is not to say that you should be using it. Does this help? Cheers, Felix

Hi Felix,

I wasn't aware that the overlapping region of R2 was removed for the M-bias plot too, I thought that R1 and R2 were processed independently.

I understand my plots now thanks!

Cheers,
Amira