[pkachroo]

Bismark output details

Dear all,

Hello. My name is Priya. I am working on RRBS data. I used BISMARK for paired-end Illumina data. I had few queries regarding the ouput. Could someone kindly look into my output file data as I am not sure if the output is not unusual and I can proceed surely with further downstream analysis

1. I am attaching normal text document also containing mapping report for:

B1929_NoIndex_L008_R1_001_val_1.fq_Bismark_paired-end_mapping_report.txt

and few lines of CpG Methylation ratio and bedgraph file which do not seem properly aligned to me so I am a bit confused like mentioned by Felix Krueger above:
<Chromosome> <Start Position> <End Position> <Methylation Percentage>

2. Can methylation percentage be zero like I have in my text document ? In this case do I need to filter ratios with value < 5 or 10 < before proceeding further.

3. Is is better/compulsory to use --no_overlap option for paired end data because I did not use it for my results ?

4. Also I am not sure how to understand if C's are methylated by just looking at the Methylation ratio file as mentioned in the thread written by user "wilhelml":

http://seqanswers.com/forums/showthread.php?p=46285

read_id:2:82:3541:4552#0 - 1 1343509 z : C : un-methylated C on top strand
read_id:2:82:3541:4552#0 + 1 1343515 Z : C : methylated C on top strand

Is there any Key in BISMARK for understanding this ?

I would be grateful if I get some help from this forum for my further analysis.

Thanks & sincere regards,

Priya

[kuul2jai]

Hi,
I tried to used the option to align multiple fastq.gz files at once.

bismark $refgenome_dir -1 lane1_r1.fastq.gz,lane2_r1.fastq.gz -2 lane1_r2.fastq.gz,lane2_r2.fastq.gz

Only lane1_r1 & r2 were recognized, but not lane2 r1 & r2. For now I think I can fix this by catting the fastq files together before running bismark.

PS: I tried different seperators: space, comma, comma with space

[fkrueger]

Hi,
I tried to used the option to align multiple fastq.gz files at once.

bismark $refgenome_dir -1 lane1_r1.fastq.gz,lane2_r1.fastq.gz -2 lane1_r2.fastq.gz,lane2_r2.fastq.gz

Only lane1_r1 & r2 were recognized, but not lane2 r1 & r2. For now I think I can fix this by catting the fastq files together before running bismark.

PS: I tried different seperators: space, comma, comma with space

Paired-end mode requires files to be comma separated, so
-1 lane1_r1.fastq.gz,lane2_r1.fastq.gz -2 lane1_r2.fastq.gz,lane2_r2.fastq.gz
should have worked. Could you send me the on screen messages that you got to [email protected]

[fkrueger]

Dear all,

Hello. My name is Priya. I am working on RRBS data. I used BISMARK for paired-end Illumina data. I had few queries regarding the ouput. Could someone kindly look into my output file data as I am not sure if the output is not unusual and I can proceed surely with further downstream analysis

1. I am attaching normal text document also containing mapping report for:

B1929_NoIndex_L008_R1_001_val_1.fq_Bismark_paired-end_mapping_report.txt

and few lines of CpG Methylation ratio and bedgraph file which do not seem properly aligned to me so I am a bit confused like mentioned by Felix Krueger above:
<Chromosome> <Start Position> <End Position> <Methylation Percentage>

2. Can methylation percentage be zero like I have in my text document ? In this case do I need to filter ratios with value < 5 or 10 < before proceeding further.

3. Is is better/compulsory to use --no_overlap option for paired end data because I did not use it for my results ?

4. Also I am not sure how to understand if C's are methylated by just looking at the Methylation ratio file as mentioned in the thread written by user "wilhelml":

http://seqanswers.com/forums/showthread.php?p=46285

read_id:2:82:3541:4552#0 - 1 1343509 z : C : un-methylated C on top strand
read_id:2:82:3541:4552#0 + 1 1343515 Z : C : methylated C on top strand

Is there any Key in BISMARK for understanding this ?

I would be grateful if I get some help from this forum for my further analysis.

Thanks & sincere regards,

Priya

Hi Priya,

To 1) What exactly doesn’t seem properly aligned to you in the output files? Just as a reminder, bedGraph files are naturally 0-based, which might cause calls to look offset by 1. What is a bit puzzling is that the bedGraph file doesn’t appear to be sorted properly by chromosome and position, did you get any warning (memory etc) during the run?
As a general comment your files seem to have aligned with an alright mapping efficiency (looks a bit like ~40bp long reads?), however the methylation percentages especially in non-CG context are really quite high:
C methylated in CpG context: 73.5%
C methylated in CHG context: 29.9%
C methylated in CHH context: 29.7%

This normally indicates that something did not work out as expected with the bisulfite conversion efficiency. If this is the case then you should probably not start analyzing the results much further…

To 2) Yes, the methylation can be zero whenever you only have unmethylated but no methylated reads. You do not want to filter out on reads based on the percentage methylation, but you might want to filter position based on the read coverage. If you run the methylation extractor with the options --bedGraph and --counts you can see the exact number of methylated and unmethylated calls as well, and you could also use --cutoff_threshold to filter for a certain read depth per position. You could of course also filter later with the package you are using for the downstream analysis (we use the Bisfulfite methylation pipeline in SeqMonk for such a coverage filtering step).

To 3) We would recommend using --no_overlap for paired end data because you would get an unfair coverage of overlapping parts of reads, which might affect coverage filtering of reads. As overlapping calls are technical replicates of each other that were merely copied by PCR they should be the same anyway.

To 4) This is really simple. Cytosine positions receiving a '+' are methylated and '–' means unmethylated. The last column only indicates the cytosine context which may seem redundant but if you write several outputs into the same file (e.g. CHH and CHG context) it is useful if you want to separate them out later. So,
read_id:2:82:3541:4552#0 - 1 1343509 z is unmethylated (CpG)
read_id:2:82:3541:4552#0 + 1 1343515 Z is methylated (CpG)

It is not immediately obvious just from this position to tell the strand on which this C resides, however the methylation extractor output sorts the calls into files that carry the strand name in their file names.

If you have any further questions don’t hesitate to contact me further via email.

Best,
Felix

[pkachroo]

Thanks a lot Felix for such elaborate & clear reply. It surely helped me to clear these initial doubts and i can proceed accordingly

PS: Regarding alignment i was confused before but now it is clear to me & there was no memory warning after running the cluster job for Bismark.

[mxenoph]

Low mapping efficiency and bedGraph error

Hi,

I just started using bismark (v0.7.7) for a set of publicly available dataset (GSE30199) from mouse and I am running on some errors and low mapping efficiency for paired-end reads.

The dataset contains both single-end and paired-end reads. Bismark maps the SE reads with an efficiency of about 55% for most of the samples, but for the PE data the mapping efficiency for all the samples is really low (less than 10% nearly in all samples). I initially ran bismark with the -q -n 1 -k 2 --best and the -I and -X as defined by the PE library and --directional, with the bowtie 1 aligner; however the mapping efficiency only increased by ~5% when I ran it again changing only the options -n to 2 and -I and -X set to the default values. A colleague of mine ran an older version of bismark on the same data, but treated the PE files as SE and got mapping efficiencies of ~60%. Also, the bisulfite genome was prepared for MM9 in that case whereas I am using MM10 (but don't think that caused the poor mapping).

I am wondering if treating PE as SE would be a good () solution to my problem or if I could use a different set of options to increase the mapping efficiency for the PE reads.

Now, I also get an error (see below) when running bismark_methylation_extractor on the SE sam files.

No such file or directory at /homes/mxenoph/bin/bismark_methylation_extractor line 175, <IN> line 35808595[/FONT]

the script has already processed 35808571 lines of the sam file and produced the output files except the bedGraph one. I briefly checked the bismark_methylation_extractor and on line 175 it tries to create the bedGraph file. What I don't get is why it tries to set the output directory for that file as the directory of the script and not the specified (by -o) one.

Any input would be very much appreciated. Thanks for bearing with me and my long post!

[fkrueger]

Hi,

I just started using bismark (v0.7.7) for a set of publicly available dataset (GSE30199) from mouse and I am running on some errors and low mapping efficiency for paired-end reads.

The dataset contains both single-end and paired-end reads. Bismark maps the SE reads with an efficiency of about 55% for most of the samples, but for the PE data the mapping efficiency for all the samples is really low (less than 10% nearly in all samples). I initially ran bismark with the -q -n 1 -k 2 --best and the -I and -X as defined by the PE library and --directional, with the bowtie 1 aligner; however the mapping efficiency only increased by ~5% when I ran it again changing only the options -n to 2 and -I and -X set to the default values. A colleague of mine ran an older version of bismark on the same data, but treated the PE files as SE and got mapping efficiencies of ~60%. Also, the bisulfite genome was prepared for MM9 in that case whereas I am using MM10 (but don't think that caused the poor mapping).

I am wondering if treating PE as SE would be a good () solution to my problem or if I could use a different set of options to increase the mapping efficiency for the PE reads.

Now, I also get an error (see below) when running bismark_methylation_extractor on the SE sam files.

No such file or directory at /homes/mxenoph/bin/bismark_methylation_extractor line 175, <IN> line 35808595[/FONT]

the script has already processed 35808571 lines of the sam file and produced the output files except the bedGraph one. I briefly checked the bismark_methylation_extractor and on line 175 it tries to create the bedGraph file. What I don't get is why it tries to set the output directory for that file as the directory of the script and not the specified (by -o) one.

Any input would be very much appreciated. Thanks for bearing with me and my long post!

Hi there,

If you see a much lower mapping efficiency for paired end files compared to single end mode it is usually caused by one of the following factors:
- Adapter contamination or poor qualities, possibly especially in read 2
- The fragment size was larger than the allowed insert size (which is obviously not the problem here)
- The fragment length was so short that you started getting paired end reads that completely overlap

To rectify that I would suggest that you get Trim Galore and simply run it with the default parameters, like so:

trim_galore --paired --trim1 read_1.fastq read_2.fastq

This will
- remove poor quality base calls from the 3’ end of reads
- remove Illumina adapters from both reads
- trim 1 extra bp from the 3’ end so that completely overlapping fragments do still align with Bowtie 1

The Trim Galore documentation gives you more details about why these overlapping reads are a problem, and you can also take a look at this protocol of how we treat our samples internally. I am convinced that following these steps will fix the low mapping efficiencies you are seeing.

In general, using PE as SE reads would work, but it is obviously not making the most out of the data. In any case, for paired-end data you should use the option --no_overlap for the methylation extractor to disregard redundant information of overlapping reads.

About the error message in the methylation extractor: I would probably need some more information about the command you used, line 176 in my current development version (0.7.8) reads:
open (OUT,'>',$output_dir.$bedGraph_output) …
which should attempt to write the bedGraph output to the directory specified with –o. Please find the current version of the methylation extractor attached, and just send me and email if there are any further problems.

[mxenoph]

Thank you for your quick reply.

Based on the fastqc reports there wasn't any adapter contamination problem. However the PHRED scores were quite low towards the read end..didn't think that would affect the mapping so much though and was counting on -n set to a high number to offset any issues.
Is there an empirical rule that tells you that you should trim the data if a certain percentage of the bases in the read have low PHRED or should I always perform the trimming? I currently run bismark on the trimmed reads and hopefully that will increase the mapping efficiency.

About the error, I called the command from a bash script:

bismark_methylation_extractor -o "$out_dir/$exp/" --comprehensive --report --bedGraph -s "$infile"

and I am now running the same command from the terminal with absolute and relative path for the output directory, just to check if that's what causes the problem. Hopefully I am not doing more stupid than that!

[fkrueger]

I would always run adapter trimming, simply because it removes incorrect portions of the read that you can't spot between millions of other reads and it also allows you to use fairly strict mapping settings (such as the default parameters), which is a) much quicker and b) less prone to mismapping events. I believe that increasing -n from 1 to 2 for the test dataset in the Bismark paper would increase the alignment time from ~30 mins to 150 mins (!) without noticably improving the mapping results....

Can you please keep me posted whether the methylation extractor works with your two approaches and/or if it was the evil bash script...

Cheers,
Felix

[fkrueger]

We have just released a new version of Bismark (v0.7.8) which adds a few more functions to Bismark and the methylation extractor. Here are the main changes:

• Bismark: Added new option '--non_bs_mm' which prints an extra column at the end of SAM files showing the number of non-nisulfite mismatches of a read. This option is not available in '--vanilla' format.
Format for single-end reads: "XA:Z:mismatches".
Format for paired-end reads: read 1: "XA:Z:mismatches", read 2: "XB:Z:mismatches"
• Bismark: The mapping report file names were changed to _bismark_(SE/PE)_report.txt (Bowtie 1) or _bt2_bismark_(SE/PE)_report.txt (Bowtie 2) to keep it more uniform
• Methylation extractor: The input file(s) may now be specified with a file path which abrogates the need to be in the same directory as the input file(s) when calling the methylation extractor
• Methylation extractor: Added new function '--buffer_size' to increase the physical memory used for the sorting the output by chromosomal positions (only needed for bedGraph output)
• Methylation extractor: Reference sequence files containing pipe ('|') characters were found to crash the methylation extractor as the chromosome name was used for filenames. These characters are now replaced with underscores when the reads are sorted during the bedGraph step
• Updated the Bismark User Guide with sections for the bedGraph and genome-wide methylation report outputs, and Appendix IV is now showing alignment stats for the test data

Bismark can be downloaded here: http://www.bioinformatics.babraham.a...jects/bismark/

[fkrueger]

We have just released a new version of Bismark (v0.7.9) which fixes the following three issues for the methylation extractor:

- The new function '--buffer_size <string>' for the bedGraph sort command is set to the new default value of 2G (sort would die if this option was not set)
- The replacement of pipe ('|') characters in the name of reference chromosomes for the bedGraph sorting step should now work as intended
- If multiple files are to be processed for genome-wide methylation reports in a single command, the reference genome is only read in once (instead of stopping due to the presence of several chromosomes with the same name...)

Bismark is available here.

[ruggedtextile]

I'm running Bismark on some PE RRBS (non-directional) data that we've just recieved. This our first time looking at RRBS data and running Bismark and I'm surprised at the speed.

Running with bowtie2 specified and default parameters otherwise:

bismark --bowtie2 --non_directional --phred64-quals /mnt/nas_omics/shared_data/iGenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/WholeGenomeFasta -1 4_lib1_lane1_1.fq.gz -2 4_lib2_lane1_2.fq.gz

And it's taken 3 hours to process 5M sequences. I.e this point in the output has been reached:

...
Processed 4600000 sequences so far
Processed 4700000 sequences so far
Processed 4800000 sequences so far
Processed 4900000 sequences so far
Processed 5000000 sequences so far
Processed 5100000 sequences so far
Processed 5200000 sequences so far
...

This machine has 96GB RAM and 8 CPU cores, so I feel it should be faster than this? Strangely when I run top I see 4 bowtie2-align processes (as expected), but with two of them on 100% and the other two on ~20%. Is it normal for two of the alignment processes to be throttled in this way or is this indicative of a bottleneck somewhere else (IO maybe?).

Other tools like BSMap were orders of magnitude faster than this, so I feel I must be missing something or have something misconfigured somewhere. Any advice welcome.

[fkrueger]

This speed seems very slow indeed, I'll give a few ideas (in no particular order) what I could think of:

- has the data been trimmed to remove adapter, poor quality and the RRBS bias? I would say a typical directional RRBS run with trimmed reads with Bowtie1 should align ~15-20 million reads/h, for non-directional mode probably a bit less. I don't know any numbers for non-directional Bowtie2 RRBS though.

- The use of Bowtie 2 in Bismark is in my experience indeed quite a bit slower than Bowtie 1, especially for shorter reads. How long is your read length, and do you really need to run alignments with Bowtie2?

- Is the data really non-directional? There are not many protocols that generate this type of data, and this can indeed slow down the alignments noticably. That is because if the data is not non-directional then the 2 threads running against the complementary strands are trying very hard but can't find many alignments. This might also explain why 2 threads are working very hard, while the other 2 are already idling as they are waiting for the other threads to catch up.

- is the data really encoded with the old Illumina qualities (phred64). Not sure if Bowtie2 would complain about it but the wrong quality encoding scheme could make alignments much slower.

If you could send me a FastQC report (the zip file) of your FastQ files via email I could take a quick look and advise further.

[ruggedtextile]

This speed seems very slow indeed, I'll give a few ideas (in no particular order) what I could think of:

- has the data been trimmed to remove adapter, poor quality and the RRBS bias? I would say a typical directional RRBS run with trimmed reads with Bowtie1 should align ~15-20 million reads/h, for non-directional mode probably a bit less. I don't know any numbers for non-directional Bowtie2 RRBS though.

- The use of Bowtie 2 in Bismark is in my experience indeed quite a bit slower than Bowtie 1, especially for shorter reads. How long is your read length, and do you really need to run alignments with Bowtie2?

- Is the data really non-directional? There are not many protocols that generate this type of data, and this can indeed slow down the alignments noticably. That is because if the data is not non-directional then the 2 threads running against the complementary strands are trying very hard but can't find many alignments. This might also explain why 2 threads are working very hard, while the other 2 are already idling as they are waiting for the other threads to catch up.

- is the data really encoded with the old Illumina qualities (phred64). Not sure if Bowtie2 would complain about it but the wrong quality encoding scheme could make alignments much slower.

If you could send me a FastQC report (the zip file) of your FastQ files via email I could take a quick look and advise further.

Thanks for the fast reply Felix. You've given me a few things to try straightaway. In answer to your questions:

1. Yes the data was trimmed, but by the sequencing center not by me, so I don't know the precise details. I will circle back and check with them, but the data seems to align fine with BSMap so I suspect this is not the problem.

2. These are 50bp PE reads. I have no reason to use bowtie2 over bowtie1. I had naively assumed that 2 was always better than 1, but clearly not in this case. I will rerun with bowtie1 and see if that resolves anything.

3. Hmmm. I will recheck about the directionality.

4. Again, I will recheck with the sequencing center. FastQC reports the encoding as 'Illumina 1.5', which I thought was Phred64.

I'll email you the FastQC report now. It all looked fine to me, but this is my first time with any form of BS data.

[oria34]

HI all,


This is my first time using Bismark (actually, my first time in this NextGen stuff) and I really think is a very complete software that makes things pretty easy to the newbies. So, many thanks for you work Felix

I have some Fish BS converted DNA set to compare methylation state at some genes of interest so the BedGraph option is a really useful tools to my purposes. However I am having some problems when the script tries to write down the file.

My data is: paired-en, directional (I guess, not sure), It was trimmed and the encoding is Illumina 1.5 (Phred 64???). Roughly 50 million reads.



Warning: I am really new with all this stuff (coming from classical Pop. Genetics) so please people, be gentle with bioinformatic terms

Firstly I tried to run the scripts in two different files (in two different remote servers) but apparently the temporal files that the script creates have the same name so the program failed. Thas was pretty easy to fix by running the scripts in two different folders, one for each data set so the temporal files with the same name are now in different folder.

But...I will crash again but for different reasons:


Processed 32126386 lines from BPF1_methyl_1.fq_bismark_bt2_pe.sam in total
Total number of methylation call strings processed: 64252772

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

Total methylated C's in CpG context: 121971031
Total methylated C's in CHG context: 869279
Total methylated C's in CHH context: 2305201

Total C to T conversions in CpG context: 71035093
Total C to T conversions in CHG context: 301484428
Total C to T conversions in CHH context: 799786286

C methylated in CpG context: 63.2%
C methylated in CHG context: 0.3%
C methylated in CHH context: 0.3%


================================================================
Methylation information will now be written into a bedGraph file
================================================================

Using the following files as Input:
CpG_context_BPF1_methyl_1.fq_bismark_bt2_pe.txt

Now writing methylation information for file CpG_context_BPF1_methyl_1.fq_bismark_bt2_pe.txt to individual files for each chromosome
Failed to open filehandle: Too many open files at bismark_methylation_extractor line 2921, <IN> line 65658042.
32584.414u 1126.687s 9:56:48.66 94.1% 0+0k 0+164137336io 0pf+0w

Any ideas? I will try the separate script from the Bismark site in the meanwhile.

For the record: Alignment took something like 20h and methylation_extractor around 12h before crash. I used Botwie 2 - I also naively assumed that was better to use a newer version.

Here is the aligment report:

Bowtie was run against the bisulfite genome of /fs/lustre/wrk/user/Bisulfite/Bismark_Runs/ with the specified options: -q --phred64 --score-min L,0,-0.2 --ignore-quals --no-mixed --no-discordant --maxins 500

Option '--directional' specified: alignments to complementary strands will be ignored (i.e. not performed)
Final Alignment report
======================
Sequence pairs analysed in total: 62785587
Number of paired-end alignments with a unique best hit: 40067266
Mapping efficiency: 63.8%
Sequence pairs with no alignments under any condition: 18464624
Sequence pairs did not map uniquely: 4253697
Sequence pairs which were discarded because genomic sequence could not be extracted: 281

Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT: 20027636 ((converted) top strand)
GA/CT/CT: 0 (complementary to (converted) top strand)
GA/CT/GA: 0 (complementary to (converted) bottom strand)
CT/GA/GA: 20039349 ((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: 1631659739

Total methylated C's in CpG context: 131857571
Total methylated C's in CHG context: 1142574
Total methylated C's in CHH context: 3268321

Total C to T conversions in CpG context: 80295548
Total C to T conversions in CHG context: 361837285
Total C to T conversions in CHH context: 1053258440

C methylated in CpG context: 62.2%
C methylated in CHG context: 0.3%
C methylated in CHH context: 0.3%