N.B.: There is now a tutorial available here. This tutorial largely replaces this README file and users are encouraged to read it.
If you use Bison in your work please site the following: Ryan D.P. and Ehninger D. Bison: bisulfite alignment on nodes of a cluster. BMC Bioinformatics 2014, Oct 18;15(1):337
One can index all fasta files (files with extension .fa or .fasta) in a directory as follows:
bison_index [OPTIONS] directory/Options that are not specific to bison are simply passed to bowtie2, which must
be in your PATH. The output is placed under directory/bisulfite_genome.
Alignment can be performed as follows (bison_herd is the same):
mpiexec bison [OPTIONS] -g directory/ {-1 fastq_1.gz -2 fastq_2.gz | -U fastq.fq}"directory" is identical to that used for indexing. For further details type "bison -h". For non-directional libraries, "mpiexec -N 5" should be used, otherwise "mpiexec -N 3". Resource managers, such as slurm, should work in an equivalent manner. All options not explicitly mentioned by typing "bison -h" are passed to bowtie2. Consequently, using the --very-sensitive or --dovetail options will work as expected. Bison already passes the following flags to bowtie2:
-q --reorder --no-mixed --no-discordantbison_herd is equivalent, except that you can specify more nodes. You may also
input multiple files (comma-separated, no spaces) to align, in which case
alignments will be printed to multiples files. Furthermore, you may use
wild-cards in your file list. For example:
mpiexec -N 17 bison_herd -o Alignments -g directory/ -1 exp1/sample*_1.fq.gz,/some/other/path/foo*_1.fq.gz -2 exp1/sample*_2.fq.gz,/some/other/path/foo*_2.fq.gzMake sure to not have multiple input files with the same name
(e.g., sample*/read1.fastq), as they will all be written to the same file
(overwriting any subsequent alignments)!
In non-targeted sequencing experiments, it is usually wise to mark likely PCR
duplicates. These are then ignored by the methylation extractor so as to not
increase the error rate of any particular position. bison_markduplicates and
able to process a BAM file generated by bison/bison_herd and produce an
identical BAM file with the 0x400 bit set in the FLAG field for reads that are
likely duplicates. This step is not required and should be avoided if you are
performing RRBS or other targeted sequencing, as the false-positive rate of will
then be too high.
There is also a methylation extractor that produces a bedGraph file, called
bison_methylation_extractor. Note, coordinate-sorted BAM files should not
be used! The methylation extractor can be told to ignore certain parts of each
read. This is particularly useful in cases where there is methylation bias
across the length of reads (i.e., if one plots the average methylation
percentage summed per position over all reads, the value goes up/down toward the
5' or 3' end). It is recommended to always run bison_mbias (with the -pdf option
if you have R and ggplot2 installed) to generate the required information for
constructing an M-bias plot. The bison_mbias2pdf script can convert this to a
PDF file (or a series of PNG files) and will also suggest what, if any, regions
should be ignored. These regions are strand and read number (in the case of
paired-end reads) dependent. While the suggested regions are often good, the
should not be blindly accepted (just look at the graph and use your best
judgement).
See the "Auxiliary files" section, below, for additional files.
The following programs and scripts will be available if you type "make auxiliary":
This python script can accept a filename prefix and the names of at least 2 bedGraph files and output 3 files for input into BSseq. A single chromosome can be processed at a time, if desired, by using the -chr option. The output files will be named $prefix.M, $prefix.Cov, and $prefix.gr. $prefix.M is a matrix with a header line that lists the number of reads supporting methylation at each site in the bedGraph files. If there is no coverage in a given sample, the value is set to 0. $prefix.Cov is the analogous file listing coverage in each sample (again, 0 denotes no coverage). $prefix.gr lists the coordinates for each line in the .Cov and .M files. Loading these files into R would be performed as follows (in this example "Chr17" was the prefix):
M <- as.matrix(read.delim("Chr17.M", header=T))
Cov <- as.matrix(read.delim("Chr17.Cov", header=T))
bed <- read.delim("Chr17.bed", header=F)
#Remember that BED and bedGraph files are 0-based!
gr <- GRanges(seqnames=Rle(bed$V1),ranges=IRanges(start=bed$V2+1, end=bed$V3), strand=Rle("*", nrow(bed)))
groups <- data.frame(row.names=colnames(M),
var1 <- c(1,1,1,1,2,2,2,2)) #A very simple experiment with 2 groups of 4 samples
BS1 <- BSseq(M=M, Cov=Cov, gr=gr, pData=groups, sampleNames=colnames(M)) #You'll want to set some of the additional options!bedGraph2methylKitAs above, but each bedGraph file is converted to a .methylKit file. The bedGraphs should be of CpGs and not have had the strands merged (i.e., don't run the merge_CpGs command below).
bedGraph2MOABSLike bedGraph2methylKit, but each bedGraph file is converted to a .moabs file.
The bedGraph files should ideally contain single-C metrics rather than having
been merged to form CpG metrics, though both are supported. The resulting .moabs
files can then be used by mcomp in the MOABS package.
bedGraph2MethylSeekRAs above, but each bedGraph file is converted into a .MethylSeekR file. The
bedGraphs MUST be merged before-hand with bison_merge_CpGs to create per-CpG
metrics, as this is what MethylSeekR is expecting. Input is performed with the
readMethylome() function. Chromosome lengths can be computed with:
samtools faidx genome.fa
where genome.fa is a multifasta file containing all of the chromosomes. The
resulting .fai file is simply a text file and can be loaded into R with:
fai <- read.delim("genome.fa.fai", header=F)
chromosome_lengths <- fai$V2
names(chromosome_lengths) <- fai$V1
d <- readMethylome("file.MethylSeekR", chromosome_lengths)make_reduced_genomeCreate a reduced representation genome appropriate for reads of a given size ($size, default is 36bp). MspI and TaqI libraries are supported. Nucleotides greater than $size+10% are converted to N.
merge_bedGraphs.pyThis will merge bedGraphs from technical replicates of a single sample into a single bedGraph file, summing the methylation metrics as it goes. The output, like the input is coordinate sorted.
bison_merge_CpGsMethylation is usually symmetric at CpG sites. While the output bedGraph files have a single-C resolution, this will convert that to single-CpG resolution by summing Cs in the same CpG from opposite strands. This saves space and will often speed up downstream statistics.
While there are helper scripts, mentioned above, for a number of packages, other packages either do not require a helper script or can use one of the aforementioned scripts. Import instructions for such packages are mentioned below.
BiSeq requires input in an identical format as BSseq. Consequently, just use the bedGraph2BSseq.py helper script. The following example commands should then suffice to load everything into R:
exptData <- SimpleList(Sequencer="Some sequencer", Year="2014") #This is just descriptive information
M <- as.matrix(read.delim("chr17.M", header=T))
Cov <- as.matrix(read.delim("chr17.Cov", header=T))
bed <- read.delim("chr17.bed", header=F)
gr <- GRanges(seqnames=Rle(bed$V1),ranges=IRanges(start=bed$V2+1, end=bed$V3), strand=Rle("*", nrow(bed)))
groups <- DataFrame(row.names=colnames(M),
group = c(1,1,1,1,2,2,2,2)) #A very simple experiment with 2 groups of 4 samples
d <- BSraw(exptData=exptData, rowData=gr, colData=groups, totalReads=Cov, methReads=M)The BEAT Bioconductor package conveniently expects per-sample position and methylation information in a format already present in bedGraph files. However, this information is in a slightly different format than bedGraph, so the following awk script can be used. Note that BEAT expects files named as sample_name.positions.csv.
awk '{if(NR>1){printf("%s,%i,%i,%i\n",$1,$2+1,$5,$6)}else{printf("chr,pos,meth,unmeth\n")}}' sample.bedGraph > sample.positions.csvbison_herd has the ability to use a semi-arbitrary number of nodes. In practice,
if bison is given N nodes, it will effectively use 2*((N-1)/2)+1 or
4*((N-1)/4)+1 nodes, for directional and non-directional libraries,
respectively. As an example, if you allot 20 nodes for a directional library,
bison_herd will only use 19 of them (17 for non-directional reads). The excess
nodes will exit properly and, unless you specify --quiet, produce an error
message.
The options -mp, -queue-size, and -@ are bison_herd-specific and deserve further
description.
-mp sets the number of threads that the master node will use to process alignments produced by the worker nodes. Worker nodes are grouped into twos or fours, where each group has the a number of nodes equal to the number of possible bisulfite converted strands. As the number of allocated nodes increases, a point is eventually reached where a single thread on the master node is unable to keep up with the workers. In my experience, for directional libraries, one thread can handle approximately 130 bowtie2 threads (i.e., if using -p 11, -mp should be increased once ~12 worker nodes are allocated, since that would equate to 132 threads in use by bowtie2). One should keep in mind that there are already at least 3 other threads concurrently running on the master node (sending and storing fastq reads, receiving alignments, and writing alignments). Consequently, there is a practical limit to the number of nodes is determined by how many cores are available on each node.
-queue-size determines the maximum difference between reads sent for alignment
and reads processed. This option is unavailable if bison_herd was compiled with
-DNOTHROTTLE. By default, the thread that sends reads for alignment will pause
if it has sent more than ~1 million reads than have been processed. The purpose
of this is to prevent overwhelming of the MPI unexpected message buffer, since
the thread on the master node that sends reads can generally process reads
faster than all of the worker nodes combined can align them. Setting this value
too high may result in bison_herd crashing with otherwise cryptic messages
involving MPI_Send. In such cases, decreasing the value used by -queue-size
should resolve the problem. On the other hand, setting this value too low can
result in a deadlocks, due to buffering at various levels. The default value
hasn't resulted in deadlocking or crashes on our cluster, but yours may be
different! This difference is checked every 100000 reads, which can changed by
editting the THROTTLE_CHECK_INTERVAL value in bison.h prior to compilation.
-@ specifies the number of compression threads used for writing the output BAM file. In practice, a single compression thread can write ~80 million paired-end reads per hour (depending on CPU speed). I routinely use -@ 4 when using more than ~9 nodes as this allows writing to occur as quickly as reads are processed. To determine if the number of compression threads should be increased, not the time difference (especially early on) between when each master processor thread has processed 100000 reads and when those reads have been written to a file. Even when --reorder is used, if there is >1 second between these, then you may benefit from increasing the number of compression threads. For those curious, this option is identical to that used in samtools.
bison_herd generally uses blocking, but not synchronous sends. What this means
in practice is that many reads will be queued by the master node for sending to
the worker nodes. Likewise, many alignments can be queued by the worker nodes
for sending back to the master node. The queue that many MPI implementations use
for this is relatively small and immutable. While a full queue should cause
MPI_Send to block until there is sufficient space, occasionally a constellation
of events can occur that cause this queue to overflow and the master node to
then crash. This can be alleviated by limiting the possible number of reads that
could ever possibly be in the queue at any single time. As the queue is not
directly pollable, the difference between the number of reads sent and written
is used as a surrogate. The maximum number of reads in the wild is then either
2x or 4x this difference (since a read is queued per worker node). In reality,
the queue should be emptier than this as there are normally reads buffered on
the worker nodes (being fed to bowtie2, being aligned or being sent) and
elsewhere on the master node (being received, waiting to be processed, being
processed, waiting to be written, or being written).
Throttling is not always required, particularly as an increasing number of nodes are used. Throttling can be disabled altogether by compiling with -DNOTHROTTLE, which will remove all related components.
For debugging, a special debug mode is available for both bison and bison_herd
by compiling with -DDEBUG. Instead of running of needing multiple nodes, both
programs will then run as if they were just a single node. Compiling with this
option adds the -taskid option to both programs. The taskid is equivalent to the
node number in the bison (or bison_herd) hierarchy. Node 0 is the master node
and performs the final file writing. For bison, nodes 1-4 are equivalent to the
worker nodes that align reads to the original top, original bottom,
complementary to original top and complementary to original bottom strands,
respectively. For directional libraries, only the first 2 are used. These will
write alignments to a file for final processing when run as taskid 0. This is
useful when odd alignments are being output and the source of the error needs to
be tracked down. The mode for bison_herd is similar, except there are always 8
theoretical worker nodes (i.e., taskid 1-8 need to be run prior to taskid 0).
This allows testing multiple master processor threads with both directional and
non-directional reads.
In general, this mode should not be used unless you are running into extremely odd bugs.
Bison is generally similar to bismark, however the indexes are incompatible,
due to bismark renaming contigs. Also, the two will not produce identical
output, due to algorithmic differences. Running bison_methylation_extractor
on the output of bismark will also produce different results, again due to
algorithmic differences. In addition, bison always outputs BAM files directly.
Bison needn't be run on multiple computers. You can also use a single
computer for all compute nodes (e.g. mpiexec -n 5 bison ...). The same holds
true for bison_herd. Both bison and bison_herd seem to be faster than bismark,
even when limited to the same resources.
Allow lower case reads in fastq files (previously, this would result in corrupt BAM files.
HTSlib is now a submodule in the Github repository. This simplifies compilation. Further, that is the only supported compilation method now (samtools-0.1.19 is no longer supported).
Somehow, the methylation extractor was still defaulting to a minimum phred score of 10, when the documentation said it was defaulting to 5.
CRAM files can now be produced and processed. Both bison and bison_herd
will output in CRAM format if the -C option is given.
The header @PG line is now rewritten to contain the actual command executed and the bison/herd version.
Excess space allocated to hold the genome is now returned.
Output BAM/CRAM files can now be sorted on the fly. The method for this is
similar to that used by samtools, where temporary files are written and
then merged. This merge step is performed in parallel if multiple output
files are being written by bison_herd.
Fixed a bug in bison_CpG_coverage, where previously only the first chromosome was used.
Added bedGraph2MOABS to convert bedGraph files for use by MOABS. See usage above.
Added support for HTSlib.
Fixed a small bug wherein --reorder wasn't being invoked when multiple output BAM files were to be used.
Fixed a small bug that only manifested in DEBUG mode.
There is now a tutorial.
The default minimum MAPQ and Phred scores used by bison_mbias have been
updated to match bison_methylation_extractor.
The various bedGraph files didn't previously have a "track" line. The UCSC Genome Browser requires this, so bedGraph files produced will now contain it. It should be noted that this is the very minimal line required. Bison does not provide facilities for making these changes, users need to edit things manually or use external programs for this. It should also be noted that any changes to the "track" or other header lines should be made after all processing with Bison is complete.
Add conversion scripts for import into MethylSeekR, BiSeq, and BEAT.
Revamped how bison_markduplicates works. The 3' coordinates are now
ignored, soft-clipped bases on the 5' end are now incorporated in
determining the 5' coordinate and methylation calls are also used in
determining if reads/pairs are duplicates. This should be a much more
robust (though more resource intensive) method than that previously used.
Whereas the previous version kept unmarked the read/pair with the highest
MAPQ, this one will do that for the read/pair with the highest summed phred
score (a la picard).
Note: The indices produced by previous versions are not guaranteed to be
compatible unless you used a multi-fasta file. There was a serious
implementation problem with how bison_index worked when given multiple
files as input and how multiple files were read into memory in previous
versions. If you used a multi-fasta file, then everything will continue
to work correctly. However, if you used multiple fasta files in a list
then I strongly encourage you to delete the previous indices (just remove
the bisulfite_genome directory) and reindex. The technical reasons for this
issue are that when the bison tools previously read multiple fasta files
into memory, they would do so in whatever order they appeared in the
directory structure, which can change over time and isn't guaranteed to
match the order of files someone specified during indexing. While the
alignments wouldn't be affected by this, the methylation calls could have
been seriously compromised. In this version, bison_index will only accept a
directory, not a list of files, and it will always alphasort() the list of
files in that directory prior to processing. This should eliminate this
problem. My apologies to anyone affected by this.
Added --genome-size option to a number of the tools. Many of the bison programs need to read the genome into memory. By default, 3 gigabases worth of memory are allocated for that and the size increased as needed. For smaller genomes, this wasted space. For larger genomes, the constant reallocation of space could seriously slow things down. Consequently, this option was added to any tool that reads the genome into memory. It's convenient to overestimate this slightly, so if your genome is 3.8 gigabases, then just use 4000000000 as the genome size.
bison_merge_CpGs can now take multiple input files at once.
A number of small bug fixes, such as when "genome_dir" doesn't end in a /.
Fixed an off-by-one error in bison_mbias. Also, at some point 1-methylation percentage started getting calculated. That's been fixed.
Added bison_markduplicates, which, as the name implies, marks apparent PCR duplicates. The methylation extractor and m-bias calculator have also been updated to ignore marked duplicates.
Fixed a bug in the CpG coverage program, which wasn't properly handling single-C bedGraph files before (if they were merged, then they were being handled correctly).
Fix how hard and soft-clipped bases are dealt with (previously, soft- clipped bases resulted in an error and hard-clipped bases in incorrect position assignments!).
Multiple bug fixes related to local alignment, which previously didn't work correctly. These issues seem to generally now be resolved. May thanks to user mvijayen on seqanswers for providing a perfect usage example for testing (see thread http://seqanswers.com/forums/showthread.php?t=39914).
The maximum length of a single contig is now (2^64)-1 (instead of the previous 2^64). I don't think bowtie2 would even support something that long, but if it did then bison wouldn't (internally, a position of 2^64 means a base is inserted, soft, or hard-clipped).
A previously missing "*" caused Bison to use the entirety of the description line in the fasta file as the chromosome name. This caused errors since bowtie2 only uses every before the first space (the proper method). Bison now does the same.
A note about creating methylation-bias metrics with locally aligned reads is in order. If a read is soft-clipped, that portion is still included in the M-bias metrics. Likewise, if you pass -OT X,X,X,X or similar parameters to the methylation extractor, the soft-clipped area is also included in there.
Another note regarding local alignments is that the XX auxiliary tag (effectively the more verbose version of the MD tag) contains soft-clipped sequences. I could probably have these removed if someone would like.
Properly fixed some wording on the textual output (i.e., removed the word "unique").
Lowered the default MAPQ and Phred thresholds used by the methylation extractor to 10 each. That the MAPQ threshold was originally 20 was an error on my part.
Added support for file globbing in bison_herd. You may now input multiple files using a combination of wild-cards (, ?, etc.) and commas. Remember to put these in quotes (e.g., "foo/1.fq.gz","bar/*1.fq.gz") so the shell doesn't perform the expansion!). As before, specifying multiple inputs with the same file name (e.g., sample1/reads.fq,sample2/reads.fq) will cause the output from the first reads.fq alignment to be over-written by the second.
Fixed the text output, since "unique alignments" isn't really correct, given that alignments with scores of 0 or 1 can be output but aren't unique.
Added information in the Makefile and above about compiling with openmpi.
Fixed a bug in bison_herd wherein the -upto option wasn't being handled properly. -upto now accepts an unsigned long in bison_herd.
Fixed a bug in bison_herd when paired-end reads were used. This was due to how bowtie2 reads from FIFOs. Changing how things were written to the FIFOs on the worker nodes resolved the problem.
The bison_mbias program has been heavily revamped. It still outputs the number of methylated or unmethylated CpG calls per position, but now keeps the metrics for each strand (and read, when paired-end reads are used) separate. If R and the ggplot2 library are installed, the program can also run the bison_mbias2pdf program (see below).
Created an bison_mbias2pdf Rscript that will read in the output of bison_mbias and plot the results, indicating the region of each read that should be included in methylation extraction. This script also print these suggestions in the format used by bison_methylation_extractor, for convenience.
The methylation extractor can now be told to only include certain regions of each read in the output methylation metrics. This is needed when there is apparent bias in the methylation at one or both ends of a read.
Previously, the recalculated MAPQ was incorrect when only 1 read in a pair had a valid secondary alignment. This has been fixed.
Fixed another MAPQ recalculation bug, affecting reads with MAPQ 2 that have MAPQ=6.
Fixed a bug in writing unmapped reads.
Fixed a bug in bison_herd that allowed early termination without warning.
Added a note to the methylation summary statistics output at the end of a run that the numbers will include double counting of any site covered by both mates in a pair. These metrics are only meant for general information and not further analysis, so I don't consider that a bug (it's actually a design decision for the sake of performance).
--ignore-quals is no longer passed to bowtie2 by default. Specifying this will marginally decrease both correct and incorrect alignments. It will also generally decrease the alignment rate.
Fixed --unmapped, which are now written to the directory specified by -o
--maxins was already 500 by default, so it is no longer set by default.
Added bison_herd, see above for usage
The methylation extractor now has a -phred option, to exclude methylation calls from low confidence base-calls. The default threshold is 20.
Added a script to convert bedGraph files to a format suitable for BSseq.
Fixed a bug in bison_merge_CpGs
Both bison and bison_herd now check to ensure that the MPI implementation actually supports the level of thread support requested (previously, this was just assumed).
Fixed a number of minor bugs.
Added support for uncompressed fastq files, as well as bzipped files (previously, only gzipped fastq files worked properly).
--score-min is now parsed by bison prior to being sent to bowtie2, read MAPQ scores are recalculated accordingly by the same algorithm used by bowtie2 (N.B., this only bears a vague correspondence to -10*log10(probability the mapping position is wrong)!).
Added a bison_mbias function, to process the aligned BAM file and create a text file containing the percentage of methylated C's as a function of read position. For the utility of this, see: Hansen KD, Langmead B and Irizarry RA, BSmooth: from whole genome bisulfite sequencing to differentially methylated regions. Genome Biol 2012; 13(10):R83.
The methylation extractor now accepts the -q options, which sets the MAPQ threshold for a read to be included in the methylation results. The default is a minimum MAPQ of 20, which seems to be a reasonable threshold from a few simulations.
In DEBUG mode, the output BAM files used to have fixed names. This was a problem in cases where debugging was being performed on multiple input files. Now, the OT/OB/CTOT/CTOB.bam filename is prepended with an appropriate prefix (extracted from the input file name). In addition, the output directory is now respected in DEBUG mode.
Included an "auxiliary" directory, that includes functions for making an RRBS genome and other possibly useful functions.
Initial release