view release on metacpan or  search on metacpan

Bio/Search/Tiling/MapTiling.pm  view on Meta::CPAN

#
# BioPerl module for Bio::Search::Tiling::MapTiling
#
# Please direct questions and support issues to <bioperl-l@bioperl.org>
#
# Cared for by Mark A. Jensen <maj@fortinbras.us>
#
# Copyright Mark A. Jensen
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

=head1 NAME

Bio::Search::Tiling::MapTiling - An implementation of an HSP tiling
algorithm, with methods to obtain frequently-requested statistics

=head1 SYNOPSIS

 # get a BLAST $hit from somewhere, then
 $tiling = Bio::Search::Tiling::MapTiling->new($hit);

 # stats
 $numID = $tiling->identities();
 $numCons = $tiling->conserved();
 $query_length = $tiling->length('query');
 $subject_length = $tiling->length('subject'); # or...
 $subject_length = $tiling->length('hit');

 # get a visual on the coverage map
 print $tiling->coverage_map_as_text('query',$context,'LEGEND');

 # tilings
 $context = $tiling->_context( -type => 'subject', -strand=> 1, -frame=>1);
 @covering_hsps_for_subject = $tiling->next_tiling('subject',$context);
 $context = $tiling->_context( -type => 'query', -strand=> -1, -frame=>0);
 @covering_hsps_for_query   = $tiling->next_tiling('query', $context);

=head1 DESCRIPTION

Frequently, users want to use a set of high-scoring pairs (HSPs)
obtained from a BLAST or other search to assess the overall level of
identity, conservation, or coverage represented by matches between a
subject and a query sequence. Because a set of HSPs frequently
describes multiple overlapping sequence fragments, a simple summation of
statistics over the HSPs will generally overestimate those
statistics. To obtain an accurate estimate of global hit statistics, a
'tiling' of HSPs onto either the subject or the query sequence must be
performed, in order to properly correct for this. 

This module will execute a tiling algorithm on a given hit based on an
interval decomposition I'm calling the "coverage map". Internal object
methods compute the various statistics, which are then stored in
appropriately-named public object attributes. See
L<Bio::Search::Tiling::MapTileUtils> for more info on the algorithm. 

=head2 STRAND/FRAME CONTEXTS

In BLASTX, TBLASTN, and TBLASTX reports, strand and frame information
are reported for the query, subject, or query and subject,
respectively, for each HSP. Tilings for these sequence types are only
meaningful when they include HSPs in the same strand and frame, or 
"context". So, in these situations, the context must be specified
in the method calls or the methods will throw. 

Contexts are specified as strings: C<[ 'all' | [m|p][_|0|1|2] ]>, where
C<all> = all HSPs (will throw if context must be specified), C<m> = minus
strand, C<p> = plus strand, and C<_> = no frame info, C<0,1,2> = respective
(absolute) frame. The L</_make_context_key> method will convert a (strand,
frame) specification to a context string, e.g.:

    $context = $self->_context(-type=>'query', -strand=>-1, -frame=>-2);

returns C<m2>.

The contexts present among the HSPs in a hit are identified and stored
for convenience upon object construction. These are accessed off the
object with the L</contexts> method. If contexts don't apply for the
given report, this returns C<('all')>. 

=head1 TILED ALIGNMENTS

The experimental method L<ALIGNMENTS/get_tiled_alns> will use a tiling
to concatenate tiled hsps into a series of L<Bio::SimpleAlign>
objects:

 @alns = $tiling->get_tiled_alns($type, $context);

Each alignment contains two sequences with ids 'query' and 'subject',
and consists of a concatenation of tiling HSPs which overlap or are
directly adjacent. The alignment are returned in C<$type> sequence
order. When HSPs overlap, the alignment sequence is taken from the HSP
which comes first in the coverage map array.

The sequences in each alignment contain features (even though they are
L<Bio::LocatableSeq> objects) which map the original query/subject
coordinates to the new alignment sequence coordinates. You can
determine the original BLAST fragments this way:

 $aln = ($tiling->get_tiled_alns)[0];
 $qseq = $aln->get_seq_by_id('query');
 $hseq = $aln->get_seq_by_id('subject');
 foreach my $feat ($qseq->get_SeqFeatures) {
    $org_start = ($feat->get_tag_values('query_start'))[0];
    $org_end = ($feat->get_tag_values('query_end'))[0];
    # original fragment as represented in the tiled alignment:
    $org_fragment = $feat->seq;
 }
 foreach my $feat ($hseq->get_SeqFeatures) {
    $org_start = ($feat->get_tag_values('subject_start'))[0];
    $org_end = ($feat->get_tag_values('subject_end'))[0];
    # original fragment as represented in the tiled alignment:
    $org_fragment = $feat->seq;
 }

=head1 DESIGN NOTE

The major calculations are made just-in-time, and then memoized. So,
for example, for a given MapTiling object, a coverage map would
usually be calculated only once (for the query), and at most twice (if
the subject perspective is also desired), and then only when a
statistic is first accessed. Afterward, the map and/or any statistic
is read from storage. So feel free to call the statistic methods
frequently if it suits you.

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to
the Bioperl mailing list.  Your participation is much appreciated.

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

=head2 Support

Please direct usage questions or support issues to the mailing list:

I<bioperl-l@bioperl.org>

rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
of the bugs and their resolution. Bug reports can be submitted via
the web:

  https://github.com/bioperl/bioperl-live/issues

=head1 AUTHOR - Mark A. Jensen

Email maj -at- fortinbras -dot- us

=head1 APPENDIX

The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _

=cut

# Let the code begin...

package Bio::Search::Tiling::MapTiling;
use strict;
use warnings;

# Object preamble - inherits from Bio::Root::Root