bam2wig.pl
A program to convert Bam alignments into a wig representation file.
bam2wig.pl [--options...] <file.bam>
bam2wig.pl --extend --rpm --mean --out file --bw file1.bam file2.bam
Required options: -i --in <filename.bam> repeat if multiple bams, or comma-delimited list Reporting options (pick one): -s --start record at 5' position -d --mid record at midpoint of alignment or pair -a --span record across entire alignment or pair -e --extend extend alignment (record predicted fragment) --cspan record a span centered on midpoint --smartcov record paired coverage without overlaps, splices --ends record paired endpoints --coverage raw alignment coverage Alignment reporting options: -l --splice split alignment at N splices -t --strand record separate strands as two wig files --flip flip the strands for convenience Paired-end alignments: -p --pe process paired-end alignments, both are checked -P --fastpe process paired-end alignments, only F are checked --minsize <integer> minimum allowed insertion size (30) --maxsize <integer> maximum allowed insertion size (600) --first only process paired first read (0x40) as single-end --second only process paired second read (0x80) as single-end Alignment filtering options: -K --chrskip <regex> regular expression to skip chromosomes -B --blacklist <file> interval file of regions to skip (bed, gff, txt) -q --qual <integer> minimum mapping quality (0) -S --nosecondary ignore secondary (0x100) alignments (false) -D --noduplicate ignore duplicate (0x400) alignments (false) -U --nosupplementary ignore supplementary (0x800) alignments (false) --intron <integer> maximum allowed gap (intron) size in bp (none) Shift options: -I --shift shift reads in the 3' direction -x --extval <integer> explicit extension size in bp (default is to calculate) -H --shiftval <integer> explicit shift value in bp (default is to calculate) --chrom <integer> number of chromosomes to sample (4) --minr <float> minimum pearson correlation to calculate shift (0.5) --zmin <float> minimum z-score from average to test peak for shift (3) --zmax <float> maximum z-score from average to test peak for shift (10) -M --model write peak shift model file for graphing Score options: -r --rpm scale depth to Reads Per Million mapped -m --mean average multiple bams (default is addition) --scale <float> explicit scaling factor, repeat for each bam file --fraction assign fractional counts to multi-mapped alignments --splfrac assign fractional count to each spliced segment --format <integer> number of decimal positions (4) --chrnorm <float> use chromosome-specific normalization factor --chrapply <regex> regular expression to apply chromosome-specific factor Wig format: --bin <integer> bin size for span or extend mode (10) --bdg bedGraph, default for span and extend at bin 1 --fix fixedStep, default for bin > 1 --var varStep, default for start, mid --nozero do not write zero score intervals in bedGraph Output options: -o --out <filename> output file name, default is bam file basename -b --bw convert to bigWig format (supports bdg, fix, var) --bwapp /path/to/wigToBigWig path to external converter (default searches \$PATH) -z --gz gzip compress text output General options: -c --cpu <integer> number of parallel processes (4) --temp <directory> directory to write temporary files (output path) -V --verbose report additional information -v --version print version information -h --help show full documentation
The command line flags and descriptions:
Specify the input Bam alignment file. More than one file may be specified, either with repeated options, a comma-delimited list, or simply appended to the command. Bam files will be automatically indexed if necessary.
Specify that the 5' position should be recorded in the wig file.
Specify that the midpoint of the alignment (single-end) or fragment (paired-end) will be recorded in the wig file.
Specify that the entire span of the alignment (single-end) or fragment (paired-end) will be recorded in the wig file.
Specify that the alignment should be extended in the 3' direction and that the entire length of the extension be recorded in the wig file. The extension may be defined by the user or empirically determined.
Specify that a defined span centered at the alignment (single-end) or fragment (paired-end) midpoint will be recorded in the wig file. The span is defined by the extension value.
Smart alignment coverage of paired-end alignments without double-counting overlaps or recording gaps (intron splices).
Record both endpoints of paired-end fragments, i.e. the outermost or 5' ends of properly paired fragments. This may be useful with ATAC-Seq, Cut&Run-Seq, or other cleavage experiments where you want to record the locations of cutting yet retain the ability to filter paired-end fragment sizes.
Specify that the raw alignment coverage be calculated and reported in the wig file. This utilizes a special low-level operation and precludes any alignment filtering or post-normalization methods. Counting overlapping bases in paired-end alignments are dependent on the bam adapter (older versions would double-count).
Legacy option for supporting previous versions of bam2wig.
Indicate that the bam file contains alignments with splices, such as from RNASeq experiments. Alignments will be split on cigar N operations and each sub fragment will be recorded. This only works with single-end alignments, and is disabled for paired-end reads (just treat as single-end). Only start and span recording options are supported.
Indicate that separate wig files should be written for each strand. The output file basename is appended with either '_f' or '_r' for both files. Strand for paired-end alignments are determined by the strand of the first read.
Flip the strand of the output files when generating stranded wig files. Do this when RNA-Seq alignments map to the opposite strand of the coding sequence, depending on the library preparation method.
The Bam file consists of paired-end alignments, and only properly mapped pairs of alignments will be counted. Properly mapped pairs include FR reads on the same chromosome, and not FF, RR, RF, or pairs aligning to separate chromosomes. Both alignments are required to be present before the pair is counted. The default is to treat all alignments as single-end.
The Bam file consists of paired-end alignments, but to increase processing time and be more tolerant of weird pairings, only the forward alignment is required and considered; all reverse alignments are ignored. The default is to treat all alignments as single-end.
Specify the minimum paired-end fragment size in bp to accept for recording. Default is 30 bp.
Specify the maximum paired-end fragment size in bp to accept for recording. Default is 600 bp.
Take only the first read of a pair, indicated by flag 0x40, and record as a single-end alignment. No test of insert size or proper pair status is made.
Take only the second read of a pair, indicated by flag 0x80, and record as a single-end alignment. No test of insert size or proper pair status is made.
Set a minimum mapping quality score of alignments to count. The mapping quality is a range from 0-255, with higher numbers indicating lower probability of a mapping error. Multi-mapping alignments often have a map quality of 0. The default is 0 (accept everything).
Boolean flag to skip secondary alignments, indicated by the alignment bit flag 0x100. Secondary alignments typically represent alternative mapping locations, or multi-mapping events. By default, secondary alignments are included.
Boolean flag to skip duplicate alignments, indicated by the alignment bit flag 0x400. Duplicates alignments may represent a PCR or optical duplication. By default, duplicate alignments are included.
Boolean flag to skip supplementary alignments, indicated by the alignment bit flag 0x800. Supplementary alignments are typically associated with chimeric fragments. By default, supplementary alignments are included.
Provide a regular expression to skip certain chromosomes. Perl-based regular expressions are employed. Expressions should be quoted or properly escaped on the command line. Examples might be
'chrM' 'scaffold.+' 'chr.+alt|chrUn.+|chr.+_random'
Provide a file of genomic intervals from which to exclude alignments. Examples might include repeats, ribosomal RNA, or heterochromatic regions. The file should be any text file interpretable by Bio::ToolBox::Data with chromosome, start, and stop coordinates, including BED and GFF formats. Note that this only excludes overlapping alignments, and does not include extended alignments.
Provide a positive integer as the maximum intron size allowed in an alignment when splitting on splices. If an N operation in the CIGAR string exceeds this limit, the alignment is skipped. Default is 0 (no filtering).
Specify that the positions of the alignment should be shifted towards the 3' end. Useful for ChIP-Seq applications, where only the ends of the fragments are counted and often seen as separated discrete peaks on opposite strands flanking the true target site. This option is disabled with paired-end and spliced reads (where it is not needed).
Provide the value in bp that the recorded position should be shifted. The value should be 1/2 the average length of the library insert size. The default is to automatically and empirically determine the appropriate shift value using cross-strand correlation (recommended).
Manually set the length for reads to be extended. By default, the shift value is determined empirically and extension is set to 2X the shift value. This is also used for the cspan mode.
Indicate the number of sequences or chromosomes to sample when empirically determining the shift value. The reference sequences listed in the Bam file header are taken in order of decreasing length, and one or more are taken as a representative sample of the genome. The default value is 4.
Provide the minimum Pearson correlation value to accept a shift value when empirically determining the shift value. Enter a decimal value between 0 and 1. Higher values are more stringent. The default is 0.5.
Specify the minimum z-score (or number of standard deviations) from the chromosomal mean depth to test for a peak shift. Increase this number to test for strong robust peaks, which give a better estimations of the shift value. Default is 3.
Specify the maximum z-score (or number of standard deviations) from the chromosomal mean depth to test for a peak shift. This excludes erroneous peaks due to repetitive sequence alignments with high coverage. Increase this number to include more robust peaks that can give a better estimation of the shift value. Default is 10.
Indicate that the shift model profile data should be written to file for examination. The average profile, including for each sampled chromosome, are reported for the forward and reverse strands, as well as the shifted profile. A standard text file is generated using the output base name. The default is to not write the model shift data.
Convert the data to Reads (or Fragments) Per Million mapped. This is useful for comparing read coverage between different datasets. The default is no RPM conversion.
Optionally provide your own scaling factor. This will be multiplied with every position when generating the wig file. This may be combined with the rpm factor as well. When combining multiple bam files, either a single scale factor may be supplied for all files, or individual scale factors may be supplied for each bam file. If supplying multiple, use the option multiple times or give a comma-delimited list. The values should be in the same order as the bam files.
When processing multiple bam files, this option will take the mean or average across all bam files. Without this option, the bam files are simply added. When combined with the rpm option, each bam file will be scaled separately before taking the average.
Indicate that multi-mapping alignments should be given fractional counts instead of full counts. The number of alignments is determined using the NH alignment tag. If a read has 10 alignments, then each alignment is given a count of 0.1.
Indicate that spliced segments should be given a fractional count. This allows a count to be assigned to each spliced segment while avoiding double-counting. Best used with RNASeq spliced point data (--start or --mid); not recommended for --span.
Indicate the number of decimal postions reported in the wig file. This is only applicable when rpm, scale, or fraction options are provided. The default value is 4 decimal positions.
Apply a normalization factor to the counts on specific chromosomes only. Usually this is to normalize, for example, variable copy-number chromosomes, such as transfected vector sequences, or haploid sex chromosomes.
Specify the Perl-based regular expression to match the chromosome(s) to apply the specific normalization factor. For example, 'chrX$' to specify the X chromosome only.
Specify the bin size in bp for the output wig file. In general, specifying a larger bin size will decrease the run time and memory requirements in exchange for loss of resolution. The default for span, center span, or extend mode is 10 bp; all other modes is 1 bp.
Specify that the output wig format is a bedGraph-style wig format. This is the default format for extend, span, and cspan modes of operation.
Specify that the output wig format is in fixedStep wig format. This is the default format for coverage mode of operation.
Specify that the output wig format is in variableStep wig format. This is the default format for start and midpoint modes of operation.
When writing bedGraph format, skip (do not write) intervals with a value of zero. Does not apply to fixedStep or variableStep formats.
Specify the output base filename. An appropriate extension will be added automatically. By default it uses the base name of the input file.
Specify whether or not the wig file should be further converted into an indexed, compressed, binary BigWig file. The default is false.
Optionally specify the full path to the UCSC wigToBigWig conversion utility. The application path may be set in the .biotoolbox.cfg file or found in the default environment $PATH, which makes this option mostly unnecessary.
$PATH
Specify whether (or not) the output text file should be compressed with gzip. Disable with --nogz. Does not apply to bigWig format.
--nogz
Specify the number of parallel instances to run simultaneously. This requires the installation of Parallel::ForkManager. With support enabled, the default is 4. Disable multi-threaded execution by setting to 1.
Optionally specify an alternate temporary directory path where the temporary files will be written. The default is the specified output file path, or the current directory. Temporary files will always be written in a subdirectory of the path specified with the template "bam2wigTEMP_XXXX".
Print extra informational statements during processing. The default is false.
Print the version number.
Display this POD documentation.
This program will enumerate aligned sequence tags and generate a wig, or optionally BigWig, file. Alignments may be counted and recorded in several different ways. Strict enumeration may be performed and recorded at either the alignment's start or midpoint position. Alternatively, either the alignment or fragment may be recorded across its span. Finally, a basic unstranded, unshifted, and non-transformed alignment coverage may be generated.
Both paired-end and single-end alignments may be counted. Alignments with splices (e.g. RNA-Seq) may be counted singly or separately. Alignment counts may be separated by strand, facilitating analysis of RNA-Seq experiments.
For ChIP-Seq experiments, the alignment position may be shifted in the 3 prime direction. This effectively merges the separate peaks (representing the ends of the enriched fragments) on each strand into a single peak centered over the target locus. Alternatively, the entire predicted fragment may be recorded across its span. This extended method of recording infers the mean size of the library fragments, thereby emulating the coverage of paired-end sequencing using single-end sequence data. The shift value is empirically determined from the sequencing data or provided by the user. If requested, the shift model profile may be written to file.
The output wig file may be either a variableStep, fixedStep, or bedGraph format. The wig file may be further converted into a compressed, indexed, binary bigWig format, dependent on the availability of the appropriate conversion utilities.
The type of wig file to generate for your Bam sequencing file can vary depending on your particular experimental application. Here are a few common sequencing applications and my recommended settings for generating the wig or bigWig file.
To generate a straight-forward coverage map, similar to what most genome browsers display when using a Bam file as source. NOTE that this mode is pure raw coverage, and does not include any filtering methods. The other modes allow alignment filtering.
bam2wig.pl --coverage --in <bamfile>
When you have paired-end alignments and need explicit alignment coverage without double-counting overlaps (as would occur if you counted as single-end span) or uncovered insertion (as would occur if you counted as paired-end span) and not counting gaps (e.g. intron splices, as would occur with span mode), use the smart paired-end coverage mode. This properly assembles coverage from paired-end alignments taking into account overlaps and gaps.
bam2wig --smartcov --in <bamfile>
When sequencing Chromatin Immuno-Precipitation products, one generally performs a 3 prime shift adjustment to center the fragment's end reads over the predicted center and putative target. To adjust the positions of tag count peaks, let the program empirically determine the shift value from the sequence data (recommended). Otherwise, if you know the mean size of your ChIP eluate fragments, you can use the --shiftval option.
To evaluate the empirically determined shift value, be sure to include the --model option to examine the profiles of stranded and shifted read counts and the distribution of cross-strand correlations.
Depending on your downstream applications and/or preferences, you can record strict enumeration (start positions) or coverage (extend position).
Finally, to compare ChIP-Seq alignments from multiple experiments, convert your reads to Reads Per Million Mapped, which will help to normalize read counts.
bam2wig.pl --start --shift --model --rpm --in <bamfile> bam2wig.pl --extend --model --rpm --in <bamfile>
If both ends of the ChIP eluate fragments are sequenced, then we do not need to calculate a shift value. Instead, we will simply count at the midpoint of each properly-mapped sequence pair, or record the defined fragment span.
bam2wig.pl --mid --pe --rpm --in <bamfile> bam2wig.pl --span --pe --rpm --in <bamfile>
With RNA-Sequencing, we may be interested in either coverage (generating a transcriptome map) or simple tag counts (differential gene expression), so we can count in one of two ways.
To compare RNA-Seq data from different experiments, convert the read counts to Reads Per Million Mapped, which will help to normalize read counts.
bam2wig --span --splice --rpm --in <bamfile> bam2wig --mid --rpm --in <bamfile>
If the library was generated in such a way as to preserve strand, then we can separate the counts based on the strand of the alignment. Note that the reported strand may be accurate or flipped, depending upon whether first-strand or second-strand synthesized cDNA was sequenced, and whether your aligner took this into account. Check the Bam alignments in a genome browser to confirm the orientation relative to coding sequences. If alignments are opposite to the direction of transcription, you can include the --flip option to switch the output.
bam2wig --span ---splice --strand --rpm --in <bamfile> bam2wig --pos mid --strand --rpm --in <bamfile>
Use the smart paired-end coverage mode to properly record paired-end alignments with splice junctions.
bam2wig --smartcov --strand --rpm --in <bamfile>
To help users visualize how this program records alignments in a wig file, drawn below are 10 alignments, five forward and five reverse. They may be interpreted as either single-end or paired-end. Drawn below are the numbers that would be recorded in a wig file for various parameter settings. Note that alignments are not drawn to scale and are drawn for visualization purposes only. Values of X represent 10.
....>>>>>>.....................................<<<<<<............. .....>>>>>>..................................<<<<<<............... ........>>>>>>.......................................<<<<<<....... ........>>>>>>.........................................<<<<<<..... ..........>>>>>>............................................<<<<<<
....11..2.1.......................................1.1.....1.1....1
......11..2.1..................................1.1.....1.1....1...
F...11..2.1....................................................... R.................................................1.1.....1.1....1
....122244433311.............................112222111122221211111
......121.2211.................................1111....1111...11..
F...122244433311.................................................. R............................................112222111122221211111
........................1.1...11121.1..1..........................
...................11222211112344365544411........................
12223445789999XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX999887665321111
............................1...111..1............................
..........................111123333432111.........................
....12224455555555555555555555555555555555555555555443333332211111
Timothy J. Parnell, PhD Huntsman Cancer Institute University of Utah Salt Lake City, UT, 84112
This package is free software; you can redistribute it and/or modify it under the terms of the Artistic License 2.0.
To install Bio::ToolBox, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Bio::ToolBox
CPAN shell
perl -MCPAN -e shell install Bio::ToolBox
For more information on module installation, please visit the detailed CPAN module installation guide.