[jeni]

Hello,

My bismark report says:

C methylated in CpG context: 52.0%
C methylated in CHG context: 1.0%
C methylated in CHH context: 0.8%
C methylated in unknown context (CN or CHN): 20.6%

Does this means the bisulfite conversion rate is 99% because the non-CpGs methylation rate is around 1%?

Thanks.
Jeni

[fkrueger]

The orientation of a paired-end alignment is given solely by the alignment of Read 1. Since you ran your data in default = directional mode (and rightly so given that you've got a mapping efficiency of 80%!), R1 can only align to the OT and OB strands. Technically R2 aligns to the CTOT and CTOB strands, but since they are part of a paired-end alignment Bismark factors this in automatically and does a G to A lookup for methylation calls of these reads. In short: don't worry, it is all being handled for you!

[chxu02]

Hi Devon,
Thanks for your answer. Still, I'm not quite sure what happened at the methylation calling step. How does bismark process the methylation info from the #2 reads which uniformly reflected the sequence info from the CTOT or CTOB strands?

[dpryan]

Hi Youyou, I'll save Felix from needing to reply.

Internally, bismark is calling bowtie2 (or bowtie) with the --norc flag when it's aligning to the C->T converted genome (this would be the OT strand). Even though bowtie is being told to not align to the reverse complement (CTOT) it'll still reverse complement read #2 of a pair prior to alignment. This is due to the library type flag (--fr) defaulting to what's appropriate for Illumina.

[chxu02]

Hi Felix,
I'm stuck in the claim that "If a library is directional, only reads which are (bisulfite converted) versions of the original top strand (OT) or the original bottom strand (OB) will be sequenced."
I'm quite sure my library is directional. The mapping efficiency is 80% when I use --directional. But I did pair-end sequencing which means that the CTOT & CTOB strands were sequenced from another end. So what happened actually when I ran bismark --directional with my pair-end reads?

Youyou

[fkrueger]

Glad that you got it sorted. If your libraries were derived from a ChIP-ped sample I would agree that you probably don't want to deduplicate the data. Unless you used a unique molecule identifier for each read there is no way to discriminate between genuine ChIP-ped fragments and PCR amplified reads, but deduplication would completely take the depth out of your experiment.

[chxu02]

Hi Youyou, I don't think file size is the best way to judge the outcome your experiment. Rather, it might be helpful to know whether all steps finished successfully, the sequence mode and length, whether the reads were trimmed, the mapping efficiency, duplication rate etc. Ideally could you attach here or send me the reports (FastQC, trimming, mapping, dedupliciation and methylation extraction reports) via email? Cheers, Felix

Hi Felix,

Thanks for your prompt reply. Now I figured it out. The distribution of reads from the library with smaller bedgraph is further enriched at some certain genomic regions compared to another one. That's why two libraries ended up with similar amount of C mapped, while distinct numbers of lines in bedgraph generated.

I have one more concern about the deduplication step. Each of my BS-seq libraries was made from 1ng ChIP DNA. So I don't expect to see genome-wide distribution of reads, and with >200M reads, I can tolerate some level of PCR duplication which should be determined by genome size and library depth, while the bismark deduplication script seems to cut the number of duplicated reads arbitrarily to 1, which vastly reduced the depth of my libraries. I just wonder if there is any alternative for me to do the deduplication step.

Appreciate your contribution.

Youyou

[fkrueger]

Hi Youyou, I don't think file size is the best way to judge the outcome your experiment. Rather, it might be helpful to know whether all steps finished successfully, the sequence mode and length, whether the reads were trimmed, the mapping efficiency, duplication rate etc. Ideally could you attach here or send me the reports (FastQC, trimming, mapping, dedupliciation and methylation extraction reports) via email? Cheers, Felix

[chxu02]

much smaller bedgraph & cov files generated...

Hi all,

I got some results I can not explain. I have 2 human BS-seq libraries. Before sequencing with great depth, I sequenced a little bit to get 1M raw reads for each library, and analyzed the data with bismark. It all looked good. The methylation extractor files were with similar size, and bedgraph& cov files with similar size too. Now I sequenced to get 250M PE raw reads for each library. The methylation extractor files are with similar size, while the bedgraph and cov files from one library are with only 1/4 of the size of files from another library. When imported into IGV, it looks like I have much less coverage for the library with smaller files. Because of the biological information in the 2 samples, I expect to see same coverage while different methylation levels. What's going on here?

Youyou

[fkrueger]

Hopefully last release for the year: Bismark v0.13.1

We have just released a new version of Bismark (v0.13.1), which is available from the Babraham Bioinformatics website. This version adds several new useful options and bug fixes, perhaps most notably the new option --merge_CpG for the coverage2cytosine module that pools CpG information from both top and bottom strand methylation evidence and a fix for the methylation extractor option --multicore which now accurately shows the correct number of events in the M-bias plots.

Here is a more detailed list of all changes:
• Bismark Genome Preparation: Added a check for unique chromosome names to the Bismark indexer to avoid disappointments later
• Methylation Extractor: Fixed a bug for the M-bias reports when the option --multicore was used, in which case only the numbers of one core were used to constuct the report. Now every different thread writes out an individual M-bias table, and once the methylation extraction has completed all these individual files are merged into a single, cumulative table as it should be
• Methylation Extractor: Added a new option --mbias_off, which processes the files as normal but does not write out any M-bias files. This option is meant for users who run the methylation extractor two times, the first time to figure out whether there is a bias that needs to be removed, and the second time using the --ignore options, but without overwriting the already existent M-bias reports
• bismark2bedGraph: Deferred removal of the input file path information a little so that specifying file paths doesn't prevent bismark2bedGraph from finding the input files anymore
• bismark2bedGraph: If the specified output directory doesn't exist it will be created
• bismark2bedGraph: Changed the way scaffolds are sorted (with --gazillion/--scaffold specified) to -k3,3V (this was done following a suggestion by Volker Brendel, Indiana University: "The -k3,3V sort option is critical when the sequence names are numbered scaffolds (without left-buffering of zeros)
• coverage2cytosine: Added a new option --merge_CpG that will post-process the genome-wide report to write out an additional coverage file which has the top and bottom strand methylation evidence pooled into a single CpG dinucleotide entity. This may be the desirable input format for some downstream processing tools such as the R-package bsseq (by K.D. Hansen). An example would be:

  genome-wide CpG report (old)
         gi|9626372|ref|NC_001422.1|     157     +       313     156     CG
         gi|9626372|ref|NC_001422.1|     158     -       335     156     CG

  merged CpG evidence coverage file (new)
         gi|9626372|ref|NC_001422.1|     157     158     67.500000       648     312

This option is currently experimental, and only works if CpG context only and a single genome-wide report were specified (i.e. it doesn't work with the options --CX or --split_by_chromosome)
• coverage2cytosine: Changed the processing of not-covered chromosomes so that they are sorted and not processed randomly. This should make runs more reproducible

Comments or suggestions welcome. Happy New Year! Felix

[fkrueger]

Hi,

I'm wondering if I run several instances of bismark alignment separately and then later merged the bam files together. I'd like to get a report on the final bam file, like the PE report file. It's required to run bismark2report command, so I'm wondering if there is a program to run just that, say, given a bam file and get the necessary info from it as a stand-alone command, just like bismark2report? Currently, I don't see an easy way to get this.. Could it be a good addition to the next release?

Thanks.

What you are looking for would be something that reads several mapping reports, sums up the numbers and generates one overall report which I am afraid doesn't currently exist in this form. It is also not possible to get all the information just from the merged BAM file because some of the information never actually makes it into the BAM file(s).
What you could do in the meantime though is sum the numbers from 2 or more files up e.g. in Excel and then provide this 'merged' report file to bismark2report using the option --alignment_report FILE. I am still aiming to implement an option that will perform multiple alignments as well as automatic merging once everything has completed, but I wasn't planning to write something to merge reports at the current time (something that sounds very trivial but might take quite some time )

[Wonghe]

Thank you for the suggestion. We have trimmed both the 5' and 3' ends to remove base with poor quality. Hope this will help increase the overall mapping rate.

[fkrueger]

Hi,
We recently encounter very poor mapping efficiency for our RRBS libraries sequenced with Hi-Seq Rapid mode 1x51bp. Using the default parameter in Bismark we only get 9-45% mapping efficiency. We do observe that about 10% of the reads had NN bases in first 3 bases of the reads. Would this affect the mapping rate and any suggestion to increase the efficiency?

While Ns in the read do not count towards the maximum allowed error score for use with Bowtie, they do count for the mismatch parameter -n (so -n 1, which is the default, should remove reads with 2 NN in the seed region). If you use bowtie 2 (--bowtie2) the Ns should only count as relatively small penalty (of -1), which might already improve your mapping a bit. I would agree with Devon though that you should probably best trim the reads from the 5' end, in addition to the adapter, quality and RRBS trimming from the 3' end. A command for Trim Galore could look like this:

trim_galore --clip_r1 2 --rrbs file.fastq.gz

More information on the recommended RRBS trimming can be found here.

[dpryan]

@Wonghe: Try performing 5' quality trimming. That might increase the mappability a bit.

[Wonghe]

low mapping rate in SR reads

Hi,
We recently encounter very poor mapping efficiency for our RRBS libraries sequenced with Hi-Seq Rapid mode 1x51bp. Using the default parameter in Bismark we only get 9-45% mapping efficiency. We do observe that about 10% of the reads had NN bases in first 3 bases of the reads. Would this affect the mapping rate and any suggestion to increase the efficiency?