[PeteH]

Hi Felix,

Currently the SAM header produced by Bismark has the @SQ fields in a more-or-less random order. This means that the @SQ fields will almost certainly not be in the same order as the they are in the reference genome/index. It's not a huge deal but it means I have to re-header my SAM files because I rely on the order of the @SQ fields in the SAM header in downstream tools. Is it possible to change this behaviour so that the @SQ fields in the SAM header are in the same order as in the reference genome/index?

My understanding is that the @SQ information is stored in a hash ("chromosomes") and so there is no guarantee on the order that these will be printed in when called by the function generate_SAM_header(). I'm not sure how hard that is to do this in the current framework. My idea was to store an index variable as an additional value in that hash, where the index is the position of that chromosome in the reference genome/index, but my limited Perl skills couldn't get it working.

I understand if you don't consider it an issue and I will just continue to re-header my SAM files.

Thanks,
Pete

[dariober]

Regarding picard I was under a probably false impression that samtools and picard are using the same algorithm to eliminate duplicates... I had a difficulty to find how picard defines PE duplicate (which is clearly explained in deduplicate_bismark)

Hi- See this link for an explanation of MarkDuplicates. if you haven't already seen it http://sourceforge.net/apps/mediawik...icates_work.3F

Q: How does MarkDuplicates work?

A: Essentially what it does (for pairs; single-end data is also handled) is to find the 5' coordinates and mapping orientations of each read pair. When doing this it takes into account all clipping that has taking place as well as any gaps or jumps in the alignment. You can thus think of it as determining "if all the bases from the read were aligned, where would the 5' most base have been aligned". It then matches all read pairs that have identical 5' coordinates and orientations and marks as duplicates all but the "best" pair. "Best" is defined as the read pair having the highest sum of base qualities as bases with Q >= 15.

If your reads have been divided into separate BAMs by chromosome, inter-chromosomal pairs will not be identified, but MarkDuplicates will not fail due to inability to find the mate pair for a read.

And here http://www.biostars.org/p/3917/ for some comments about rmdup vs MarkDuplicates.

Dario

[linnet]

duplicates

Dear fkrueger and Dario,

Thank you both for lighting replies and suggestions! I will try picard, since I want to have more permanent solution rather than bothering fkruager every time I have a big dataset. If I am completely stuck, I may be in touch though.

Regarding picard I was under a probably false impression that samtools and picard are using the same algorithm to eliminate duplicates... I had a difficulty to find how picard defines PE duplicate (which is clearly explained in deduplicate_bismark)

Thanks again,
Skirmantas

[dariober]

Alternatively I have tried samtools sort; rmdup; sort by name, but then methylation extractor is not happy, probably due to some pair-mates removed.

Hi- I think samtools rmdup doesn't handle paired reads correctly. I would try to use picard MarkDuplicates instead.

Good luck!
Dario

[fkrueger]

Hi Skirmantas,

Deduplicating such a large file with only 48GB of memory is probably not going to work... Admittedly, Perl is probably not the best option to store and quickly look up such an extraordinary number of entries, but the deduplication procedure is storing a composite string of chromosome:start:end:strand which is sort of essential to the process. One could probably reduce the memory footprint somewhat by not including the chromosome into the string but using a different hash for each chromosome, and replacing the end coordinate by the length of the fragment, but I don't know how much this would reduce the footprint.
Alternatively I could offer you to set up an FTP site for your file and deduplicate it here overnight ...

[linnet]

deduplicate_bismark

Hello,

Does anyone knows input file size to memory footprint of deduplicate_bismark ? I have 156 GB WGBS bam file (10^9 PE reads) and I suspect 48 GB workstation is getting into HDD swap usage, which is not great... Any ideas or solutions? Alternatively I have tried samtools sort; rmdup; sort by name, but then methylation extractor is not happy, probably due to some pair-mates removed.

Thanks,
Skirmantas

[dpryan]

For Bison, I wrote a bedGraph2Methylkit program. You'd actually need to download the full bison source code (as well as samtools) for compilation, but that's probably the fastest route if you've already used the methylation extractor to generate bedGraph files (I've never tried this on Bismark-generated bedGraph files, so YMMV).

[kalyankpy]

Dear Chris,

Easiest way to move from Bismark to Methylkit is to use the .sam file from alignment as input for Methylkit.

In the other case one has to extract methylation data from .sam file and convert it into bed or bedgraph file. One has to process further if you intend to merge the strand coverage on either strands. Lastly one has to convert the merged strand coverage data into methylkit format. I have been following this method for few samples. I have a script written to do this job. This script was written for an earlier version of Bismark (0.7.5) and all the paths are hardcoded (works only on my computer). Contact me if you need some help in using this.

So, among the two methods, easiest one would be to give .sam file as input for methylkit. Somehow, I felt the methylkit tutorial/vignette is not upto the standard of other bioconductor packages. In case you are familiar with R, it will be an easy ride!

Kalyan

[cjsifuen]

Hi all,

Sorry if this question is a bit simple, I am very new to all of this. I am trying to reanalyze some RRBS data, following methods from a paper that was recently published. I used bismark to align paired-end RRBS reads, then used the bismark_methylation_extractor to make methylation calls. Now, I am trying to use methylKit to perform further analysis, but I cannot figure out how to get the methylation call file with the format that methylKit asks for. Is there a way to get the output file with only with the format below?

chrBase chr base strand coverage freqC freqT
chr21.9764539 chr21 9764539 R 12 25.00 75.00
chr21.9764513 chr21 9764513 R 12 0.00 100.00
chr21.9820622 chr21 9820622 F 13 0.00 100.00
chr21.9837545 chr21 9837545 F 11 0.00 100.00
chr21.9849022 chr21 9849022 F 124 72.58 27.42
chr21.9853326 chr21 9853326 F 17 70.59 29.41

I am trying to repeat what the group in the paper did, so I would prefer to read in the bismark methylation call files and not the bismark alignment files.

Thanks,

Chris

[fkrueger]

Oh, too late...

[fkrueger]

Hi Brent,

Using --bam will output a bam file, or, if there is no samtools available, make a sam.gz file.

Cheers,
Felix

[brentp]

... and now I see that it already does write directly to bam. thanks.

[brentp]

comments most welcome.

Hi, Bismark is creating a lot of plain-text output--e.g. 16GB of SAM from 4GB of fastq.gz. It would be nice if it output BAM format instead of SAM. Even if it could be forced to output .sam.gz, that would help those of us on clusters with limited disk space.

[fkrueger]

We have just made available a new version of Bismark (v0.10.1) which fixes a few issues and adds some useful support for unfinished genomes with lots of scaffolds instead of just a handful of chromosomes. Here are all the changes in more detail:

• Bismark methylation extractor: The methylation extractor does now detect automatically whether Bismark alignment file(s) were run in single-end or paired-end mode. The automatic detection can be overridden by manually specifying -s or -p and this option is only available for SAM/BAM files

• bismark2bedGraph: When run in stand-alone mode, the coverage file will replace 'bedGraph' as the file ending with 'bismark.cov'. If the output filename is anything other than 'bedGraph', '.bismark.cov' will be appended to the filename
• bismark2bedGraph: When run in stand-alone mode, '--counts' will be enabled by default for the coverage output
• bismark2bedGraph: Added a new option '--scaffolds/--gazillion' for users working with unfinished genomes sporting tens or even hundreds of thousands of scaffolds/contigs/chromosomes. Such a large number of reference sequences frequently resulted in errors with pre-sorting reads to individual chromosome files because of the operating system's limitation of the number of filehandles that can be written to at any one time (typically this limit is anything between 128 and 1024 filehandles; to find out this limit on Linux, type: ulimit -a). To bypass the limitation of open filehandles, the option '--scaffolds' does not pre-sort methylation calls into individual chromosome files. Instead, all input files are temporarily merged into a single file (unless there is only a single file), and this file will then be sorted by both chromosome AND position using the UNIX sort command. Please be aware that this option might take a looooong time to complete, depending on the size of the input files, and the memory you allocate to this process (see '--buffer_size')
• bismark2bedGraph: Added a new option '--ample_memory'. Using this option will not sort chromosomal positions using the UNIX sort command, but will instead use two arrays to sort methylated and unmethylated calls, respectively. This may result in a faster sorting process for very large files, but this comes at the cost of a larger memory footprint (as an estimate, two arrays of the length of the largest human chromosome 1 (~250 million bp) consume around 16GB of RAM). Note however that due to the overhead of creating and looping through huge arrays this option might in fact be *slower* for small-ish files (up to a few million alignments). Note also that this option is not currently compatible with options '--scaffolds/--gazillion'. This option still needs some efficiency testing as to when it actually makes sense to use it, but it produces identical results to the default sort option. Thanks to Yi-Shiou Chen for contributing this twist

• deduplicate_bismark: The deduplication script does now detect automatically whether a Bismark alignment file was run in single-end or paired-end mode (this happens separately for every file analysed). The automatic detection can be overridden by manually specifying -s or -p and this option is only available for SAM/BAM files

• bismark2report: Specifying a single file for each of the optional reports does now will now work as intended, instead of being skipped

• coverage2cytosine: Added some counting and statements to indicate when the run finished successfully (it proved to be difficult to follow the report process for a genome with nearly half a million scaffolds...)

Bismark is available from http://www.bioinformatics.babraham.a...jects/bismark/. Bug reports or comments most welcome.

[fkrueger]

We have just released a new version of Bismark (v0.10.0) that adds a variety of convenient features and bug fixes. The changes in detail are:

• Bismark: The option '--prefix' does now also work for the C->T and G->A transcribed temporary files to allow multiple instances of Bismark to be run on the same file(s) in the same folder (e.g. using Bowtie and Bowtie 2 or some stricter and laxer parameters concurrently)

• bismark2report: Changed the behavior of this module to automatically find all Bismark mapping reports in the current working directory, and to try and work out whether the optional reports are present as well (i.e. deduplication, splitting and M-bias reports). This uses the file basename and will fail if the files have been renamed at any stage
• bismark2report: Added commas as separator for large numbers to improve readability

• Bismark methylation extractor: will now delete unsused methylation context files (e.g. CTOT and CTOB files for a directional library)

• bismark2bedGraph: Dropped the option -k3,3 from the sort command to result in a dramatic speed increase while sorting. This option had been used previously to enable sorting by chromosome in addition to position, but should no longer be needed because the files are being read in sorted by chromosome already
• bismark2bedGraph: This module does now produce these two output files:
(1) A bedGraph file, which now contains a header line: 'track type=bedGraph'. The genomic start coords are 0-based, the end coords are 1-based. These changes should make the file truly compatible with the UCSC genome browser.
(2) A coverage file ending in .cov. This file replaces the former 'bedGraph --counts' file and is required to proceed with the subsequent step to generate a genome-wide cytosine report (the module doing this has been renamed to coverage2cytosine to reflect this file name change)

• coverage2cytosine: Changed the name of this module from 'bedGraph2cytosine' to 'coverage2cytosine' to reflect the change that this module now requires the methylation coverage file produced by the bismark2bedGraph module (this coverage file replaces the former "bedGraph --counts" output)
• coverage2cytosine: Previously, the cytosine report would always report every C position in any context, even though the default should have reported CpG positions only. This has now been fixed

• Bismark genome preparation: Made a couple of changes to make the genome preparation fully non-interactive. This means that the path to the genome folder and to Bowtie (1/2) have to be specified up front (for Bowtie (1/2) it is otherwise assumed that it is in the PATH). Furthermore, already existing bisulfite indices in the target folder will be overwritten and the user is no longer prompted if he agrees to this. We got rid of this because creating a second index (Bowtie 1 as well as 2) in the same folder in non-interactive mode got stuck in loops asking whether it is alright to proceed or not, generating therabyte sized log files without ever starting doing anything useful...)

• deduplicate_bismark: Renamed the rather long deduplication script to this slightly shorter one. Also added some filehandle closing statements that might have caused buffering issues under certain circumstances

Bismark is available from https://www.bioinformatics.babraham....jects/bismark/