SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
CisGenome -- an integrated tool for ChIP-seq data analysis hji Bioinformatics 66 12-30-2014 02:55 PM
Bismark Bisulfite Aligner - Now supporting CpG, CHG and CHH context fkrueger Bioinformatics 27 10-11-2013 06:40 AM
Bismark v0.6.beta1: Now supporting gapped Bisulfite-Seq alignments fkrueger Bioinformatics 6 03-19-2012 06:06 AM
Apparent duplication levels incongruence between bismark and fastqc with BS-Seq data gcarbajosa Bioinformatics 2 12-13-2011 09:43 AM

Reply
 
Thread Tools
Old 09-04-2014, 01:38 PM   #421
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

Quote:
Originally Posted by gene_x View Post
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.
Quote:
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.
Quote:
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.
fkrueger is offline   Reply With Quote
Old 09-08-2014, 02:02 PM   #422
christophe
Junior Member
 
Location: Phoenix

Join Date: May 2010
Posts: 7
Default

Quote:
Originally Posted by adeirossi View Post
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
christophe is offline   Reply With Quote
Old 09-10-2014, 09:30 AM   #423
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default 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?
gene_x is offline   Reply With Quote
Old 09-11-2014, 02:41 AM   #424
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by AndrewMartins View Post
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
dpryan is offline   Reply With Quote
Old 09-11-2014, 04:46 AM   #425
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

Quote:
Originally Posted by gene_x View Post
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.
fkrueger is offline   Reply With Quote
Old 09-11-2014, 06:41 AM   #426
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default

Quote:
Originally Posted by fkrueger View Post
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
Code:
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!
gene_x is offline   Reply With Quote
Old 09-11-2014, 06:47 AM   #427
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

Indeed, you should not sort the files by coordinates at all before running the deduplication.
fkrueger is offline   Reply With Quote
Old 09-11-2014, 08:36 AM   #428
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default

Quote:
Originally Posted by fkrueger View Post
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?

Last edited by gene_x; 09-11-2014 at 08:43 AM.
gene_x is offline   Reply With Quote
Old 09-11-2014, 08:45 AM   #429
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

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.
fkrueger is offline   Reply With Quote
Old 09-11-2014, 08:51 AM   #430
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default

Quote:
Originally Posted by fkrueger View Post
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

Code:
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:

Quote:
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?
gene_x is offline   Reply With Quote
Old 09-12-2014, 06:02 AM   #431
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

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....
fkrueger is offline   Reply With Quote
Old 09-12-2014, 08:33 AM   #432
jnhutchinson
Junior Member
 
Location: Boston

Join Date: Feb 2012
Posts: 6
Default 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
jnhutchinson is offline   Reply With Quote
Old 09-12-2014, 12:30 PM   #433
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

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.
fkrueger is offline   Reply With Quote
Old 09-12-2014, 01:51 PM   #434
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default

Quote:
Originally Posted by fkrueger View Post
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?
gene_x is offline   Reply With Quote
Old 09-13-2014, 03:20 AM   #435
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default

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.
fkrueger is offline   Reply With Quote
Old 09-15-2014, 07:37 AM   #436
jnhutchinson
Junior Member
 
Location: Boston

Join Date: Feb 2012
Posts: 6
Default

Quote:
Originally Posted by fkrueger View Post
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.
jnhutchinson is offline   Reply With Quote
Old 09-15-2014, 08:38 AM   #437
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default

Quote:
Originally Posted by fkrueger View Post
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

Code:
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
Code:
samtools sort -n lane1.bam lane1_sort
Do step 1 and step 2 for all lanes.

3. Merge all sorted bam files

Code:
samtools merge all_lanes.bam lane1_sort.bam lane2_sort.bam ...
4. remove duplication

Code:
deduplicate_bismark -p --bam all_lanes.bam
5. extract methylation calls

Code:
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!
gene_x is offline   Reply With Quote
Old 09-15-2014, 11:09 AM   #438
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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.
dpryan is offline   Reply With Quote
Old 09-15-2014, 11:13 AM   #439
gene_x
Senior Member
 
Location: MO

Join Date: May 2010
Posts: 108
Default

Quote:
Originally Posted by dpryan View Post
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!
gene_x is offline   Reply With Quote
Old 10-01-2014, 05:42 AM   #440
fkrueger
Senior Member
 
Location: Cambridge, UK

Join Date: Sep 2009
Posts: 625
Default 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
fkrueger is offline   Reply With Quote
Reply

Tags
alignment, bisulfite, bisulphite, methylation, sequencing

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:18 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO