[jnhutchinson]

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]

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]

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]

Multithreading the methylation extractor

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

[AnneH]

Hi,

I ran Bismark on some of my data and tried subsequently to import the SAM files into CLC genomics workbench.
I encountered some problems since the length of the sequences in the reference genome (which I also have in the CLC and which should be coupled with the SAM file that is to be imported) doesn't match with the length of the sequences reported in the SAM file.

I guess this has something to do with the genome preparation step.
Is there some way to avoid the sequence length to be changed in this process?

[fkrueger]

Hi Anne,

The length of the genome sequences is not changed at any point during the genome preparation, the only thing that does happen is that nucleotides get replaced. The length of the chromosome sequences should be printed to the SAM header (the @SQ lines, eg:
@SQ SN:1 LN:197195432
@SQ SN:10 LN:129993255
@SQ SN:11 LN:121843856
@SQ SN:12 LN:121257530
...).

If these sequences don't match up in CLC, could it be that it is using slighlty different reference sequences and/or genome build?

[AnneH]

Hi Felix,

now that is embarassing, it's true, I used a slightly different reference genome...
My bad!

Thanks for the fast reply though!

[fkrueger]

That's a perfectly fine answer as far as I am concerned, doesn't require me to fix anything :P Good luck!

[gene_x]

bismark report

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.

[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?

[dpryan]

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

[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.

[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,

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 )