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

[fkrueger]

Originally Posted by gene_x

You mean "This should affect all C contexts and not only CHH?
How do you click the CHH and CHG away? I can only zoom in certain bases...

If you click on e.g. 'CHH total calls' in the legend that track will disappear/show again and the scale adjust accordingly.

Can you help me clarify one thing? In the M-bias plot for read2, is position 1 the start of the read2? I'm a bit confused about the orientation of the read2. I thought in M-bias plot, position 1 is the the end of read2 (after read 2 is reverse complemented)? probably my understanding was wrong...

The first position does indeed correspond to the first base or read 2 as it is read by the sequencer. If this was the reverse complemented end of the sequenced fragment most of the artefacts would be completely obscured because these would likely happen at various positions within the read.

Also, in M-bias plot for read2, the first 2 bases have much lower average CpG methylation, that suggests they should be ignored in the downstream analysis? I could use the -ignore_r2 option in methylation extractor to do it?
Thanks!

The first positions of read 2 are the result of carrying out the end-repair/A-tailing procedure with unmethylated cytosines. Sonication is quite likely to leave lagging ends, so these would be filled in with unmethylated Cs (of course only at cytosine positions) which would subsequently be converted to Ts in the bisulfite reaction. If you then sequence Read 2, which is the reverse complement of such a T, it would appear as A and be called as unmethylated. As you rightly said this can be remedied by using --ignore_r2 2.

[christophe]

Originally Posted by adeirossi

Hi Christophe,

I have also seen a falling number of methylation calls for read 2 in my data. I have always assumed this was because I ran bismark_methylation_extractor with the --no_overlap option and many of the read pairs had overlaps (i.e. insert length less than twice the read length). With --no_overlap, when the extractor encounters an overlap between read 1 and read 2, I believe it uses the read 1 calls only and discards the read 2 calls, which would explain why the number of read 2 calls decreases as you go toward its 3' end. Because this "no overlap" filtering occurs after mapping, I assume there is no impact on the alignment of read 2.

If, on the other hand, you are not using the --no_overlap option or you don't have many reads pairs with inserts <200 bp, then I have no idea what's going on.

Andrew

Hi Andrew and Felix,

Sorry for this late message.

Thanks Andrew for your reply to my question regarding M-bias.
I effectively ran the Extractor with the --no_overlap option on before. As suggested I just rerun Bismark Methyl Extractor without the --no_overlap option and it worked. I do not see the number of calls falling anymore for CHH and CHG.

Thanks Felix for confirming Andrew's answer. Yes, you are right, the CpG calls was also falling slightly but better now.

I also noticed that the R2 calls counts are more heterogeneous across the read length than the number of calls across R1. I do not know if anyone else noticed that.

Thanks,
Best,
Christophe

[gene_x]

Hi, Felix,
I have another question about the alignment files. I just checked the bam file after running bismark alignment and it is unsorted. Then using the bam output to do deduplication, the resulting deduplicated.bam is also unsorted. My first question is deduplication step doesn't require the input bam to be sorted?

I was trying to merge several deduplicated.bam files from different lanes into one big final bam and run methylation extractor on it and the resulting coverage file (.zero.cov with --zero_based specified) has methylation calls on non-Cs locations like chr1 133 134. That's how I traced it back to the issue of unsorted bam..

So it's probably better to implement sorting into the bismark pipeline?

[dpryan]

Originally Posted by AndrewMartins

Your post is about a tool that can be used for to would be completely obscured because these would likely happen at various positions within the read. i can told you for use another site of tool its mainly supply pipeline on customers.

Gotta love incoherent spam

[fkrueger]

Originally Posted by gene_x

Hi, Felix,
I have another question about the alignment files. I just checked the bam file after running bismark alignment and it is unsorted. Then using the bam output to do deduplication, the resulting deduplicated.bam is also unsorted. My first question is deduplication step doesn't require the input bam to be sorted?

I was trying to merge several deduplicated.bam files from different lanes into one big final bam and run methylation extractor on it and the resulting coverage file (.zero.cov with --zero_based specified) has methylation calls on non-Cs locations like chr1 133 134. That's how I traced it back to the issue of unsorted bam..

So it's probably better to implement sorting into the bismark pipeline?

It should be generally possible to merge several bam files before running the extractor without a problem, and files do not have to be sorted for this. You only need to ensure that read 1 and read 2 always follow each other in the bam file, else you would run into problems with stand identity and the no_overlap function. Some merging tools do not guarantee this order, in which case it might help to sort the reads by read id to team pairs up again. I'm currently at a conference so I'm somewhat slow replying.

[gene_x]

Originally Posted by fkrueger

It should be generally possible to merge several bam files before running the extractor without a problem, and files do not have to be sorted for this. You only need to ensure that read 1 and read 2 always follow each other in the bam file, else you would run into problems with stand identity and the no_overlap function. Some merging tools do not guarantee this order, in which case it might help to sort the reads by read id to team pairs up again. I'm currently at a conference so I'm somewhat slow replying.

I indeed run into problem after sorting bam files based on chromosomal coordinates.. the methylation extractor report an error saying that read1 and read2 are not the same and suggested me to use an unsorted bam..

I then tried to sort the deduplicated bam with the command

samtools sort -n deduplicated.bam deduplicated_sort

And then run methylation extractor on the deduplicated_sort.bam, however, the same error message came up again and methylation extractor couldn't proceed.. Is it because after deduplication, read1 and read2 within one bam file don't always match even after I sort them based on read names? Should I not sort anything from the beginning at all?

Thanks!

[fkrueger]

Indeed, you should not sort the files by coordinates at all before running the deduplication.

[gene_x]

Originally Posted by fkrueger

Indeed, you should not sort the files by coordinates at all before running the deduplication.

Should I sort it by read names using -n in samtools sort just like I listed before doing deduplication? Or should I not sort it at all?

[fkrueger]

Of you can, just use the Bismark files as they are generated. If you have to
using samtools sort -n should do the trick as well.

[gene_x]

Originally Posted by fkrueger

Of you can, just use the Bismark files as they are generated. If you have to
using samtools sort -n should do the trick as well.

OK. I'm running that now.

Another question, I'm not sure if my previous run without doing any sorting was correct.

Here is my command

samtools merge input.bam plate1/plate1_all_sort.bam plate2/plate2_all_sort.bam
deduplicate_bismark -p --bam input.bam
bismark_methylation_extractor -p --no_overlap --ignore 3 --ignore_r2 3 --bedGraph --counts --zero_based --report input.deduplicated.bam 2> input.meth_extractor_log.txt

The resulting input.deduplicated.zero.cov file have methylation coverage on bases that are not Cs. Here are the top rows:

chr2 133 134 0 0 8
chr2 134 135 0 0 1
chr2 228 229 0 0 9
chr2 229 230 0 0 2
chr2 262 263 0 0 15
chr2 263 264 0 0 2
chr2 304 305 0 0 13
chr2 305 306 0 0 1
chr2 316 317 0 0 11
chr2 317 318 0 0 1
chr2 318 319 0 0 12
chr2 319 320 0 0 1
chr2 326 327 0 0 11

Do you have any suggestions what could go wrong?

[fkrueger]

Hi gene_x,

not quite sure what is going here or which genome you are aligning your reads to, but for at least human or mouse chromosome 2 the first couple of megabases are masked by Ns, and there is no way that Bismark would map any reads to these sequences or extract methylation calls from it....

[jnhutchinson]

Has anyone here actually tried this?

I'm seeing some strange RRBS results with high non-CpG methylation levels (~4%) in my samples. I suspect it's due to poor(er) conversion, but my client thinks it's due to the non-standard areas we selected for RRBS (larger fragments targetted to non-CpG island areas). I could use a more objective measure of conversion to help decide the issue.

Would appreciate any and all pointers/scripts to get this to work.

Best,
John

[fkrueger]

Hi John,
To see if it is the difference of the captured regions I would suggest looking at non-CG context only in CpG islands. This could be done with a subset of your reads (maybe methylation calls from some 10M reads), import them into SeqMonk, design probes over CGIs and then look at the average methylation (shouldn't take more than 5 mins to find out). Alternatively you could try to look at overall methylation levels on the mitochondrium which is normally very lowly methylated, I am not sure however how well the MT is covered in an RRBS setting... Basically whatever lowest average methylation level you find in any larger aggregate of regions can be considered the upper limit of bisulfite conversion error.

[gene_x]

Originally Posted by fkrueger

Hi gene_x,

not quite sure what is going here or which genome you are aligning your reads to, but for at least human or mouse chromosome 2 the first couple of megabases are masked by Ns, and there is no way that Bismark would map any reads to these sequences or extract methylation calls from it....

Right, there should not be methylation calls on those regions.

A related question, if you have multiple lanes of data (fastq), do you run alignment on each individual lane first and them merge the resulting bam file or do you merge all fastq first and then do a big alignment on the merged fast file?

[fkrueger]

If done correctly either way is fine. I find it easier to merge FastQ up front because that way you can set off a pipeline without having to intervene until you get the final reports.

[jnhutchinson]

Originally Posted by fkrueger

Hi John,
To see if it is the difference of the captured regions I would suggest looking at non-CG context only in CpG islands. This could be done with a subset of your reads (maybe methylation calls from some 10M reads), import them into SeqMonk, design probes over CGIs and then look at the average methylation (shouldn't take more than 5 mins to find out). Alternatively you could try to look at overall methylation levels on the mitochondrium which is normally very lowly methylated, I am not sure however how well the MT is covered in an RRBS setting... Basically whatever lowest average methylation level you find in any larger aggregate of regions can be considered the upper limit of bisulfite conversion error.

Thanks Felix, I'll give those a try. Our RRBS protocol is modified to (ideally) capture non-CpG island regions so I'll look into the chrM methods first.

[gene_x]

Originally Posted by fkrueger

If done correctly either way is fine. I find it easier to merge FastQ up front because that way you can set off a pipeline without having to intervene until you get the final reports.

Hi, Felix,

I've tried a few different ways of doing it but the resulting .zero.cov file still has methylation calls on non-CpG sites. I couldn't figure out how this could happen..

Can you take a look at my pipeline and see if there is anything suspicious?

1. Run adapter trimming and PE bismark(bowtie1) mapping on individual lanes

bismark -n 2 -l 50 --chunkmbs 1024 -X 800 -un --ambiguous --bam $index -1 lane1_read1.fastq -2 lane1_read2.fastq

2. Sort bam files based on read name

samtools sort -n lane1.bam lane1_sort

Do step 1 and step 2 for all lanes.

3. Merge all sorted bam files

samtools merge all_lanes.bam lane1_sort.bam lane2_sort.bam ...

4. remove duplication

deduplicate_bismark -p --bam all_lanes.bam

5. extract methylation calls

bismark_methylation_extractor -p --no_overlap --ignore 3 --ignore_r2 3 --bedGraph --counts --zero_based --report all_lanes.deduplicated.bam

Does this look right to you?

Also I'm wondering if you actually have tried to map individual lanes first and them merge the bam files etc? Is it possible that there might be some hidden bug in the process?

Thank you very much!

[dpryan]

You can probably skip the sorting/merging and just use "samtools cat". The alignments will then be in the same order as they would have had you merged prior to alignment and this will be significantly faster.

[gene_x]

Originally Posted by dpryan

You can probably skip the sorting/merging and just use "samtools cat". The alignments will then be in the same order as they would have had you merged prior to alignment and this will be significantly faster.

Good to know that. Thanks!

[fkrueger]

We have just released a new version of Bismark (v0.13.0), which is available from the Babraham Bioinformatics website. This version adds a couple of useful options and changes some default behavior. Perhaps most notably the methylation extractor may now optionally be run in a multithreaded manner which greatly reduces its processing time (in a preliminary benchmark the elapsed time went down almost linearly when more cores were being used for this process, see below for more details). Here is a list of all changes:

o Bismark: Fixed renaming issue for SAM to BAM files (which would have replaced any occurrence of sam in the file name, e.g. sample1_... instead of the file extension .sam)


o Methylation Extractor: Added new option '--multicore INT' to set the number of cores to be used for the methylation extraction process. If system resources are plentiful this is a viable option to speed up the extraction process (we observed a near linear speed increase for up to 10 cores specified). Please note that a typical process of extracting a BAM file and writing out '.gz' output streams will in fact use ~3 cores per value of --multicore INT specified (1 for the methylation extractor itself, 1 for a Samtools stream, 1 for a GZIP stream), so --multicore 10 is likely to use around 30 cores of system resources. This option has no bearing on the speed of the bismark2bedGraph or genome-wide cytosine report processes

o Methylation Extractor: Added two new options '--ignore_3prime INT' (for single-end alignments and Read 1 of paired-end alignments) and '--ignore_3prime_r2 INT' (for Read 2 of paired-end alignments) to remove positions that display a methylation call bias on the 3' end of reads

o Methylation Extractor: The option --no_overlap is now the default for paired-end data. One may explicitly choose to include overlapping data with the option '--include_overlap'

o Methylation Extractor: The splitting report will now be written out by default (previously optional --report)

o Methylation Extractor: In paired-end mode, read-pairs which had been skipped because either read was shorter than a specified (very high) value of '--ignore' or '--ignore_r2' will now have the information of the other read extracted if it meets the length criteria (if applicable). Thanks to Andrew Dei Rossi for contributing a patch


o bismark2bedGraph: Fixed the location of the sorting directory which could have failed if an output directory had been specified