Bio-ToolBox
view release on metacpan or search on metacpan
scripts/bam2wig.pl view on Meta::CPAN
'minsize=i' => \$min_isize, # minimum paired insert size to accept
'first!' => \$first_read, # only take first read
'second!' => \$second_read, # only take second read
'fraction!' => \$multi_hit_scale, # scale by number of hits
'splfrac!' => \$splice_scale, # divide counts by number of spliced segments
'r|rpm!' => \$rpm, # calculate reads per million
'm|separate|mean!' => \$do_mean, # rpm scale separately
'scale=s' => \@scale_values, # user specified scale value
'chrnorm=f' => \$chrnorm, # chromosome-specific normalization
'chrapply=s' => \$chrapply, # chromosome-specific normalization regex
'K|chrskip=s' => \$chr_exclude, # regex for skipping chromosomes
'exclude|blacklist=s' => \$exclude_file, # file for skipping regions
'bin=i' => \$bin_size, # size for binning the data
'format=i' => \$dec_precison, # format to number of decimal positions
'b|bw!' => \$bigwig, # generate bigwig file
'bwapp=s' => \$bwapp, # utility to generate a bigwig file
'bdg!' => \$do_bedgraph, # write a bedgraph output
'fix!' => \$do_fixstep, # write a fixedStep output
'var!' => \$do_varstep, # write a varStep output
'nozero!' => \$no_zero, # do not write zero coverage
'z|gz!' => \$gz, # compress text output
'c|cpu=i' => \$cpu, # number of cpu cores to use
'intron=i' => \$max_intron, # maximum intron size to allow
'window=i' => \$window, # window size to control memory usage
'V|verbose!' => \$verbose, # print sample correlations
'temp=s' => \$tempdir, # directory to write temp files
'adapter=s' => \$BAM_ADAPTER, # explicitly set the adapter version
'h|help' => \$help, # request help
'v|version' => \$print_version, # print the version
) or die " unrecognized option(s)!! please refer to the help documentation\n\n";
# Print help
if ($help) {
# print entire POD
pod2usage(
{
'-verbose' => 2,
'-exitval' => 1,
}
);
}
# Print version
if ($print_version) {
print " Biotoolbox script bam2wig.pl, version $VERSION\n";
my $v = Bio::ToolBox->VERSION;
print " Biotoolbox package version $v\n";
exit;
}
### Check for requirements and set defaults
# more global variables
my (
$unwanted_flags, $main_callback, $callback, $wig_writer,
$outbase, $chromo_file, $binpack, $buflength,
$coverage_dump, $coverage_sub, $post_bw_convert
);
check_defaults();
print " Writing temp files to $tempdir\n" if $verbose;
my $items = $paired ? 'fragments' : 'alignments';
# record start time
my $start_time = time;
### Open files
my @sams;
printf " Processing files %s...\n", join( ", ", @bamfiles );
foreach (@bamfiles) {
# this will open each bam file using the high level API
# with the appropriate installed adapter
my $sam = open_db_connection($_) or die " unable to open bam file '$_'!\n";
$sam->split_splices(1) if $splice;
push @sams, $sam;
}
# generate the chromosome name list
# we generate this from the first bam file only, on the assumption that they
# all have the same sequences
# record the chromosome name, rather than internal tid, just in the off chance
# that they are not in the same order (!!!???)
my @seq_list;
my %seq_name2length;
for my $tid ( 0 .. $sams[0]->n_targets - 1 ) {
my $chr = $sams[0]->target_name($tid);
if ( $chr_exclude and $chr =~ /$chr_exclude/xi ) {
print " skipping sequence $chr\n" if $verbose;
next;
}
push @seq_list, $chr;
$seq_name2length{$chr} = $sams[0]->target_len($tid);
}
# set the wrapper reference
# this depends on which adapter was opened
my $wrapper_ref;
if ( $splice or $use_smartpe ) {
$wrapper_ref =
ref( $sams[0] ) eq 'Bio::DB::Sam' ? 'Bio::DB::Bam::AlignWrapper'
: ref( $sams[0] ) eq 'Bio::DB::HTS' ? 'Bio::DB::HTS::AlignWrapper'
: 'none';
eval { require $wrapper_ref; 1 };
}
printf " Using the %s Bam adapter and align wrapper $wrapper_ref\n", ref( $sams[0] )
if $verbose;
### Process user provided exclusion lists
my $exclude_hash = process_exclusion_list();
### Calculate shift value
if ( $shift or $use_extend ) {
unless ( ( $shift and $shift_value ) or ( $use_extend and $extend_value ) ) {
print " Calculating 3' shift value...\n";
$shift_value = determine_shift_value();
}
# quietly exit here after determining shift value if no wig is to be generated
exit
unless ( $use_start
or $use_extend
scripts/bam2wig.pl view on Meta::CPAN
unless ( $shift_value or $extend_value ) {
# not provided by user, empirical calculation required
my $stat = 0;
eval {
# required for calculating shift
require Statistics::Descriptive;
$stat = 1;
};
unless ($stat) {
print STDERR
" FATAL: Provide a shift value or install the Perl module Statistics::Descriptive\n"
. " to empirically determine the shift value.\n";
exit 1;
}
}
if ( defined $correlation_min ) {
if ( $correlation_min <= 0 or $correlation_min >= 1 ) {
print STDERR
" FATAL: cannot use minimum correlation value of $correlation_min!\n use --help for more information\n";
exit 1;
}
}
else {
$correlation_min = 0.5;
}
$chr_number ||= 4;
$zmin ||= 3;
$zmax ||= 10;
}
# center span should have extend value
if ( $use_cspan and not $extend_value ) {
print STDERR
" FATAL: please use the --extval option to define an extend value when using center span\n";
exit 1;
}
# check mapping quality
if ( defined $min_mapq ) {
if ( $min_mapq > 255 or $min_mapq < 0 ) {
print STDERR " FATAL: quality score must be 0-255!\n";
exit 1;
}
}
else {
$min_mapq = 0;
}
# check paired-end insert size
unless ( defined $max_isize ) {
if ($use_smartpe) {
# Smart paired-end coverage doesn't care at all about the size of
# the insertion, but it is nevertheless tested in the pe_callback to
# accomodate other functions. So set this to a reasonably really high number.
$max_isize = 100000; # 100 kb should be sufficiently high
}
else {
# vast majority of paired-end fragments are less than 600 bp
# really, really big fragments are most likely mapping errors
$max_isize = 600;
}
}
unless ( defined $min_isize ) {
$min_isize = 30;
}
# chromosome-specific normalization
if ( $chrnorm or $chrapply ) {
# must have both of these parameters
unless ($chrnorm) {
print STDERR
" FATAL: missing --chrnorm value a for specific-chromosome normalization!\n";
exit 1;
}
unless ($chrapply) {
print STDERR
"FATAL: missing --chrapply regex for specific-chromosome normalization!\n";
exit 1;
}
}
# check flag parameters
if ($verbose) {
printf " %s secondary 0x100 reads\n", $nosecondary ? 'Skipping' : 'Including';
printf " %s duplicate 0x400 reads\n", $noduplicate ? 'Skipping' : 'Including';
printf " %s supplementary 0x800 reads\n",
$nosupplementary ? 'Skipping' : 'Including';
}
# set bin size
if ($bin_size) {
if ( $bin_size < 0 ) {
print STDERR " FATAL: bin size cannot be negative!\n";
exit 1;
}
}
else {
# set default to 10 bp for any span or coverage, or 1 bp for point data
$bin_size = ( $use_start or $use_mid or $use_ends ) ? 1 : 10;
}
# determine binary file packing and length
if ( $multi_hit_scale
or $splice_scale
or $chrnorm
or ( $use_coverage and $bin_size > 1 ) )
{
# pack as floating point values, this is 32 bit
$binpack = 'f';
}
else {
# dealing only with integers here
# yes we do occasionally have depth greater than 65,536
# originally short (16 bits), now long (32 bits)
$binpack = 'L';
}
$buflength = length( pack( $binpack, 1 ) );
scripts/bam2wig.pl view on Meta::CPAN
print " Recording single-end shifted, stranded alignment midpoints\n";
}
# span
elsif ( not $paired and not $do_strand and not $shift and $use_span ) {
$callback = \&se_span;
print " Recording single-end alignment span\n";
}
elsif ( not $paired and not $do_strand and $shift and $use_span ) {
$callback = \&se_shift_span;
print " Recording single-end shifted alignment span\n";
}
elsif ( not $paired and $do_strand and not $shift and $use_span ) {
$callback = \&se_strand_span;
print " Recording single-end stranded alignment span\n";
}
elsif ( not $paired and $do_strand and $shift and $use_span ) {
$callback = \&se_shift_strand_span;
print " Recording single-end shifted, stranded alignment span\n";
}
# center span
elsif ( not $paired and not $do_strand and not $shift and $use_cspan ) {
$callback = \&se_center_span;
print " Recording single-end alignment center-span\n";
}
elsif ( not $paired and not $do_strand and $shift and $use_cspan ) {
$callback = \&se_shift_center_span;
print " Recording single-end shifted alignment center-span\n";
}
elsif ( not $paired and $do_strand and not $shift and $use_cspan ) {
$callback = \&se_strand_center_span;
print " Recording single-end stranded alignment center-span\n";
}
elsif ( not $paired and $do_strand and $shift and $use_cspan ) {
$callback = \&se_shift_strand_center_span;
print " Recording single-end shifted, stranded alignment center-span\n";
}
# extend
elsif ( not $paired and not $do_strand and not $shift and $use_extend ) {
$callback = \&se_extend;
print " Recording single-end extended alignment span\n";
}
elsif ( not $paired and not $do_strand and $shift and $use_extend ) {
$callback = \&se_shift_extend;
print " Recording single-end shifted, extended alignment span\n";
}
elsif ( not $paired and $do_strand and not $shift and $use_extend ) {
$callback = \&se_strand_extend;
print " Recording single-end stranded, extended alignment span\n";
}
elsif ( not $paired and $do_strand and $shift and $use_extend ) {
$callback = \&se_shift_strand_extend;
print " Recording single-end shifted, stranded, extended alignment span\n";
}
# paired-end start
elsif ( $paired and not $do_strand and $use_start ) {
$callback = \&pe_start;
print " Recording paired-end fragment start\n";
}
elsif ( $paired and $do_strand and $use_start ) {
$callback = \&pe_strand_start;
print " Recording paired-end stranded fragment start\n";
}
# paired-end midpoint
elsif ( $paired and not $do_strand and $use_mid ) {
$callback = \&pe_mid;
print " Recording paired-end fragment midpoint\n";
}
elsif ( $paired and $do_strand and $use_mid ) {
$callback = \&pe_strand_mid;
print " Recording paired-end stranded fragment midpoint\n";
}
# paired-end span
elsif ( $paired and not $do_strand and $use_span ) {
$callback = \&pe_span;
print " Recording paired-end fragment span\n";
}
elsif ( $paired and $do_strand and $use_span ) {
$callback = \&pe_strand_span;
print " Recording paired-end stranded, fragment span\n";
}
# paired-end center span
elsif ( $paired and not $do_strand and $use_cspan ) {
$callback = \&pe_center_span;
print " Recording paired-end fragment center-span\n";
}
elsif ( $paired and $do_strand and $use_cspan ) {
$callback = \&pe_strand_center_span;
print " Recording paired-end stranded, fragment center-span\n";
}
# paired-end smart coverage
elsif ( $paired and not $do_strand and $use_smartpe ) {
$callback = \&smart_pe;
print " Recording smart paired-end coverage\n";
}
elsif ( $paired and $do_strand and $use_smartpe ) {
$callback = \&smart_stranded_pe;
print " Recording stranded, smart paired-end coverage\n";
}
elsif ( $paired and not $do_strand and $use_ends ) {
$callback = \&pe_ends;
print " Recording paired-end fragment endpoints\n";
}
elsif ( $paired and $do_strand and $use_ends ) {
$callback = \&pe_strand_ends;
print " Recording stranded, paired-end fragment endpoints\n";
}
else {
die "programmer error!\n" unless $shift; # special exception
}
# summary of wig file being written
if ($do_bedgraph) {
print " Writing bedGraph format in $bin_size bp increments\n";
}
elsif ($do_varstep) {
printf " Writing variableStep format in $bin_size bp bins\n";
}
elsif ($do_fixstep) {
printf " Writing fixedStep format in $bin_size bp bins\n";
}
}
sub process_exclusion_list {
if ( $exclude_file and -e $exclude_file ) {
my $i = 0;
eval { require Set::IntervalTree; $i = 1; };
unless ($i) {
print " WARNING! Please install Set::IntervalTree to use exclusion lists\n";
undef $exclude_file;
return;
}
my %list_hash = map { $_ => [] } @seq_list;
my $Data = Bio::ToolBox->load_file($exclude_file)
or die "unable to read exclusion list file '$exclude_file'\n";
unless ( $Data->feature_type eq 'coordinate' ) {
printf
" WARNING! Exclusion list file '%s' does not have coordinates! Ignoring.\n",
$exclude_file;
return \%list_hash;
}
$Data->iterate(
sub {
my $row = shift;
push @{ $list_hash{ $row->seq_id } }, [ $row->start - 1, $row->end ]
if exists $list_hash{ $row->seq_id };
}
);
printf " Loaded %s exclusion list regions\n",
format_with_commas( $Data->number_rows );
return \%list_hash;
}
return;
}
### Determine the shift value
sub determine_shift_value {
# identify top regions to score
# we will walk through the largest chromosome(s) looking for the top
# 500 bp regions containing the highest unstranded coverage to use
print " sampling the top coverage regions "
. "on the largest $chr_number chromosomes...\n";
# first sort the chromosomes by size
# this is assuming all the chromosomes have different sizes ;-)
scripts/bam2wig.pl view on Meta::CPAN
$data->{r}->[$_] += $score;
}
}
else {
die sprintf(
" Paired-end flags are set incorrectly; neither 0x040 or 0x080 are set for paired-end read %s.",
$a->query->name );
}
}
sub smart_pe {
my ( $f, $data, $score, $r ) = @_;
my $set = Set::IntSpan::Fast->new;
# process both reads, adding to the integer set
foreach my $a ( $f, $r ) {
if ( $a->cigar_str =~ /N/ ) {
my $aw = $wrapper_ref->new( $a, $data->{sam} );
# check intron size from the cigar string
if ($max_intron) {
my $size = 1;
my $cigars = $aw->cigar_array;
foreach my $c ( @{$cigars} ) {
# each element is [operation, size]
$size = $c->[1] if ( $c->[0] eq 'N' and $c->[1] > $size );
}
return if $size > $max_intron; # exceed maximum intron size
}
# record
foreach my $segment ( $aw->get_SeqFeatures ) {
$set->add_range(
int( ( $segment->start - 1 ) / $bin_size ),
int( ( $segment->end - 1 ) / $bin_size )
);
}
}
else {
$set->add_range( int( $a->pos / $bin_size ),
int( ( $a->calend - 1 ) / $bin_size ) );
}
}
# record
foreach my $pos ( $set->as_array ) {
$data->{f}->[ $pos - $data->{f_offset} ] += $score;
}
}
sub smart_stranded_pe {
my ( $f, $data, $score, $r ) = @_;
my $set = Set::IntSpan::Fast->new;
# determine strand based on the forward alignment
my ( $strand, $offset );
my $flag = $f->flag;
if ( $flag & 0x0040 ) {
# it's the first read, therefore fragment is forward
$strand = 'f';
$offset = 'f_offset';
}
elsif ( $flag & 0x0080 ) {
# it's the second read, therefore fragment is reverse
$strand = 'r';
$offset = 'r_offset';
}
# process both reads, adding to the integer set
foreach my $a ( $f, $r ) {
if ( $a->cigar_str =~ /N/ ) {
my $aw = $wrapper_ref->new( $a, $data->{sam} );
# check intron size from the cigar string
if ($max_intron) {
my $size = 1;
my $cigars = $aw->cigar_array;
foreach my $c ( @{$cigars} ) {
# each element is [operation, size]
$size = $c->[1] if ( $c->[0] eq 'N' and $c->[1] > $size );
}
return if $size > $max_intron; # exceed maximum intron size
}
# record
foreach my $segment ( $aw->get_SeqFeatures ) {
$set->add_range(
int( ( $segment->start - 1 ) / $bin_size ),
int( ( $segment->end - 1 ) / $bin_size )
);
}
}
else {
$set->add_range( int( $a->pos / $bin_size ),
int( ( $a->calend - 1 ) / $bin_size ) );
}
}
# record
foreach my $pos ( $set->as_array ) {
$data->{$strand}->[ $pos - $data->{$offset} ] += $score;
}
}
sub pe_ends {
my ( $a, $data, $score ) = @_;
# we always receive the forward read, never reverse, from pe_callback
# first position
$data->{f}->[ int( $a->pos / $bin_size ) - $data->{f_offset} ] += $score;
# second position
my $pos = int( ( $a->pos + $a->isize ) / $bin_size );
$data->{f}->[ $pos - $data->{f_offset} ] += $score;
}
sub pe_strand_ends {
my ( $a, $data, $score ) = @_;
# we always receive the forward read, never reverse, from pe_callback
# therefore the first and second read flag indicates orientation
# calculate both positions
my $flag = $a->flag;
if ( $flag & 0x0040 ) {
# first read implies forward orientation
# first position is forward
$data->{f}->[ int( $a->pos / $bin_size ) - $data->{f_offset} ] += $score;
# second position is reverse
my $pos = int( ( $a->pos + $a->isize ) / $bin_size );
$data->{r}->[ $pos - $data->{r_offset} ] += $score;
}
elsif ( $flag & 0x0080 ) {
# second read implies reverse orientation
# first position is reverse
$data->{r}->[ int( $a->pos / $bin_size ) - $data->{r_offset} ] += $score;
# second position is forward
my $pos = int( ( $a->pos + $a->isize ) / $bin_size );
$data->{f}->[ $pos - $data->{f_offset} ] += $score;
}
else {
die sprintf(
" Paired-end flags are set incorrectly; neither 0x040 or 0x080 are set for paired-end read %s.",
$a->query->name );
}
}
__END__
=head1 NAME
bam2wig.pl
A program to convert Bam alignments into a wig representation file.
=head1 SYNOPSIS
bam2wig.pl [--options...] E<lt>file.bamE<gt>
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 --exclude <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
=head1 OPTIONS
The command line flags and descriptions:
=head2 Input
=over 4
=item --in E<lt>filenameE<gt>
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.
=back
=head2 Reporting Options
=over 4
=item --start
Specify that the 5' position should be recorded in the wig file.
=item --mid
Specify that the midpoint of the alignment (single-end) or fragment
(paired-end) will be recorded in the wig file.
=item --span
Specify that the entire span of the alignment (single-end) or
fragment (paired-end) will be recorded in the wig file.
=item --extend
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.
=item --cspan
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.
=item --smartcov
Smart alignment coverage of paired-end alignments without
double-counting overlaps or recording gaps (intron splices).
=item --ends
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.
=item --coverage
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).
=item --position [start|mid|span|extend|cspan|coverage]
Legacy option for supporting previous versions of bam2wig.
=back
=head2 Alignment reporting options
=over 4
=item --splice
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.
=item --strand
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.
=item --flip
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.
=back
=head2 Paired-end alignments
=over 4
=item --pe
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.
=item --fastpe
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.
=item --minsize E<lt>integerE<gt>
Specify the minimum paired-end fragment size in bp to accept for recording.
Default is 30 bp.
=item --maxsize E<lt>integerE<gt>
Specify the maximum paired-end fragment size in bp to accept for recording.
Default is 600 bp.
=item --first
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.
=item --second
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.
=back
=head2 Alignment filtering options:
=over 4
=item --qual E<lt>integerE<gt>
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).
=item --nosecondary
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.
=item ---noduplicate
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.
=item --nosupplementary
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.
=item --chrskip E<lt>regexE<gt>
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'
=item --exclude E<lt>fileE<gt>
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 L<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.
=item --intron E<lt>integerE<gt>
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).
=back
=head2 Shift options
=over 4
=item --shift
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).
=item --shiftval E<lt>integerE<gt>
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).
=item --extval E<lt>integerE<gt>
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.
=item --chrom E<lt>integerE<gt>
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.
=item --minr E<lt>floatE<gt>
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.
=item --zmin E<lt>floatE<gt>
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.
=item --zmax E<lt>floatE<gt>
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.
=item --model
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.
=back
=head2 Score Options
=over 4
scripts/bam2wig.pl view on Meta::CPAN
Specify the output base filename. An appropriate extension will be
added automatically. By default it uses the base name of the
input file.
=item --bw
Specify whether or not the wig file should be further converted into
an indexed, compressed, binary BigWig file. The default is false.
=item --bwapp /path/to/wigToBigWig
Optionally specify the full path to the UCSC I<wigToBigWig> conversion
utility. The application path may be set in the F<.biotoolbox.cfg> file
or found in the default environment C<$PATH>, which makes this option
mostly unnecessary.
=item --gz
Specify whether (or not) the output text file should be compressed with
gzip. Disable with C<--nogz>. Does not apply to bigWig format.
=back
=head2 General options
=over 4
=item --cpu E<lt>integerE<gt>
Specify the number of parallel instances to run simultaneously. This requires
the installation of L<Parallel::ForkManager>. With support enabled, the
default is 4. Disable multi-threaded execution by setting to 1.
=item --temp E<lt>directoryE<gt>
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".
=item --verbose
Print extra informational statements during processing. The default is false.
=item --version
Print the version number.
=item --help
Display this POD documentation.
=back
=head1 DESCRIPTION
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.
=head1 RECOMMENDED SETTINGS
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.
=over
=item Straight coverage
To generate a straight-forward coverage map, similar to what most genome
browsers display when using a Bam file as source. B<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>
=item Smart paired-end coverage
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>
=item Single-end ChIP-Seq
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>
=item Paired-end ChIP-Seq
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>
=item Unstranded RNA-Seq
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>
=item Stranded, single-end RNA-Seq
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>
=item Paired-end RNA-Seq
Use the smart paired-end coverage mode to properly record paired-end
alignments with splice junctions.
bam2wig --smartcov --strand --rpm --in <bamfile>
=back
=head1 TEXT REPRESENTATION OF RECORDING ALIGNMENTS
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.
=over 4
=item Alignments
....>>>>>>.....................................<<<<<<.............
.....>>>>>>..................................<<<<<<...............
........>>>>>>.......................................<<<<<<.......
........>>>>>>.........................................<<<<<<.....
( run in 1.371 second using v1.01-cache-2.11-cpan-364913b4093 )