[fkrueger]

You mean "This should affect all C contexts and not only CHH?
How do you click the CHH and CHG away? I can only zoom in certain bases...

If you click on e.g. 'CHH total calls' in the legend that track will disappear/show again and the scale adjust accordingly.

Can you help me clarify one thing? In the M-bias plot for read2, is position 1 the start of the read2? I'm a bit confused about the orientation of the read2. I thought in M-bias plot, position 1 is the the end of read2 (after read 2 is reverse complemented)? probably my understanding was wrong...

The first position does indeed correspond to the first base or read 2 as it is read by the sequencer. If this was the reverse complemented end of the sequenced fragment most of the artefacts would be completely obscured because these would likely happen at various positions within the read.

Also, in M-bias plot for read2, the first 2 bases have much lower average CpG methylation, that suggests they should be ignored in the downstream analysis? I could use the -ignore_r2 option in methylation extractor to do it?
Thanks!

The first positions of read 2 are the result of carrying out the end-repair/A-tailing procedure with unmethylated cytosines. Sonication is quite likely to leave lagging ends, so these would be filled in with unmethylated Cs (of course only at cytosine positions) which would subsequently be converted to Ts in the bisulfite reaction. If you then sequence Read 2, which is the reverse complement of such a T, it would appear as A and be called as unmethylated. As you rightly said this can be remedied by using --ignore_r2 2.

[christophe]

Hi Christophe,

I have also seen a falling number of methylation calls for read 2 in my data. I have always assumed this was because I ran bismark_methylation_extractor with the --no_overlap option and many of the read pairs had overlaps (i.e. insert length less than twice the read length). With --no_overlap, when the extractor encounters an overlap between read 1 and read 2, I believe it uses the read 1 calls only and discards the read 2 calls, which would explain why the number of read 2 calls decreases as you go toward its 3' end. Because this "no overlap" filtering occurs after mapping, I assume there is no impact on the alignment of read 2.

If, on the other hand, you are not using the --no_overlap option or you don't have many reads pairs with inserts <200 bp, then I have no idea what's going on.

Andrew

Hi Andrew and Felix,

Sorry for this late message.

Thanks Andrew for your reply to my question regarding M-bias.
I effectively ran the Extractor with the --no_overlap option on before. As suggested I just rerun Bismark Methyl Extractor without the --no_overlap option and it worked. I do not see the number of calls falling anymore for CHH and CHG.

Thanks Felix for confirming Andrew's answer. Yes, you are right, the CpG calls was also falling slightly but better now.

I also noticed that the R2 calls counts are more heterogeneous across the read length than the number of calls across R1. I do not know if anyone else noticed that.

Thanks,
Best,
Christophe

[gene_x]

sorted bam

Hi, Felix,
I have another question about the alignment files. I just checked the bam file after running bismark alignment and it is unsorted. Then using the bam output to do deduplication, the resulting deduplicated.bam is also unsorted. My first question is deduplication step doesn't require the input bam to be sorted?

I was trying to merge several deduplicated.bam files from different lanes into one big final bam and run methylation extractor on it and the resulting coverage file (.zero.cov with --zero_based specified) has methylation calls on non-Cs locations like chr1 133 134. That's how I traced it back to the issue of unsorted bam..

So it's probably better to implement sorting into the bismark pipeline?

[dpryan]

Your post is about a tool that can be used for to would be completely obscured because these would likely happen at various positions within the read. i can told you for use another site of tool its mainly supply pipeline on customers.

Gotta love incoherent spam

[fkrueger]

Hi, Felix,
I have another question about the alignment files. I just checked the bam file after running bismark alignment and it is unsorted. Then using the bam output to do deduplication, the resulting deduplicated.bam is also unsorted. My first question is deduplication step doesn't require the input bam to be sorted?

I was trying to merge several deduplicated.bam files from different lanes into one big final bam and run methylation extractor on it and the resulting coverage file (.zero.cov with --zero_based specified) has methylation calls on non-Cs locations like chr1 133 134. That's how I traced it back to the issue of unsorted bam..

So it's probably better to implement sorting into the bismark pipeline?

It should be generally possible to merge several bam files before running the extractor without a problem, and files do not have to be sorted for this. You only need to ensure that read 1 and read 2 always follow each other in the bam file, else you would run into problems with stand identity and the no_overlap function. Some merging tools do not guarantee this order, in which case it might help to sort the reads by read id to team pairs up again. I'm currently at a conference so I'm somewhat slow replying.

[gene_x]

It should be generally possible to merge several bam files before running the extractor without a problem, and files do not have to be sorted for this. You only need to ensure that read 1 and read 2 always follow each other in the bam file, else you would run into problems with stand identity and the no_overlap function. Some merging tools do not guarantee this order, in which case it might help to sort the reads by read id to team pairs up again. I'm currently at a conference so I'm somewhat slow replying.

I indeed run into problem after sorting bam files based on chromosomal coordinates.. the methylation extractor report an error saying that read1 and read2 are not the same and suggested me to use an unsorted bam..

I then tried to sort the deduplicated bam with the command

samtools sort -n deduplicated.bam deduplicated_sort

And then run methylation extractor on the deduplicated_sort.bam, however, the same error message came up again and methylation extractor couldn't proceed.. Is it because after deduplication, read1 and read2 within one bam file don't always match even after I sort them based on read names? Should I not sort anything from the beginning at all?

Thanks!

[fkrueger]

Indeed, you should not sort the files by coordinates at all before running the deduplication.

[gene_x]

Indeed, you should not sort the files by coordinates at all before running the deduplication.

Should I sort it by read names using -n in samtools sort just like I listed before doing deduplication? Or should I not sort it at all?

[fkrueger]

Of you can, just use the Bismark files as they are generated. If you have to
using samtools sort -n should do the trick as well.

[gene_x]

Of you can, just use the Bismark files as they are generated. If you have to
using samtools sort -n should do the trick as well.

OK. I'm running that now.

Another question, I'm not sure if my previous run without doing any sorting was correct.

Here is my command

samtools merge input.bam plate1/plate1_all_sort.bam plate2/plate2_all_sort.bam
deduplicate_bismark -p --bam input.bam
bismark_methylation_extractor -p --no_overlap --ignore 3 --ignore_r2 3 --bedGraph --counts --zero_based --report input.deduplicated.bam 2> input.meth_extractor_log.txt

The resulting input.deduplicated.zero.cov file have methylation coverage on bases that are not Cs. Here are the top rows:

chr2 133 134 0 0 8
chr2 134 135 0 0 1
chr2 228 229 0 0 9
chr2 229 230 0 0 2
chr2 262 263 0 0 15
chr2 263 264 0 0 2
chr2 304 305 0 0 13
chr2 305 306 0 0 1
chr2 316 317 0 0 11
chr2 317 318 0 0 1
chr2 318 319 0 0 12
chr2 319 320 0 0 1
chr2 326 327 0 0 11

Do you have any suggestions what could go wrong?

[fkrueger]

Hi gene_x,

not quite sure what is going here or which genome you are aligning your reads to, but for at least human or mouse chromosome 2 the first couple of megabases are masked by Ns, and there is no way that Bismark would map any reads to these sequences or extract methylation calls from it....

[jnhutchinson]

Using fill-in position to determine bisulfite conversion efficiency?

Has anyone here actually tried this?

I'm seeing some strange RRBS results with high non-CpG methylation levels (~4%) in my samples. I suspect it's due to poor(er) conversion, but my client thinks it's due to the non-standard areas we selected for RRBS (larger fragments targetted to non-CpG island areas). I could use a more objective measure of conversion to help decide the issue.

Would appreciate any and all pointers/scripts to get this to work.

Best,
John

[fkrueger]

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.

[gene_x]

Hi gene_x,

not quite sure what is going here or which genome you are aligning your reads to, but for at least human or mouse chromosome 2 the first couple of megabases are masked by Ns, and there is no way that Bismark would map any reads to these sequences or extract methylation calls from it....

Right, there should not be methylation calls on those regions.

A related question, if you have multiple lanes of data (fastq), do you run alignment on each individual lane first and them merge the resulting bam file or do you merge all fastq first and then do a big alignment on the merged fast file?

[fkrueger]

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.