[fkrueger]

Hmm not quite sure what the problem is but to me it looks like all the information is right there:

HWI-ST571:431:C43FMACXX:8:2316:9928:82689_1:N:0:GTCCGC - chr9 23913539 z

HWI-ST571:431:C43FMACXX:8:2316:9928:82689_1:N:0:GTCCGC is the Read ID
- is the methylation state (+ is methylated, - is unmethylated)
chr9 is the chromosome, here chromosome 9
23913539 is the mapping position of the C in question
z is the context, so here CpG context

Or am I missing something?

[jcrow06]

Hi there,

So I have recieved some Bismark output files:

Bismark (SAM/BAM)

Mbias files & Plots

Methylation Extractor reports



The documentation is very good and I have a good understanding of all the files however in my data from the Methylation extractor there are several sites that are not marked with the chromosome header.

The report files according to methylationExractor should have the following headers:

(1) seq-ID (2) methylation state (3) chromosome (4) start position (= end position) (5) methylation call

In my files I have as follows:

HWI-ST571:431:C43FMACXX:8:2316:9928:82689_1:N:0:GTCCGC - chr9 23913539 z
HWI-ST571:431:C43FMACXX:8:2316:9928:82689_1:N:0:GTCCGC + chr9 23913531 Z
HWI-ST571:431:C43FMACXX:8:2316:9928:82689_1:N:0:GTCCGC + chr9 23913529 Z
HWI-ST571:431:C43FMACXX:8:2316:9928:82689_1:N:0:GTCCGC - chr9 23913485 z
HWI-ST571:431:C43FMACXX:8:2316:9928:85484_1:N:0:GTCCGC + chr2 28594224 Z
HWI-ST571:431:C43FMACXX:8:2316:9928:85484_1:N:0:GTCCGC - chr2 28594341 z
HWI-ST571:431:C43FMACXX:8:2316:9928:85757_1:N:0:GTCCGC - 30126 6594 z
HWI-ST571:431:C43FMACXX:8:2316:9928:85757_1:N:0:GTCCGC - chr2 223724995 z
HWI-ST571:431:C43FMACXX:8:2316:9929:12206_1:N:0:GTCCGC - chr12 64474195 z
HWI-ST571:431:C43FMACXX:8:2316:9929:59562_1:N:0:GTCCGC + chr20 55837341 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:59562_1:N:0:GTCCGC - chr20 55837591 z
HWI-ST571:431:C43FMACXX:8:2316:9929:59562_1:N:0:GTCCGC + chr20 55837557 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:59562_1:N:0:GTCCGC + chr20 55837519 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:89987_1:N:0:GTCCGC + chr5 92603043 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:89987_1:N:0:GTCCGC + 26965 45121 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:89987_1:N:0:GTCCGC + 26965 301 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:89987_1:N:0:GTCCGC + chr5 92603089 Z
HWI-ST571:431:C43FMACXX:8:2316:9929:89987_1:N:0:GTCCGC + chr5 92603122 Z

I am not sure about this as they seem to be a strange format. Initially I though maybe there was a formatting issue and it skipped the "chr" but none of the numbers correspond to a CG site.

When I run bismark2bedgraph the job completes and these files are also listed but without the chromosomal location I cannot locate the methylation site.

Have any of you observed this before or know what has happened? Or could it be a small bug?

Cheers

[fkrueger]

Hi! I work in the Jaffe Lab, and I'm using the new version of bismark that was attached in the message above. I'm having problems with it. I sampled the first million reads, and I'm aligning them to the human genome. I'm using --multicore 4, and submitted the job to 20 cores, 2Gb each (max allowable memory per core -- 3Gb). The program has been running for over 19 hours, its using 61Gb of scratch space, and has been on the following stage for the past 16 hours:

Now starting a Bowtie paired-end alignment for GAread1CTread2GAgenome (reading in sequences from /scratch/temp/4549202.1.shared.q/WGC025593L_combined_R1.fastq.gz.temp.4.gz.GA_plus_CT.fastq.gz, with the options: -q -n 1 -k 2 --best --maxins 500 --chunkmbs 512 --norc)

This seems slow. I was wondering if you could help. Thanks in advance!

Yes this is unacceptably slow, this process shouldn't take more than 20 seconds or so... I am not quite sure if the memory allocation is the problem, but I wouldn't be surprised. Just as a guideline: Each instance of Bismark will need ~3GB to hold the genome in memory, and around another 3-4GB per alignment thread, so 2 times this for default directional, and 4 times this for non-directional mapping. The other processes (i.e. Samtools or gzip streams etc) shouldn't need much extra memory. Since you specified --multi 4 you will need this memory times four again. so roughly 40-50GB for directional or ~60GB for non-directional alignments. Is there are a chance that your cluster just won't execute the commands properly?
As a solution you could either try to ask for enough resources or go down with --multi, but each instance of Bismark will need ~10-16GB either way. Cheers, Felix

[nivanov]

Hi! I work in the Jaffe Lab, and I'm using the new version of bismark that was attached in the message above. I'm having problems with it. I sampled the first million reads, and I'm aligning them to the human genome. I'm using --multicore 4, and submitted the job to 20 cores, 2Gb each (max allowable memory per core -- 3Gb). The program has been running for over 19 hours, its using 61Gb of scratch space, and has been on the following stage for the past 16 hours:

Now starting a Bowtie paired-end alignment for GAread1CTread2GAgenome (reading in sequences from /scratch/temp/4549202.1.shared.q/WGC025593L_combined_R1.fastq.gz.temp.4.gz.GA_plus_CT.fastq.gz, with the options: -q -n 1 -k 2 --best --maxins 500 --chunkmbs 512 --norc)

This seems slow. I was wondering if you could help. Thanks in advance!

[fkrueger]

Again the problem here was that the filenames were specified with /full/path/information/ that broke things. I had fixed this for FastQ files but you were using FastA files... Please find attached a version that should be working now.

Regarding your other question about parallelisation: I haven't done extensive testing on this but the -p parameter only allows Bowtie 2 to produce the alignments somewhat faster, whereas --multicore will should almost fully parallel. The -p speed increase seems to be plateauing fairly quickly, possibly also because of the unique alignment determination and methylation calling processes need some time, too. If you've got enough system resources I would certainly recommend using --multi over -p (or potentially a combination thereof) in terms of speed increase.

[ajaffe]

--multicore and --temp_dir still crashing

Hi,

We are still having similar issues when specifying "--multicore" and "--temp_dir" in the same run. Note this is not a problem when specifying "--multicore" but not "--temp_dir" or "--temp_dir" but not "--multicore". Eg:

> bismark --temp_dir Pilot/tmp --non_directional --bowtie2 --multicore 4 -f -o Pilot/tmp lambda/ -1 /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta -2 /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_2.fasta

...
The provided filenames for paired-end alignments are /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta and /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_2.fasta
The provided filenames for paired-end alignments are /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta and /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_2.fasta
Failed to write output /dcs01/ajaffe/Brain/DNAm/WGBS/Pilot/tmp//dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta.temp.1: No such file or directory
Failed to write output /dcs01/ajaffe/Brain/DNAm/WGBS/Pilot/tmp//dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta.temp.4: No such file or directory
Failed to write output /dcs01/ajaffe/Brain/DNAm/WGBS/Pilot/tmp//dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta.temp.2: No such file or directory
Failed to write output /dcs01/ajaffe/Brain/DNAm/WGBS/Pilot/tmp//dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta.temp.3: No such file or directory

> bismark --non_directional --bowtie2 --multicore 4 -f -o Pilot/tmp lambda/ -1 /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta -2 /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_2.fasta

...

====================
Bismark run complete
====================

> bismark --non_directional --bowtie2 --temp_dir Pilot/tmp/ -f -o Pilot/tmp lambda/ -1 /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_1.fasta -2 /dcs01/ajaffe/junctionCountsData/sim1/Reads/run1_cov1/sample_01_2.fasta

...

====================
Bismark run complete
====================


>

I was also wondering how --multicore differed from the -p flag with bowtie2, e.g. if I specify -p 4 keeping bismark itself on 1 core, will this have the same linear speed up as specifying --multicore 4 and -p 1?

Thanks!!

[bruce01]

Thanks Felix, I have tried BisSNP but it runs very slowly (so slow as to be unusable, trying to figure out why currently) and have read the MethylExtract paper but not yet gotten around to running it.

[fkrueger]

It appears that the coverage in BS-Seq can be affected by several factors in addition it simply being a harsh chemical treatment, e.g. the conditions used for the bisulfite-treatment, polymerase bias during amplification and the sequence composition of the single-stranded DNA. So yes maybe this is not the ideal way of looking at CNVs in these samples.

There are to my knowledge currently two tools out there that attempt to look at SNVs in bisulfite data and correct for this. One is called Bis-SNP (http://epigenome.usc.edu/publicationdata/bissnp2011/), the other one MethylExtract (http://bioinfo2.ugr.es/MethylExtract/). Maybe it would be worth looking into them?

[bruce01]

Hi Felix,

thanks for the reply, sorry I missed the section in the manual.

I have tried using methylation data for CNV, and comparing this to shallow WGS and also exome for same samples. Exome, WGS are consistent, but methylseq is not (it is whole genome, from a MiSeq QC run to determine conversion efficiency). We also have captured methylseq (CpGiant) which is also inconsistent with WGS, exome.

The 'pre-bisulfite state' I refer to is before the DNA is bisulfite-treated, so yes, essentially the reference genome, albeit with mutations etc. My boss is thinking this can be used in SNP, indel calling, as well as CNV. Essentially instead of sequencing exome. I am unconvinced due to bisulfite treatment damaging DNA, but said I would try.

Anyway, thanks for your help and sorry if this is obscure, and my nomenclature incorrect/confusing,

Bruce.

[fkrueger]

Hi Felix,

firstly, thanks for Bismark, it is a great tool. And secondly, I have read the entire thread here and searched elsewhere, so apologies if this is discussed elsewhere.

My issue is understanding the XM:Z tag in Bismark output which I believe indicates methylation status (from a line in the change log). I would like to be able to interpret this tag (i.e. what do x, h indicate?), to the end that possibly the sequence could be 'reconverted' to pre-bisulfite conversion (PI is trying to 'add value' to our methylseq...). Your thoughts on how useful this would be would also be of interest, specifically for copy number analysis (using off-target sequence with CopywriteR package, for example).

Thanks,

Bruce.

Hi Bruce,

First allow me to paste a passage from the User Guide:

The methylation call string contains a dot ‘.’ for every position in the BS-read not involving a cytosine, or contains one of the following letters for the three different cytosine methylation contexts (UPPER CASE = METHYLATED, lower case = unmethylated):

z	unmethylated C in CpG context
Z	methylated C in CpG context 
x	unmethylated C in CHG context
X	methylated C in CHG context
h	unmethylated C in CHH context
H	methylated C in CHH context
u	unmethylated C in Unknown context (CN or CHN)
U	methylated C in Unknown context (CN or CHN)

I am not quite sure I understand your other question fully though, what exactly do you mean with a pre-bisulfite state? In theory it is of course possible to do various things with the read, for which you need to take the reference sequence, the CIGAR string as well as the mismatch string into consideration.

Can't you just take the BAM file and use this to do copy number analyses? Using the coordinates of the read alignmnents you can relatively quickly extract the fully unconverted sequence from the reference genome (is this the pre-bisulfite state?), or go back and extract the bisulfite sequence as it comes from the sequencer out from the FastQ file (e.g. based on the read ID). Does this answer your questions?

[bruce01]

Hi Felix,

firstly, thanks for Bismark, it is a great tool. And secondly, I have read the entire thread here and searched elsewhere, so apologies if this is discussed elsewhere.

My issue is understanding the XM:Z tag in Bismark output which I believe indicates methylation status (from a line in the change log). I would like to be able to interpret this tag (i.e. what do x, h indicate?), to the end that possibly the sequence could be 'reconverted' to pre-bisulfite conversion (PI is trying to 'add value' to our methylseq...). Your thoughts on how useful this would be would also be of interest, specifically for copy number analysis (using off-target sequence with CopywriteR package, for example).

Thanks,

Bruce.

[fkrueger]

New Bismark release: v0.14.3

We have just released a new version of Bismark (v0.14.3) to address the following issues:

• Bismark: Changed the renaming settings for paired-end files so that 'sam' within the filename no longer gets renamed to 'bam' (e.g. smallsample.sam > smallbample.sam)

• Bismark: Input files specified with filepath information are now handled properly in --multicore runs

• Bismark: The --multicore option currently requires the output files to be in BAM format, so specifying --sam at the same time has been disallowed

• Methylation Extractor: fixed another bug for the same issue as in 0.14.1 that had crept into the 0.14.2 release (to do with --ignore_3prime)

• coverage2cytosine: Changed the option --merge_CpG so that CGs starting at position 1 are not considered (since the 3-base sequence context of the bottom strand C at position 2 can not be determined)

Bismark is available for download from the Babraham Bioinformatics Project website.

[enrico16]

Bismark v0.14.2 fails to delete and merge files in --multicore mode

Hi Felix,

I think I might have found a bug in bismark v0.14.2 when run in --multicore mode at the cleaning up and merging stages.

Bismark fails to:

• delete the temporary sequence files (*.fastq.gz.temp.N)
• merge the BAM files (*.fastq.gz.temp.N_bismark_bt2_pe.bam)
• delete temporary BAM files (*.fastq.gz.temp.N_bismark_bt2_pe.bam)
• write the final report text file (*.fastq.gz_bismark_bt2_PE_report.txt)

Please see the full log here: http://pastebin.com/wP3SphjK

I think this happens when the --temp_dir argument is specified. From what I can tell Bismark tries to find the above files in the wrong location and mixes the input (02trimming) and the temp/output (05alignment) folders up a bit by concatenating them.

For example, when running bismark in this way:

bismark --multicore 2 --bowtie2 -p 2 --temp_dir /home/enrico16/analysis/05alignment/SampleX --output_dir /home/enrico16/analysis/05alignment/SampleX /home/enrico16/analysis/00genome -1 /home/enrico16/analysis/02trimming/SampleX/SampleX-trimmed-pair1.fastq.gz -2 /home/enrico16/analysis/02trimming/SampleX/SampleX-trimmed-pair2.fastq.gz

I get this type of errors:

/home/enrico16/analysis/02trimming/SampleX/SampleX-trimmed-pair1.fastq.gz.temp.1	Failed to delete temporary FastQ file 

/home/enrico16/analysis/05alignment/SampleX//home/enrico16/analysis/02trimming/SampleX/SampleX-trimmed-pair1.fastq.gz.temp.1: No such file or directory

Can you help me figure out what's happening here? I might well be doing something wrong. If not, can you please fix this in your code?

Thanks a lot!

Best,
Enrico

[fkrueger]

Thanks for that, I am sure some renaming settings are going wild here. I'll take a look once I am back from holiday. Cheers, Felix

[jhooge]

bismark renames alignment file

Hello Felix,

I found a glitch in bismark 0.14.2, which made me smile a bit.
I ran bismark with the following arguments:

bismark --quiet --bowtie2 -N 1 -p 32 -X 1000 -o . -B smallsample1 /custom_indexes/hg19canon_bsgenome -1 ./smallsample1_R1_val_1.fq.gz -2 ./smallsample1_R2_val_2.fq.gz

Given the basename "smallsample1" I would expect the two following output files:

smallsample1_pe.sam
smallsample1_PE_report.txt

What I get however is:

smallbample1_pe.sam
smallsample1_PE_report.txt

I have now idea, but this looks like a string formatting mistake somewhere in your code. Or is there something I'm missing?

If it helps, here is the log of my bismark call:

Path to Bowtie 2 specified as: bowtie2
Output format is BAM (default)
Alignments will be written out in BAM format. Samtools found here: 'software/SAMtools/1.1-foss-2014b/bin/samtools'
Reference genome folder provided is /custom_indexes/hg19canon_bsgenome/   (absolute path is '/custom_indexes/hg19canon_bsgenome/)'
FastQ format assumed (by default)
Attention: using more than 4 cores per alignment thread has been reported to have diminishing returns. If possible try to limit -p to a value of 4
Each Bowtie 2 instance is going to be run with 32 threads. Please monitor performance closely and tune down if needed!
Library is assumed to be strand-specific (directional), alignments to strands complementary to the original top or bottom strands will be ignored (i.e. not performed!)
Output will be written into the directory: /pipelines/trimmed/
Setting parallelization to single-threaded (default)

Single-core mode: setting pid to 1

Paired-end alignments will be performed
=======================================

The provided filenames for paired-end alignments are ./smallsample1_R1_val_1.fq.gz and ./smallsample1_R2_val_2.fq.gz
Input files are in FastQ format
Writing a C -> T converted version of the input file smallsample1_R1_val_1.fq.gz to smallsample1_R1_val_1.fq.gz_C_to_T.fastq

Created C -> T converted version of the FastQ file smallsample1_R1_val_1.fq.gz (47 sequences in total)

Writing a G -> A converted version of the input file smallsample1_R2_val_2.fq.gz to smallsample1_R2_val_2.fq.gz_G_to_A.fastq

Created G -> A converted version of the FastQ file smallsample1_R2_val_2.fq.gz (47 sequences in total)

Input files are smallsample1_R1_val_1.fq.gz_C_to_T.fastq and smallsample1_R2_val_2.fq.gz_G_to_A.fastq (FastQ)
Now running 2 instances of Bowtie 2 against the bisulfite genome of /custom_indexes/hg19canon_bsgenome/ with the specified options: -q -N 1 --score-min L,0,-0.2 -p 32 --reorder --ignore-quals --no-mixed --no-discordant --maxins 1000 --quiet

Now starting a Bowtie 2 paired-end alignment for CTread1GAread2CTgenome (reading in sequences from smallsample1_R1_val_1.fq.gz_C_to_T.fastq and smallsample1_R2_val_2.fq.gz_G_to_A.fastq, with the options: -q -N 1 --score-min L,0,-0.2 -p 32 --reorder --ignore-quals --no-mixed --no-discordant --maxins 1000 --quiet --norc))
Found first alignment:
M00978:20:000000000-A5TV0:1:1101:15681:1332_1:N:0:3/1   99      chr16_CT_converted      70485744        42      26M     =       70485880        162  TTATTGGTGGTTGGGTATTGTGGTTT       >AABBFFFDBBB222EEFGGGFFGHF      AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:26 YS:i:0  YT:Z:CP
M00978:20:000000000-A5TV0:1:1101:15681:1332_2:N:0:3/2   147     chr16_CT_converted      70485880        42      26M     =       70485744        -162 AAAAATTAGTTGGGTGTGGTGGTGGG       GFG3FFGGGGGGFEA1AFADDA1A1>      AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:26 YS:i:0  YT:Z:CP
Now starting a Bowtie 2 paired-end alignment for CTread1GAread2GAgenome (reading in sequences from smallsample1_R1_val_1.fq.gz_C_to_T.fastq and smallsample1_R2_val_2.fq.gz_G_to_A.fastq, with the options: -q -N 1 --score-min L,0,-0.2 -p 32 --reorder --ignore-quals --no-mixed --no-discordant --maxins 1000 --quiet --nofw))
Found first alignment:
M00978:20:000000000-A5TV0:1:1101:15681:1332_1:N:0:3/1   77      *       0       0       *       *       0       0       TTATTGGTGGTTGGGTATTGTGGTTT   >AABBFFFDBBB222EEFGGGFFGHF       YT:Z:UP
M00978:20:000000000-A5TV0:1:1101:15681:1332_2:N:0:3/2   141     *       0       0       *       *       0       0       CCCACCACCACACCCAACTAATTTTT   >1A1ADDAFA1AEFGGGGGGFF3GFG       YT:Z:UP

[B]>>> Writing bisulfite mapping results to [COLOR="Red"]smallbample1_pe.sam[/COLOR] <<<[/B]

Current working directory is: /pipelines/trimmed

Now reading in and storing sequence information of the genome specified in: /custom_indexes/hg19canon_bsgenome/

chr chr1 (249250621 bp)
chr chr2 (243199373 bp)
chr chr3 (198022430 bp)
chr chr4 (191154276 bp)
chr chr5 (180915260 bp)
chr chr6 (171115067 bp)
chr chr7 (159138663 bp)
chr chr8 (146364022 bp)
chr chr9 (141213431 bp)
chr chr10 (135534747 bp)
chr chr11 (135006516 bp)
chr chr12 (133851895 bp)
chr chr13 (115169878 bp)
chr chr14 (107349540 bp)
chr chr15 (102531392 bp)
chr chr16 (90354753 bp)
chr chr17 (81195210 bp)
chr chr18 (78077248 bp)
chr chr19 (59128983 bp)
chr chr20 (63025520 bp)
chr chr21 (48129895 bp)
chr chr22 (51304566 bp)
chr chrX (155270560 bp)
chr chrY (59373566 bp)
chr chrM (16571 bp)


Reading in the sequence files ./smallsample1_R1_val_1.fq.gz and ./smallsample1_R2_val_2.fq.gz
Processed 47 sequences in total


Successfully deleted the temporary files smallsample1_R1_val_1.fq.gz_C_to_T.fastq and smallsample1_R2_val_2.fq.gz_G_to_A.fastq

Final Alignment report
======================
Sequence pairs analysed in total:       47
Number of paired-end alignments with a unique best hit: 28
Mapping efficiency:     59.6%

Sequence pairs with no alignments under any condition:  18
Sequence pairs did not map uniquely:    1
Sequence pairs which were discarded because genomic sequence could not be extracted:    0

Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT:       14      ((converted) top strand)
GA/CT/CT:       0       (complementary to (converted) top strand)
GA/CT/GA:       0       (complementary to (converted) bottom strand)
CT/GA/GA:       14      ((converted) bottom strand)

Number of alignments to (merely theoretical) complementary strands being rejected in total:     0

Final Cytosine Methylation Report
=================================
Total number of C's analysed:   326

Total methylated C's in CpG context:    24
Total methylated C's in CHG context:    71
Total methylated C's in CHH context:    231
Total methylated C's in Unknown context:        0

Total unmethylated C's in CpG context:  0
Total unmethylated C's in CHG context:  0
Total unmethylated C's in CHH context:  0
Total unmethylated C's in Unknown context:      0

C methylated in CpG context:    100.0%
C methylated in CHG context:    100.0%
C methylated in CHH context:    100.0%
Can't determine percentage of methylated Cs in unknown context (CN or CHN) if value was 0



====================
Bismark run complete
====================

Kind Regards,
Jens