Bio-ToolBox
view release on metacpan or search on metacpan
in the Data object.
- Allow non-stranded values when parsing UCSC files, including bed files.
- Optimize scoring subroutines
- Remove legacy subroutines from utility module
- Include new test file for utility functions
- Numerous other small changes and fixes
v1.68
- Script bam2wig.pl script can now record both ends of paired-end
fragments, rather than faking it as single-end. Paired-end start
now respectes orientation. Added new option to only record either
first or second read in a pair. Added new option to ignore
zero intervals when writing bedGraph format. Changed multi-hit
scoring to preferentially use NH instead of IH.
- Scripts get_binned_data.pl and get_relative_data.pl now
can write out column names and associated datasets in separate
groups file for use in plotting. Also specify score decimal format.
- Script get_features.pl has new option to only keep features with
explicit tag value.
- High level ToolBox convenience function parse_file() now includes
from being written to the output file.
- Made the warning about unavailable files on the UCSC FTP server less
scary in the script ucsc_table2gff3.pl.
- Updated and clarified some script documentation.
v.1.10 (svn 559)
- Significantly improved performance when collecting data from Bam files
by using a low level API. Improvements of at least 2X may be realized.
- Significantly improved the performance of the bam2wig.pl script by at
least 2X. Added a new option of recording extended regions across the
predicted fragment based on empirically determined shift values.
Sampling to determine shift values has been increased. BedGraph files
are now written more efficiently. Maximum number of identical reads are
now enforced.
- Significantly improved the performance of the split_bam_by_isize.pl
script to increase speed by at least 2X. Added an option to skip
checking of mates. Improved reporting of results.
- Added a filter option to remove overlapping nucleosomes in script
verify_nucleosome_mapping.pl; also fixed bugs in reporting offset
distances and improved output reporting.
- Removed confusing separate scan and tag datasets required for script
lib/Bio/ToolBox.pm view on Meta::CPAN
=back
=head2 Data conversion
Convert from generic tables to specific bioinformatic file types.
=over 4
=item L<bam2wig.pl>
Generate read or fragment coverage or point data representations of alignments.
=item L<data2bed.pl>
Convert a table containing coordinates into a properly formatted BED file.
=item L<data2wig.pl>
Convert a table of coordinates and values into a properly formatted WIG file,
including bigWig.
lib/Bio/ToolBox/db_helper/alignment_callbacks.pm view on Meta::CPAN
return if $a->qual < $MAPQ;
my $s = $a->pos + 1;
my $e = $a->calend;
return
unless ( ( $s >= $data->{start} and $s <= $data->{stop} )
or ( $e >= $data->{start} and $e <= $data->{stop} ) );
if ( $flag & 0x10 ) {
# reversed
# since we're working with names, only record the 5' end of fragment
return if ( ( $flag & 0x1 ) and ( $flag & 0x2 ) );
# paired and proper
# not the end we're looking for
$data->{'index'}{$e} ||= [];
push @{ $data->{'index'}{$e} }, $a->qname;
}
else {
$data->{'index'}{$s} ||= [];
push @{ $data->{'index'}{$s} }, $a->qname;
lib/Bio/ToolBox/db_helper/bam.pm view on Meta::CPAN
return if $flag & 0xf04; # UNMAP,SECONDARY,QCFAIL,DUP,SUPPLEMENTARY
# check mapping quality based on global variable
return if $a->qual < $MAPQ;
# count properly paired alignments only once by skipping reverse alignment
if ( $flag & 0x01 and $flag & 0x2 ) {
return if ( $flag & 0x10 );
}
# count this fragment
${$n}++;
};
# prepare the counting subroutine
my $counter = sub {
my $tid = shift;
my $number = 0;
$sam->bam_index->fetch( $sam->bam, $tid, 0, $sam->target_len($tid),
$callback, \$number );
return $number;
lib/Bio/ToolBox/db_helper/bam.pm view on Meta::CPAN
The subroutine returns a hash or hash reference of the defined dataset values
found within the region of interest keyed by position. The feature midpoint
is used as the key position. When multiple features are found at the same
position, a simple mean (for length data methods) or sum
(for count methods) is returned. The ncount value type is not supported
with positioned scores.
=item sum_total_bam_alignments
This subroutine will sum the total number of mapped fragments
in a bam file. Properly paired alignments are counted once, and
supplemental, secondary, and marked duplicate alignments are ignored.
Alignments may be filtered by mapping quality. Use the
L<Bio::ToolBox::db_helper> C<use_minimum_mapq()> function to set
the minimum allowed mapping quality.
Pass the subroutine one or optionally two arguments in the following order.
=over 4
lib/Bio/ToolBox/db_helper/bigbed.pm view on Meta::CPAN
The subroutine returns a hash of the defined dataset values found within
the region of interest keyed by position. The feature midpoint is used
as the key position. When multiple features are found at the same
position, a simple mean (for score or length data methods) or sum
(for count methods) is returned.
=item sum_total_bigbed_features
This subroutine will sum the total number of bed features present in a
BigBed file. This may be useful, for example, in calculating fragments
(reads) per million mapped values when the bigbed file represents
sequence alignments.
Pass either the name of a bigBed file (F<.bb>), either local or remote, or an
opened BigBed database object. A scalar value of the total number of features
is returned.
=back
=head2 Data Collection Parameters Reference
lib/Bio/ToolBox/db_helper/hts.pm view on Meta::CPAN
return if $flag & 0xf04; # UNMAP,SECONDARY,QCFAIL,DUP,SUPPLEMENTARY
# check mapping quality based on global variable
return if $a->qual < $MAPQ;
# count properly paired alignments only once by skipping reverse alignment
if ( $flag & 0x01 and $flag & 0x2 ) {
return if ( $flag & 0x10 );
}
# count this fragment
${$n}++;
};
# prepare the counting subroutine
my $counter = sub {
my $tid = shift;
my $number = 0;
$bam->hts_index->fetch( $bam->hts_file, $tid, 0, $bam->target_len($tid),
$callback, \$number );
return $number;
lib/Bio/ToolBox/db_helper/hts.pm view on Meta::CPAN
The subroutine returns a hash or hash reference of the defined dataset values
found within the region of interest keyed by position. The feature midpoint
is used as the key position. When multiple features are found at the same
position, a simple mean (for length data methods) or sum
(for count methods) is returned. The ncount value type is not supported
with positioned scores.
=item sum_total_bam_alignments
This subroutine will sum the total number of mapped fragments
in a bam file. Properly paired alignments are counted once, and
supplemental, secondary, and marked duplicate alignments are ignored.
Alignments may be filtered by mapping quality. Use the
L<Bio::ToolBox::db_helper> C<use_minimum_mapq()> function to set
the minimum allowed mapping quality.
Pass the subroutine one or optionally two arguments in the following order.
=over 4
lib/Bio/ToolBox/utility.pm view on Meta::CPAN
# remove any known file extensions
$new_name =~
s/\. (?: bw | bam | bb | useq | bigwig | bigbed | g[tf]f3? | cram | wig | bdg | bedgraph ) (?:\.gz)? $//xi;
# remove common non-useful stuff
# trying to imagine all sorts of possible things
$new_name =~
s/[_\.\-] (?: sort | sorted | dedup | dedupe | deduplicated | rmdup | mkdup | markdup | dup | unique | filt | filtered ) \b //xgi;
$new_name =~
s/[_\.\-] (?: coverage | rpm | ext\d* | extend\d* | log2fe | log\d+ | qvalue | fragment | count | lambda_control | fe | fold.?enrichment | ratio | log\d*ratio ) \b //xgi;
}
return $new_name;
}
sub sane_chromo_sort {
my @chroms = @_;
return unless scalar @chroms;
# let's try and sort in some kind of rational order
my @numeric;
scripts/bam2wig.pl view on Meta::CPAN
### 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
scripts/bam2wig.pl view on Meta::CPAN
# 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 ) {
scripts/bam2wig.pl view on Meta::CPAN
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) {
scripts/bam2wig.pl view on Meta::CPAN
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
scripts/bam2wig.pl view on Meta::CPAN
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
scripts/bam2wig.pl view on Meta::CPAN
=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]
scripts/bam2wig.pl view on Meta::CPAN
=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.
scripts/bam2wig.pl view on Meta::CPAN
=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
scripts/bam2wig.pl view on Meta::CPAN
=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.+'
scripts/bam2wig.pl view on Meta::CPAN
=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).
scripts/bam2wig.pl view on Meta::CPAN
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.
scripts/bam2wig.pl view on Meta::CPAN
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.
scripts/correlate_position_data.pl view on Meta::CPAN
Two statistics may be calculated. First, it will calculate a a Pearson
linear correlation coefficient (r value) between the datasets (default).
Additionally, an ANOVA analysis may be performed between the datasets and
generate a P-value.
By default, the correlation is determined between the data points
collected over the entire length of the feature. Alternatively, a
radius and reference point (default is midpoint) may be provided
that sets the window for collecting scores and calculating a correlation.
In general, to ensure a more reliable Pearson value, fragment ChIP or
nucleosome coverage should be used rather than point (start or midpoint)
data, as it will give more reliable results. Fragment coverage is more
akin to smoothened data and gives better results than interpolated point
data.
Normalized read-depth data should be used when possible. If necessary,
Values can be normalized using one of two methods. The values may be
converted to rank positions (compare to Kendall's tau), or scaled such
that the absolute sum values are equal (for example, when working with
sequence tag read counts).
scripts/get_datasets.pl view on Meta::CPAN
# 3' end of reverse strand
$start = $row->start - $relative_stop;
$stop = $row->start - $relative_start;
}
# midpoint is not accepted
}
else {
# length doesn't meet minimum limit
# simply take the whole fragment
$start = $row->start;
$stop = $row->end;
}
# now collect score
my $score = $row->get_score(
'start' => $start,
'stop' => $stop,
'db' => $ddb,
'dataset' => $dataset,
scripts/get_datasets.pl view on Meta::CPAN
in a BED file of intervals.
get_datasets.pl --in input.bed --data scores.bw
=item Collect normalized counts
You want to collect normalized read counts from multiple Bam files
for each feature in a BED file. This will count alignment names (safe for
paired-end alignments) over the intervals, and transform to Fragments (Reads)
Per Million, a depth-normalizing function based on the total number of
fragments counted in each dataset.
get_datasets.pl --in input.bed --method ncount --fpkm region *.bam
=item Collect stranded RNASeq data
You have stranded RNASeq data, and you would like to determine the
expression level for all genes from an annotation file. Use the
C<ncount> method to count alignment names to avoid double counting
alignments split over multiple exons.
( run in 3.049 seconds using v1.01-cache-2.11-cpan-364913b4093 )