[marcusmchale]

Thanks for the prompt reply, the manual for bismark suggests the following command:

bismark -n 1 -l 50 /data/genomes/homo_sapiens/GRCh37/ test_dataset.fastq

Which would call bowtie2 to use "-L 50".

Is there something I'm missing?

Oh, it's because of the differences in alignment strategy between bowtie1/bowtie2. Thanks for the lead!

[fkrueger]

If you type bowtie2 --help you can find the following text:

-L <int>           length of seed substrings; must be >3, <32 (22)

Obviously this is not mentioned in the manual but you need to keep the seed substrings in the range of 3 to 32. The default is 22. I hope this helps, Felix

[marcusmchale]

error with seedlen > 32

I get an error from bowtie2 when I try to define seed length of 50 in bismark. I haven't found any mention of this problem elsewhere nor mention of seedlen limits in the bowtie2 manual. Particularly, it seems strange given the recommended "typical' settings are for a seed length of 50 in the bismark manual. Can anyone help me to trace the source of this error?

Using Bowtie 2 index: /home/tair10/Bisulfite_Genome/CT_conversion/BS_CT

Error: -L argument must be <= 32; was50
Error: Encountered internal Bowtie 2 exception (#1)
Command: /cm/shared/apps/bowtie/2-2.1.0/bowtie2-align --wrapper basic-0 -q -N 1 -L 50 --score-min L,0,-0.2 --ignore-quals --norc -x /home/tair10/Bisulfite_Genome/CT_conversion/BS_CT -U 4AB_trimmed_r1.fastq_C_to_T.fastq
bowtie2-align exited with value 1

The alignment does seem to work when no seedlen is defined. Here is a sample of a read from the relevant fastq, you will note the read length is 50bp but I don't think this is relevant since the error says -L must be <= 32.

@HWI-D00458:73:C6EBDANXX:1:1101:1728:1972 1:N:0:GCTCTA_A
AGCGTGGTTTATTGATTTTTTAGATTTTCGGAATTTGAAGTTAGAGGTGT
+
CG>EFF<EFECE>D1<111@/FG>CFGGGG///0=1:FGGD1FE1FGEBG

P.S. This is my first comment in the forum (though I have been stalking this place for years) so I apologise if it is out of place.

[fkrueger]

I found one more bug in deduplicate_bismark. It is also in v0.14.4.
There are several calls of samtools directly use "samtools" instead of using $samtools_path. E.g. line 269, line 207.

Thanks for spotting that, I've fixed all these calls.

[lancelothk]

I found one more bug in deduplicate_bismark. It is also in v0.14.4.
There are several calls of samtools directly use "samtools" instead of using $samtools_path. E.g. line 269, line 207.

[fkrueger]

Bismark v0.14.4. New functionality and allele-specific alignment support

We have just released a new Bismark version (v0.14.4). This brings a few convenience features, adds some options and also fixes some bugs, further details are outlined below.

It is also worth mentioning that it should now be possible to use Bismark in conjunction with SNPsplit to align Bisulfite-Seq data in an allele-specific fashion against an N-masked genome if both genotypes are known. More information about this may be found on the SNPsplit project page.

o Bismark: Changed the FLAG values of paired-end alignments to the CTOT or CTOB strands so that reads can be properly displayed in SeqMonk when imported as BAM files. This change affects only paired-end alignments in --pbat or --non_directional mode. In detail we simply swapped the Read 1 and Read 2 FLAG values round so reads now resemble exactly concordant read pairs to the OT or OB strands. Note that results produced by the methylation extractor or further downstream of that are not affected by this change
o Bismark: Input files specified with filepath information for FastA files are now handled properly in --multicore runs (this was fixed only for FastQ files in the previous patch)
o Bismark: Unmapped and ambiguous files (options --unmapped and --ambiguous) are now written out as gzip compressed files by default
o Bismark: Changed the default mode of operation to --bowtie2. Bowtie (1) alignments may still be chosen using the option --bowtie1

o Bismark Genome Preparation: Changed the execution of the genome indexing of the parent process to system() rather than an exec() call since this seemed to lead to interesting faults when run in a pipeline setting
o Bismark Genome Preparation: Changed the default indexing mode to --bowtie2. Bowtie (1) indexing is still available via the option --bowtie1

o bismark2bedGraph: The coverage (.cov) and bedGraph (.bedGraph) files are now written out as gzip compressed files by default

o coverage2cytopsine: Added new option '--gc/--gc_context' to reprocess the genome to find methylation in GpC context. This might be useful for specialist applications where GpC methylases had been employed. The output format is exactly the same as for the normal cytosine report, and only positions covered by at least one read are reported (output file ends in .GpC_report.txt). In addition this will write out a Bismark coverage file (ending in GpC.cov)

o deduplicate_bismark: Removed redundant closing statements to get rid of warning messages
o deduplicate_bismark: The option --representative is no longer displayed in the help text. The option was once useful to determine the PCR bias that had been introduced by over digestion with bisulfite and is nearly always not what should be used for deduplication (it will be left in and is still functional for the time being though)

Bismark is available from the Babraham Bioinformatics project page.

[fkrueger]

Thanks for pointing out these issues. I have updated the manual and removed the superfluous closing statements. They will find their way into a new release which we'll be releasing later today.

Edit: Just as a quick word of warning: the --representative mode is almost certainly not what what you want to use because it will find the most highly amplified and thus biased sequence for a given position instead of a random. I will probably hide it from use in the next release...

[lancelothk]

Hi fkrueger,

I found two minor issues in bismark v0.14.3.
The deduplicate_bismark will give errors with --representative option:
Failed to close output filehandle: Bad file descriptor
Failed to close report filehandle: Bad file descriptor

I found out that it is caused by a bug in line:548. The } should be after two close lines, since OUT and REPORT have been closed in deduplicate_representative().

The -B/--basename <basename> option can be found in the script, but not in pdf version manual.

BTW, I finished extracting reports merging code into stand alone script. The most painful part is the global variable...

Thanks.

[lancelothk]

Thank you so much Felix. I will take a look at the source code.

[fkrueger]

Hi fkrueger,

I am currently trying to use bismark to analyse a huge BS-seq dataset in HPC environment. I am thinking to split the big fastq into smaller pieces, run bismark with each of them in one node, then merge back the BAM results. Do you have any suggestion how can I merge bismark reports? Do you have existing script to do this?

Thanks.

Hi Lanceloth,

As it stands there is no stand-alone script to merge the mapping reports, but the code should pretty much all be there for it is used for --multicore core runs anyway. The two subroutines

'merge_individual_mapping_reports' and
'read_alignment_report'

should contain everything. Let me know if you need help merging these into a stand-alone script.

Just as a short word of warning when you are trying to merge paired-end BAM files with samtools merge you need to make sure that the files are subsequently sorted by read name, otherwise the reads are not guaranteed to follow each other line by line. Would maybe the --multicore option be a little more feasible?

Best, Felix

[lancelothk]

Hi fkrueger,

I am currently trying to use bismark to analyse a huge BS-seq dataset in HPC environment. I am thinking to split the big fastq into smaller pieces, run bismark with each of them in one node, then merge back the BAM results. Do you have any suggestion how can I merge bismark reports? Do you have existing script to do this?

Thanks.

[jcrow06]

Hi fkrueger,

Honestly I had a feeling that this is an in-house issue.
When I approached them I was told "there is no annotation for chromosome etc and that it is per read based so I should summarize the read information to get results...."
From reading the well documented details of Bismark I did not feel this was the case

I wanted to double check before going back.
Thanks Fkrueger & dpryan for the help.
Much appreciated!
jcrow

[fkrueger]

Sorry about not getting that in the first place. I would tend to agree with Devon that someone has tampered with the chromosome names (the ones that appear as the first line in the .fa files used for the genome preparation >...). There is no place within Bismark that would affect or remove chr from a chromosome name. Maybe you could get back to the people who gave you the data and see if they could clarify this point?

[jcrow06]

Indeed as noted by dpryan
some of the hits do not have associated chr, looks like some are lost.
I looked through one of the sample files and found that there was ~ 6700 that had this missing Chr from ~6.2 Million.

sorry if my original post was not clear.
thanks for clearing it up dpryan

[dpryan]

Note the 7th, 15th and 16th lines, which lack a chromosome column.

Edit: Rather, they likely have incorrect chromosomes.