Bio::Tools::Blast
Sbjct
Toolbar
Summary
Bio::Tools::Blast::Sbjct - Bioperl BLAST "Hit" object
Package variables
Privates (from "my" definitions)
$_prog = ''
$_signif_fmt = ''
Included modules
Inherit
Synopsis
The construction of HSP objects is handled by
Bio::Tools::Blast.
You should not need to use this package directly. See
_initialize()for a description of constructor parameters.
require Bio::Tools::Blast::Sbjct;
$hit = new Bio::Tools::Blast::Sbjct (-DATA =>\@hitData,
-PARENT =>$self,
-NAME =>5,
-RANK =>5,
-RANK_BY =>'order',
-MAKE =>'query' (or 'sbjct'),
-OVERLAP =>2,
-PROGRAM =>'TBLASTN'
);
@hitData includes the summary line for the hit as element [0], plus
all lines from the HSP alignment section of the BLAST report for
the present hit.
Description
The Bio::Tools::Blast::Sbjct.pm module encapsulates data and methods for
parsing and manipulating "hits" from a BLAST report.
This module is a utility module used by the Bio::Tools::Blast.pm
and is not intended for separate use.
In Blast lingo, the "sbjct" sequences are all the sequences
in a target database which were compared against a "query" sequence.
The terms "sbjct" and "hit" will be used interchangeably in this and related modules.
This module supports BLAST versions 1.x and 2.x, gapped and ungapped.
If a Blast hit has more than one HSP, the Bio::Tools::Blast::Sbjct.pm
object has the ability to merge overlapping HSPs into contiguous
blocks. This permits the Sbjct object to sum data across all HSPs
without counting data in the overlapping regions multiple times, which
would happen if data from each overlapping HSP are simply summed. HSP
tiling is performed automatically when methods of the Sbjct object
that rely on tiled data are invoked. These include
frac_identical(),
frac_conserved(),
gaps(),
frac_aligned_query(),
frac_aligned_hit(),
num_unaligned_query(),
num_unaligned_hit().
It also permits the assessment of an "ambiguous alignment" if the
query (or sbjct) sequences from different HSPs overlap. The existence
of an overlap could indicate a biologically interesting region in the
sequence, such as a repeated domain. The Sbjct object uses the
-OVERLAP parameter to determine when two sequences overlap; if this is
set to 2 -- the default -- then any two sbjct or query HSP sequences
must overlap by more than two residues to get merged into the same
contig and counted as an overlap. See the
BUGS section below for
"issues" with HSP tiling.
The results of the HSP tiling is reported with the following ambiguity codes:
'q' = Query sequence contains multiple sub-sequences matching
a single region in the sbjct sequence.
's' = Sbjct sequence contains multiple sub-sequences matching
a single region in the query sequence.
'qs' = Both query and sbjct sequences contain more than one
sub-sequence with similarity to the other sequence.
For addition information about ambiguous BLAST alignments, see
_tile_hsps() and
http://www-genome.stanford.edu/Sacch3D/help/ambig_aln.html
Methods
Methods description
Usage : n/a; called automatically during object construction.
Purpose : Builds HSP contigs for a given BLAST hit.
: Utility method called by _tile_hsps()
Returns :
Argument :
Throws : Exceptions propagated from Bio::Tools::Blast::HSP::matches()
: for invalid sub-sequence ranges.
Status : Experimental
Comments : This method does not currently support gapped alignments.
: Also, it does not keep track of the number of HSPs that
: overlap within the amount specified by overlap().
: This will lead to significant tracking errors for large
: overlap values.
See Also : overlap(), _tile_hsps(), Bio::Tools::Blast::HSP::matches |
Usage : n/a; called automatically by display()
Purpose : Display information about each HSP in the current BLAST hit.
Example : n/a
Argument : one argument = filehandle object.
Returns : printf call.
Status : Experimental
See Also : display(), Bio::Tools::Blast::HSP::display() |
Usage : n/a; called automatically by display()
Purpose : Display information about Bio::Tools::Blast.pm data members.
: Not tab-delimited.
: Prints the rank, name, database, score, p, n, length
: of the hit sequence, length of the aligned region,
: fraction identical, fraction conserved, and the fraction aligned
: for both the query and hit sequences.
Example : n/a
Argument : one argument = filehandle object.
Returns : printf call.
Status : Deprecated, Buggy.
: Use Blast::table() or Blast::table_tiled() instead.
See Also : display() |
Usage : n/a; automatically called by Bio::Root::Object::new()
: Bio::Tools::Blast::Sbjct.pm objects are constructed
: automatically by Bio::Tools::Blast.pm, so there is no need
: for direct consumption.
Purpose : Initializes key varaiables and calls methods to parse a single Blast hit.
: Constructs Bio::Tools::Blast::HSP.pm objects for each
: high-scoring segment pair (HSP).
: Calls superclass constructor first (Bio::Root::Object.pm).
Returns : n/a
Argument : Named Parameters passed from new()
: ALL TAGS MUST BE UPPERCASE (does not call _rearrange()).
: -DATA => array reference holding all data for a single hit.
: The first element should hold the description
: line (from the desctiption section at the top of
: the report), remaining lines should hold all lines
: within the HSP alignment listing section of report.
: -PARENT => object reference to a Bio::Tools::Blast.pm object.
: -NAME => string (typically same as -RANK, just a temporary
: name to use until the actual name of hit is parsed),
: -RANK => integer,
: -RANK_BY => 'order',
: -OVERLAP => integer (maximum overlap between adjacent
: HSPs when tiling)
: -PROGRAM => string (type of Blast blastp, blastn, etc).
See Also : _set_id(), _set_hsps(), _tile_hsps(), Bio::Root::Object::new(), Bio::Tools::Blast::_set_hit_db() |
Usage : n/a; called automatically by _set_hsps()
Purpose : Sets the description of the hit sequence.
: For sequence without descriptions, sets description to "-".
Argument : Array containing description (multiple lines).
Comments : _set_hsps() calls this method with the data from the
: HSP alignment listing, which contains the complete description.
: (Formerly, this was called from the _set_desc_data() method initially.)
See Also : _set_hsps() |
Usage : n/a; called automatically during object construction.
Purpose : Creates HSP.pm objects for each HSP in a BLAST hit alignment.
: Also collects the full description of the hit from the
: HSP alignment section.
Returns : n/a
Argument : List of strings containing raw BLAST report data for
: a single hit's HSP alignment data.
Throws : Warnings for each HSP.pm object that fails to be constructed.
: Exception if no HSP.pm objects can be constructed.
: Exception if can't parse length data for hit sequence.
Comments : Requires Bio::Tools::Blast::HSP.pm.
: Sets the description using the full string present in
: the alignment data.
: Also sets Expect and P-values for the Sbjct object by
: copying from the HSP object.
: While this sacrifices some memory efficiency, it
: improves access speed for these critical data.
See Also : _initialize(), _set_desc() |
Usage : n/a; automatically called by _initialize()
Purpose : Sets the name of the Sbjct sequence from the BLAST summary line.
: The identifier is assumed to be the first
: chunk of non-whitespace characters in the description line
: Does not assume any semantics in the structure of the identifier
: (Formerly, this method attempted to extract database name from
: the seq identifiers, but this was prone to break).
Returns : n/a
Argument : String containing description line of the hit from Blast report
: or first line of an alignment section.
Throws : Warning if cannot locate sequence ID.
See Also : _initialize(), Bio::Tools::Blast::database() |
Usage : n/a; called automatically during object construction or
: as needed by methods that rely on having tiled data.
Purpose : Collect statistics about the aligned sequences in a set of HSPs.
: Calculates the following data across all HSPs:
: -- total alignment length
: -- total identical residues
: -- total conserved residues
Returns : n/a
Argument : n/a
Throws : n/a
Status : Experimental
Comments :
: This method performs more careful summing of data across
: all HSPs in the Sbjct object. Simply summing the data from all HSPs
: will overestimate the actual length of the alignment if there is
: overlap between different HSPs (often the case).
: The strategy is to tile the HSPs and sum over the
: contigs, collecting data separately from overlapping and
: non-overlapping regions of each HSP. To facilitate this, the
: HSP.pm object now permits extraction of data from sub-sections
: of an HSP.
:
: Additional useful information is collected from the results
: of the tiling. It is possible that sub-sequences in
: different HSPs will overlap significantly. In this case, it
: is impossible to create a single unambiguous alignment by
: concatenating the HSPs. The ambiguity may indicate the
: presence of multiple, similar domains in one or both of the
: aligned sequences. This ambiguity is recorded using the
: ambiguous_aln() method.
:
: This method does not attempt to discern biologically
: significant vs. insignificant overlaps. The allowable amount of
: overlap can be set with the overlap() method or with the -OVERLAP
: parameter used when constructing the Blast & Sbjct objects.
:
: For a given hit, both the query and the sbjct sequences are
: tiled independently.
:
: -- If only query sequence HSPs overlap,
: this may suggest multiple domains in the sbjct.
: -- If only sbjct sequence HSPs overlap,
: this may suggest multiple domains in the query.
: -- If both query & sbjct sequence HSPs overlap,
: this suggests multiple domains in both.
: -- If neither query & sbjct sequence HSPs overlap,
: this suggests either no multiple domains in either
: sequence OR that both sequences have the same
: distribution of multiple similar domains.
:
: This method can deal with the special case of when multiple
: HSPs exactly overlap.
:
: Efficiency concerns:
: Speed will be an issue for sequences with numerous HSPs.
:
Bugs : Currently, _tile_hsps() does not properly account for
: the number of non-tiled but overlapping HSPs, which becomes a problem
: as overlap() grows. Large values overlap() may thus lead to
: incorrect statistics for some hits. For best results, keep overlap()
: below 5 (DEFAULT IS 2). For more about this, see the "HSP Tiling and
: Ambiguous Alignments" section.
See Also : _adjust_contigs(), ambiguous_aln(), overlap(), frac_identical(), frac_conserved(), frac_aligned_query(), frac_aligned_hit(), num_unaligned_query(), num_unaligned_hit(), HSP Tiling and Ambiguous Alignments |
Usage : $ambig_code = $sbjct_object->ambiguous_aln();
Purpose : Sets/Gets ambiguity code data member.
Example : (see usage)
Returns : String = 'q', 's', 'qs', '-'
: 'q' = query sequence contains overlapping sub-sequences
: while sbjct does not.
: 's' = sbjct sequence contains overlapping sub-sequences
: while query does not.
: 'qs' = query and sbjct sequence contains overlapping sub-sequences
: relative to each other.
: '-' = query and sbjct sequence do not contains multiple domains
: relative to each other OR both contain the same distribution
: of similar domains.
Argument : n/a
Throws : n/a
Status : Experimental
See Also : _tile_hsps(), HSP Tiling and Ambiguous Alignments |
Usage : $sbjct_object->bits();
Purpose : Gets the BLAST bit score of the best HSP for the current Blast hit.
Example : $bits = $sbjct_object->bits();
Returns : Integer
Argument : n/a
Throws : Exception if bit score is not set.
Comments : For BLAST1, the non-bit score is listed in the summary line.
See Also : score() |
Usage : $sbjct_object->database();
Purpose : Get the name of the database for the hit sequence.
Example : (see usage)
Returns : String
Argument : n/a
Throws : n/a
Status : Deprecated. Use Bio::Tools::Blast::database()
Extracting database name from the seq identifier is error prone.
Comments : Database id should be the same for all hits in a given
: BLAST report, however, they do not always have the same
: name as the database name extraced by the Blast.pm object.
: The Sbjct.pm database id is obtained from the summary line. |
Usage : $sbjct_object->desc( [integer] );
Purpose : Get the description for the given BLAST hit.
Example : (see usage)
Returns : String
Argument : Integer (optional) indicating the desired length of the
: description string to be returned.
Throws : n/a
See Also : _set_desc() |
Usage : $sbjct_object->display( %named_parameters );
Purpose : Display information about Bio::Tools::Blast::Sbjct.pm data members
Example : $object->display(-SHOW=>'stats');
Argument : Named parameters: -SHOW => 'hsp',
: -WHERE => filehandle (default = STDOUT)
Returns : n/a
Status : Deprecated, Buggy.
: Use Blast::table() or Blast::table_tiled() instead.
See Also : _display_stats(), _display_hsps(), Bio::Root::Object::display |
Usage : $sbjct->end( [seq_type] );
Purpose : Gets the end coordinate for the query, sbjct, or both sequences
: in the Sbjct object. If there is more than one HSP, the largest end
: value of all HSPs is returned.
Example : $qend = $sbjct->end('query');
: $send = $sbjct->end('sbjct');
: ($qend, $send) = $sbjct->end();
Returns : scalar context: integer
: array context without args: list of two integers (queryEnd, sbjctEnd)
: Array context can be "induced" by providing an argument of 'list' or 'array'.
Argument : In scalar context: seq_type = 'query' or 'sbjct'
: (case insensitive). If not supplied, 'query' is used.
Throws : n/a
Comments : This method requires that all HSPs be tiled. If there is more than one
: HSP and they have not already been tiled, they will be tiled first.
: Remember that the start and end coordinates of all HSPs are
: normalized so that start < end. Strand information can only be
: obtained on an HSP-by-HSP basis by calling $hsp->strand().
See Also : start(), range(), HSP Tiling and Ambiguous Alignments, Bio::Tools::Blast::HSP::end() |
Usage : $sbjct_object->expect( [format] );
Purpose : Get the Expect value for the given BLAST hit.
Example : $e = $sbjct->expect;
: $e = $sbjct->expect('exp'); # get exponent only.
: ($num, $exp) = $sbjct->expect('parts'); # split sci notation into parts
Returns : Float or scientific notation number (the raw expect value, DEFAULT).
: Integer if format == 'exp' (the magnitude of the base 10 exponent).
: 2-element list (float, int) if format == 'parts' and Expect
: is in scientific notation (see Comments).
Argument : format: string of 'raw' | 'exp' | 'parts'
: 'raw' returns value given in report. Default. (1.2e-34)
: 'exp' returns exponent value only (34)
: 'parts' returns the decimal and exponent as a
: 2-element list (1.2, -34) (see Comments).
Throws : Exception if the Expect value is not defined.
Comments : Using the 'parts' argument is not recommended since it will not
: work as expected if the expect value is not in scientific notation.
: That is, floats are not converted into sci notation before
: splitting into parts.
See Also : p(), signif(), get_exponent() |
Usage : $sbjct_object->frac_aligned_hit();
Purpose : Get the fraction of the hit (sbjct) sequence which has been aligned
: across all HSPs (not including intervals between non-overlapping
: HSPs).
Example : $frac_alnq = $sbjct_object->frac_aligned_hit();
Returns : Float (2-decimal precision, e.g., 0.75).
Argument : n/a
Throws : n/a
Comments : If you need data for each HSP, use hsps() and then interate
: through the HSP objects.
: To compute the fraction aligned, the logical length of the sbjct
: sequence is used, meaning that for TBLAST[NX] reports, the
: full length of the sbjct sequence is converted into amino acids
: by dividing by 3. This is necessary because of the way
: the lengths of aligned sequences are computed.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
See Also : frac_aligned_query(), matches(), _tile_hsps(), logical_length(), length_aln() |
Usage : $sbjct_object->frac_aligned_query();
Purpose : Get the fraction of the query sequence which has been aligned
: across all HSPs (not including intervals between non-overlapping
: HSPs).
Example : $frac_alnq = $sbjct_object->frac_aligned_query();
Returns : Float (2-decimal precision, e.g., 0.75).
Argument : n/a
Throws : n/a
Comments : If you need data for each HSP, use hsps() and then interate
: through the HSP objects.
: To compute the fraction aligned, the logical length of the query
: sequence is used, meaning that for [T]BLASTX reports, the
: full length of the query sequence is converted into amino acids
: by dividing by 3. This is necessary because of the way
: the lengths of aligned sequences are computed.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
See Also : frac_aligned_hit(), _tile_hsps(), logical_length(), length_aln() |
Usage : $sbjct_object->frac_conserved( [seq_type] );
Purpose : Get the overall fraction of conserved positions across all HSPs.
: The number refers to only the aligned regions and does not
: account for unaligned regions in between the HSPs, if any.
Example : $frac_cons = $sbjct_object->frac_conserved('sbjct');
Returns : Float (2-decimal precision, e.g., 0.75).
Argument : seq_type: 'query' | 'sbjct' | 'total'
: default = 'total' (but see comments below).
Throws : n/a
Comments : Different versions of Blast report different values for the total
: length of the alignment. This is the number reported in the
: denominators in the stats section:
: "Identical = 34/120 Positives = 67/120".
: BLAST-GP uses the total length of the alignment (with gaps)
: WU-BLAST uses the length of the query sequence (without gaps).
: Therefore, when called without an argument or an argument of 'total',
: this method will report different values depending on the
: version of BLAST used.
:
: To get the fraction conserved among only the aligned residues,
: ignoring the gaps, call this method with an argument of 'query'
: or 'sbjct'.
:
: If you need data for each HSP, use hsps() and then interate
: through the HSP objects.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
See Also : frac_identical(), matches(), _tile_hsps() |
Usage : $sbjct_object->frac_identical( [seq_type] );
Purpose : Get the overall fraction of identical positions across all HSPs.
: The number refers to only the aligned regions and does not
: account for unaligned regions in between the HSPs, if any.
Example : $frac_iden = $sbjct_object->frac_identical('query');
Returns : Float (2-decimal precision, e.g., 0.75).
Argument : seq_type: 'query' | 'sbjct' | 'total'
: default = 'total' (but see comments below).
Throws : n/a
Comments : Different versions of Blast report different values for the total
: length of the alignment. This is the number reported in the
: denominators in the stats section:
: "Identical = 34/120 Positives = 67/120".
: BLAST-GP uses the total length of the alignment (with gaps)
: WU-BLAST uses the length of the query sequence (without gaps).
: Therefore, when called without an argument or an argument of 'total',
: this method will report different values depending on the
: version of BLAST used.
:
: To get the fraction identical among only the aligned residues,
: ignoring the gaps, call this method with an argument of 'query'
: or 'sbjct'.
:
: If you need data for each HSP, use hsps() and then iterate
: through the HSP objects.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
See Also : frac_conserved(), frac_aligned_query(), matches(), _tile_hsps() |
Usage : $sbjct_object->frame();
Purpose : Gets the reading frame for the hit sequence (TBLASTN/X only).
Example : $frame = $sbjct_object->frame();
Returns : Integer (-3 .. +3).
Argument : n/a
Throws : Exception if HSPs have not been set (BLAST2 reports).
See Also : hsps() |
Usage : $sbjct_object->gaps( [seq_type] );
Purpose : Get the number of gaps in the aligned query, sbjct, or both sequences.
: Data is summed across all HSPs.
Example : $qgaps = $sbjct_object->gaps('query');
: $sgaps = $sbjct_object->gaps('sbjct');
: $tgaps = $sbjct_object->gaps(); # default = total (query + sbjct)
Returns : scalar context: integer
: array context without args: two-element list of integers
: (queryGaps, sbjctGaps)
: Array context can be "induced" by providing an argument of 'list' or 'array'.
Argument : seq_type: 'query' | 'sbjct' | 'total' | 'list' (default = 'total')
Throws : n/a
Comments : If you need data for each HSP, use hsps() and then interate
: through each HSP object.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
: Not relying on wantarray since that will fail in situations
: such as printf "%d", $hit->gaps() in which you might expect to
: be printing the total gaps, but evaluates to array context.
See Also : length_aln() |
Usage : &get_exponent( number );
Purpose : Determines the power of 10 exponent of an integer, float,
: or scientific notation number.
Example : &get_exponent("4.0e-206");
: &get_exponent("0.00032");
: &get_exponent("10.");
: &get_exponent("1000.0");
: &get_exponent("e+83");
Argument : Float, Integer, or scientific notation number
Returns : Integer representing the exponent part of the number (+ or -).
: If argument == 0 (zero), return value is "-999".
Comments : Exponents are rounded up (less negative) if the mantissa is >= 5.
: Exponents are rounded down (more negative) if the mantissa is <= -5.
: This method probably belongs in a more general utility class. |
Usage : $data = $sbjct_object->homo_data( %named_params );
Purpose : Gets specific similarity data about all HSPs.
Returns : String
Argument : named parameters forwarded to Bio::Tools::Blast::HSP::homol_data().
Throws : n/a
Status : Experimental
Comments : This is an experimental method used for obtaining an
: indication of:
: 1) how many HSPs are in a Blast alignment
: 2) how strong the similarity is between sequences in the HSP
: 3) the endpoints of the alignment (sequence monomer numbers)
: "Homology data" for each HSP is in the format:
: "<integer> <start> <stop>"
: Data for different HSPs are tab-delimited.
See Also : Bio::Tools::Blast::homol_data(), Bio::Tools::Blast::HSP::homol_data() |
Usage : $sbjct_object->hsp( [string] );
Purpose : Get a single HSP.pm object for the present Sbjct.pm object.
Example : $hspObj = $sbjct_object->hsp; # same as 'best'
: $hspObj = $sbjct_object->hsp('best');
: $hspObj = $sbjct_object->hsp('worst');
Returns : Object reference for a Bio::Tools::Blast::HSP.pm object.
Argument : String (or no argument).
: No argument (default) = highest scoring HSP (same as 'best').
: 'best' or 'first' = highest scoring HSP.
: 'worst' or 'last' = lowest scoring HSP.
Throws : Exception if the HSPs have not been collected.
: Exception if an unrecognized argument is used.
See Also : hsps(), num_hsps(), _set_hsps() |
Usage : $sbjct_object->hsps();
Purpose : Get a list containing all HSP objects.
: Get the numbers of HSPs for the current hit.
Example : @hsps = $sbjct_object->hsps();
: $num = $sbjct_object->hsps(); # alternatively, use num_hsps()
Returns : Array context : list of Bio::Tools::Blast::HSP.pm objects.
: Scalar context: integer (number of HSPs).
: (Equivalent to num_hsps()).
Argument : n/a. Relies on wantarray
Throws : Exception if the HSPs have not been collected.
See Also : hsp(), num_hsps(), _set_hsps() |
Usage : $sbjct_object->is_signif();
Purpose : Determine if the given BLAST hit is significant.
Example :
Returns : Boolean
Argument : n/a
Throws : n/a
Comments : Uses criteria defined in the parent Blast.pm object
: to assess significance. Currently, only relies on
: P-value and length criteria.
: This mehtod is largely obsolete since are hits are now by
: definition significant. |
Usage : $sbjct_object->length();
Purpose : Get the total length of the hit sequence.
Example : $len = $sbjct_object->length();
Returns : Integer
Argument : n/a
Throws : n/a
Comments : Developer note: when using the built-in length function within
: this module, call it as CORE::length().
See Also : logical_length(), length_aln() |
Usage : $sbjct_object->length_aln( [seq_type] );
Purpose : Get the total length of the aligned region for query or sbjct seq.
: This number will include all HSPs
Example : $len = $sbjct_object->length_aln(); # default = query
: $lenAln = $sbjct_object->length_aln('query');
Returns : Integer
Argument : seq_Type = 'query' | 'sbjct' (Default = 'query')
Throws : Exception if the argument is not recognized.
Comments : This method will report the logical length of the alignment,
: meaning that for TBLAST[NX] reports, the length is reported
: using amino acid coordinate space (i.e., nucleotides / 3).
:
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
: If you don't want the tiled data, iterate through each HSP
: calling length() on each (use hsps() to get the HSPs).
See Also : length(), frac_aligned_query(), frac_aligned_hit(), gaps(), _tile_hsps(), Bio::Tools::Blast::HSP::length() |
Usage : $sbjct_object->logical_length( [seq_type] );
: (mostly intended for internal use).
Purpose : Get the logical length of the hit sequence.
: If the Blast is a TBLASTN or TBLASTX, the returned length
: is the length of the would-be amino acid sequence (length/3).
: For all other BLAST flavors, this function is the same as length().
Example : $len = $sbjct_object->logical_length();
Returns : Integer
Argument : seq_type = 'query' or 'sbjct' (default = 'query')
Throws : n/a
Comments : This is important for functions like frac_aligned_query()
: which need to operate in amino acid coordinate space when dealing
: with [T]BLAST[NX] type reports.
See Also : length(), frac_aligned_query(), frac_aligned_hit() |
Usage : $sbjct_object->matches( [class] );
Purpose : Get the total number of identical or conserved matches
: (or both) across all HSPs.
: (Note: 'conservative' matches are indicated as 'positives'
: in the Blast report.)
Example : ($id,$cons) = $sbjct_object->matches(); # no argument
: $id = $sbjct_object->matches('id');
: $cons = $sbjct_object->matches('cons');
Returns : Integer or a 2-element array of integers
Argument : class = 'id' | 'cons' OR none.
: If no argument is provided, both identical and conservative
: numbers are returned in a two element list.
: (Other terms can be used to refer to the conservative
: matches, e.g., 'positive'. All that is checked is whether or
: not the supplied string starts with 'id'. If not, the
: conservative matches are returned.)
Throws : Exception if the requested data cannot be obtained.
Comments : If you need data for each HSP, use hsps() and then interate
: through the HSP objects.
: Does not rely on wantarray to return a list. Only checks for
: the presence of an argument (no arg = return list).
See Also : Bio::Tools::Blast::HSP::matches(), hsps() |
Usage : $sbjct_object->n();
Purpose : Gets the N number for the current Blast hit.
: This is the number of HSPs in the set which was ascribed
: the lowest P-value (listed on the description line).
: This number is not the same as the total number of HSPs.
: To get the total number of HSPs, use num_hsps().
Example : $n = $sbjct_object->n();
Returns : Integer
Argument : n/a
Throws : Exception if HSPs have not been set (BLAST2 reports).
Comments : Note that the N parameter is not reported in gapped BLAST2.
: Calling n() on such reports will result in a call to num_hsps().
: The num_hsps() method will count the actual number of
: HSPs in the alignment listing, which may exceed N in
: some cases.
See Also : num_hsps() |
Usage : $sbjct_object->num_hsps();
Purpose : Get the number of HSPs for the present Blast hit.
Example : $nhsps = $sbjct_object->num_hsps();
Returns : Integer
Argument : n/a
Throws : Exception if the HSPs have not been collected.
See Also : hsps() |
Usage : $sbjct_object->num_unaligned_hit();
Purpose : Get the number of the unaligned residues in the hit sequence.
: Sums across all all HSPs.
Example : $num_unaln = $sbjct_object->num_unaligned_hit();
Returns : Integer
Argument : n/a
Throws : n/a
Comments : See notes regarding logical lengths in the comments for frac_aligned_hit().
: They apply here as well.
: If you need data for each HSP, use hsps() and then interate
: through the HSP objects.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
See Also : num_unaligned_query(), _tile_hsps(), frac_aligned_hit() |
Usage : $sbjct_object->num_unaligned_query();
Purpose : Get the number of the unaligned residues in the query sequence.
: Sums across all all HSPs.
Example : $num_unaln = $sbjct_object->num_unaligned_query();
Returns : Integer
Argument : n/a
Throws : n/a
Comments : See notes regarding logical lengths in the comments for frac_aligned_query().
: They apply here as well.
: If you need data for each HSP, use hsps() and then interate
: through the HSP objects.
: This method requires that all HSPs be tiled. If they have not
: already been tiled, they will be tiled first.
See Also : num_unaligned_hit(), _tile_hsps(), frac_aligned_query() |
Usage : $blast_object->overlap( [integer] );
Purpose : Gets/Sets the allowable amount overlap between different HSP sequences.
Example : $blast_object->overlap(5);
: $overlap = $blast_object->overlap();
Returns : Integer.
Argument : integer.
Throws : n/a
Status : Experimental
Comments : Any two HSPs whose sequences overlap by less than or equal
: to the overlap() number of resides will be considered separate HSPs
: and will not get tiled by _adjust_contigs().
See Also : _adjust_contigs(), BUGS |
Usage : $sbjct_object->p( [format] );
Purpose : Get the P-value for the given BLAST hit.
: (Note that P-values are not provided with NCBI Blast2 reports).
Example : $p = $sbjct->p;
: $p = $sbjct->p('exp'); # get exponent only.
: ($num, $exp) = $sbjct->p('parts'); # split sci notation into parts
Returns : Float or scientific notation number (the raw P-value, DEFAULT).
: Integer if format == 'exp' (the magnitude of the base 10 exponent).
: 2-element list (float, int) if format == 'parts' and P-value
: is in scientific notation (See Comments).
Argument : format: string of 'raw' | 'exp' | 'parts'
: 'raw' returns value given in report. Default. (1.2e-34)
: 'exp' returns exponent value only (34)
: 'parts' returns the decimal and exponent as a
: 2-element list (1.2, -34) (See Comments).
Throws : Exception if the P-value is not defined, which will occur
: with any NCBI Blast2 report.
Comments : Using the 'parts' argument is not recommended since it will not
: work as expected if the P-value is not in scientific notation.
: That is, floats are not converted into sci notation before
: splitting into parts.
See Also : expect(), signif(), get_exponent() |
Usage : $sbjct->range( [seq_type] );
Purpose : Gets the (start, end) coordinates for the query or sbjct sequence
: in the HSP alignment.
Example : ($qbeg, $qend) = $sbjct->range('query');
: ($sbeg, $send) = $sbjct->range('sbjct');
Returns : Two-element array of integers
Argument : seq_type = string, 'query' or 'sbjct' (default = 'query')
: (case insensitive).
Throws : n/a
See Also : start(), end() |
Usage : $sbjct->rank( integer or string );
Purpose : Sets/Gets the rank of the current Sbjct object relative to
: other Sbjct objects managed by a given Blast object.
Example : $sbjct->rank(1);
Returns : The current rank value.
Argument : Integer or string to be used for ranking the hit
: relative to other hits.
Throws : n/a
Comments : The rank usually corresponds to the order the listing
: of hits in the BLAST report from lowest to highest p-value.
: Rank need not be restricted to this value.
: rank() may be provided by a delegated or inherited
: iterator class in the future (such as Bio::Root::Vector.pm). |
Usage : $sbjct_object->score();
Purpose : Gets the BLAST score of the best HSP for the current Blast hit.
Example : $score = $sbjct_object->score();
Returns : Integer
Argument : n/a
Throws : n/a
See Also : bits() |
Usage : $hit->seq_inds( seq_type, class, collapse );
Purpose : Get a list of residue positions (indices) across all HSPs
: for identical or conserved residues in the query or sbjct sequence.
Example : @ind = $hit->seq_inds('query', 'identical');
: @ind = $hit->seq_inds('sbjct', 'conserved');
: @ind = $hit->seq_inds('sbjct', 'conserved', 1);
Returns : Array of integers
: May include ranges if collapse is non-zero.
Argument : seq_type = 'query' or 'sbjct' (default = query)
: class = 'identical' or 'conserved' (default = identical)
: (can be shortened to 'id' or 'cons')
: (actually, anything not 'id' will evaluate to 'conserved').
: collapse = boolean, if non-zero, consecutive positions are merged
: using a range notation, e.g., "1 2 3 4 5 7 9 10 11"
: collapses to "1-5 7 9-11". This is useful for
: consolidating long lists. Default = no collapse.
Throws : n/a.
See Also : Bio::Tools::Blast::HSP::seq_inds() |
Usage : $sbjct_object->signif( [format] );
Purpose : Get the P or Expect value for the given BLAST hit.
: The value returned is the one which is reported in the description
: section of the Blast report. For Blast1 and WU-Blast2, this
: is a P-value, for Blast2, it is an Expect value.
Example : $obj->signif() # returns 1.3e-34
: $obj->signif('exp') # returns -34
: $obj->signif('parts') # returns (1.3, -34)
Returns : Float or scientific notation number (the raw P/Expect value, DEFAULT).
: Integer if format == 'exp' (the magnitude of the base 10 exponent).
: 2-element list (float, int) if format == 'parts' and P/Expect value
: is in scientific notation (see Comments).
Argument : format: string of 'raw' | 'exp' | 'parts'
: 'raw' returns value given in report. Default. (1.2e-34)
: 'exp' returns exponent value only (34)
: 'parts' returns the decimal and exponent as a
: 2-element list (1.2, -34) (see Comments).
Throws : n/a
Status : Deprecated. Use p() or expect().
Comments : The signif() method provides a way to deal with the fact that
: Blast1 and Blast2 formats differ in what is reported in the
: description lines of each hit in the Blast report. The signif()
: method frees any client code from having to know if this is a P-value
: or an Expect value, making it easier to write code that can process
: both Blast1 and Blast2 reports. This is not necessarily a good thing, since
: one should always know when one is working with P-values or
: Expect values (hence the deprecated status).
: Use of expect() is recommended since all hits will have an Expect value.
:
: Using the 'parts' argument is not recommended since it will not
: work as expected if the expect value is not in scientific notation.
: That is, floats are not converted into sci notation before
: splitting into parts.
See Also : p(), expect(), get_exponent() |
Usage : $sbjct->start( [seq_type] );
Purpose : Gets the start coordinate for the query, sbjct, or both sequences
: in the Sbjct object. If there is more than one HSP, the lowest start
: value of all HSPs is returned.
Example : $qbeg = $sbjct->start('query');
: $sbeg = $sbjct->start('sbjct');
: ($qbeg, $sbeg) = $sbjct->start();
Returns : scalar context: integer
: array context without args: list of two integers (queryStart, sbjctStart)
: Array context can be "induced" by providing an argument of 'list' or 'array'.
Argument : In scalar context: seq_type = 'query' or 'sbjct'
: (case insensitive). If not supplied, 'query' is used.
Throws : n/a
Comments : This method requires that all HSPs be tiled. If there is more than one
: HSP and they have not already been tiled, they will be tiled first.
: Remember that the start and end coordinates of all HSPs are
: normalized so that start < end. Strand information can only be
: obtained on an HSP-by-HSP basis by calling $hsp->strand().
See Also : end(), range(), HSP Tiling and Ambiguous Alignments, Bio::Tools::Blast::HSP::start() |
Methods code
sub _adjust_contigs
{
my ($seqType, $hsp, $start, $stop, $contigs_ref, $max_overlap) = @_;
my $overlap = 0;
my ($numID, $numCons);
foreach(@$contigs_ref) {
if($start >= $_->{'start'} and $stop <= $_->{'stop'}) {
$overlap = 1;
next;
}
if($start < $_->{'start'} and $stop > ($_->{'start'} + $max_overlap)) {
eval {
($numID, $numCons) = $hsp->matches(-SEQ =>$seqType,
-START =>$start,
-STOP =>$_->{'start'}-1);
};
if($@) { warn "\a\n$@\n"; }
else {
$_->{'start'} = $start;
$_->{'iden'} += $numID;
$_->{'cons'} += $numCons;
$overlap = 1;
}
}
if($stop > $_->{'stop'} and $start < ($_->{'stop'} - $max_overlap)) {
eval {
($numID,$numCons) = $hsp->matches(-SEQ =>$seqType,
-START =>$_->{'stop'},
-STOP =>$stop);
};
if($@) { warn "\a\n$@\n"; }
else {
$_->{'stop'} = $stop;
$_->{'iden'} += $numID;
$_->{'cons'} += $numCons;
$overlap = 1;
}
}
$overlap && do {
last;
};
}
!$overlap && do {
($numID,$numCons) = $hsp->matches(-SEQ=>$seqType);
push @$contigs_ref, {'start'=>$start, 'stop'=>$stop,
'iden'=>$numID, 'cons'=>$numCons };
};
$overlap;
} | sub _display_hsps
{
my( $self, %param) = @_;
my $OUT = $self->fh();
my $hspCount = 0;
my $reply = undef;
not defined $self->{'_hsps'} and do{ print $OUT "\nHSP data not loaded.\n\n"; return; };
my($hsp);
foreach $hsp ( $self->hsps() ) {
$hspCount++;
print $OUT "\n ", '-'x25, "< HSP #$hspCount >", '-'x25, "\n";
$hsp->display( %param );
if( $hspCount < $self->num_hsps ) {
print "\n\n\t--------> <RET> FOR NEXT HSP, q TO QUIT <--------\n";
chop( $reply = <STDIN>);
$reply =~ /^q/i and return;
}
} } | sub _display_stats
{
my( $self, $OUT) = @_;
my $layout = $self->parent->_layout();
if($layout == 1) {
printf( $OUT "%-3d %-20s %-11s %-5d %-5d %-9.1e %-9.1e %-4d %-3d %-5d %-5d %-5s %-6.2f %-6.2f %-4d(%.2f) %-4d(%.2f)\n",
$self->rank(), $self->name(),
($self->database() || 'UNKNOWN DB') ,
$self->score(),$self->bits(),$self->p(),$self->expect(),
$self->gaps(), $self->n(),
$self->length(), $self->length_aln('query'),
$self->ambiguous_aln(),
$self->frac_aligned_query, $self->frac_aligned_hit,
$self->matches('iden'), $self->frac_identical('query'),
$self->matches('cons'), $self->frac_conserved('query'));
} else {
printf( $OUT "%-3d %-20s %-11s %-5d %-5d %-9.1e %-4d %-3d %-5d %-5d %-5s %-6.2f %-6.2f %-4d(%.2f) %-4d(%.2f)\n",
$self->rank(), $self->name(),
($self->database() || 'UNKNOWN DB'),
$self->score(),$self->bits(),$self->expect(),
$self->gaps(), $self->num_hsps,
$self->length(), $self->length_aln('query'),
$self->ambiguous_aln(),
$self->frac_aligned_query, $self->frac_aligned_hit,
$self->matches('iden'), $self->frac_identical('query'),
$self->matches('cons'), $self->frac_conserved('query') );
}
} | sub _initialize
{
my( $self, %param ) = @_;
my $make = $self->SUPER::_initialize( %param );
$_prog = $param{-PROGRAM} || '';
$_signif_fmt = $param{-SIGNIF_FMT};
$self->{'_rank'} = $param{-RANK} || '';
$self->_set_id( $param{-DATA}->[0]);
$self->_set_hsps( @{$param{-DATA}} );
$self->{'_overlap'} = $param{-OVERLAP} || 0;
}
} | sub _set_desc
{
my( $self, @desc ) = @_;
my( $desc);
$desc = join(" ", @desc);
if($desc) {
$desc =~ s/^\s*\S+\s+//;
$desc =~ s/^[\s!]+//;
$desc =~ s/ \d+$//;
$desc =~ s/\.+$//;
$self->{'_desc'} = $desc || '-';
} else {
$self->{'_desc'} = '-';
}
} | sub _set_hsps
{
my( $self, @data ) = @_;
my $start = 0;
my $hspCount = 0;
require Bio::Tools::Blast::HSP;
my( @hspData, @hspList, @errs, @bad_names );
my($line, $set_desc, @desc);
$set_desc = 0;
hit_loop:
foreach $line( @data ) {
if( $line =~ /^\s*Length = ([\d,]+)/ ) {
$self->_set_desc(@desc);
$set_desc = 1;
($self->{'_length'} = $1) =~ s/,//g;
next hit_loop;
} elsif( !$set_desc) {
$line =~ s/^\s+|\s+$//g;
push @desc, $line;
next hit_loop;
} elsif( $line =~ /^\s*Score/ ) {
if( not scalar @hspData ) {
$start = 1;
push @hspData, $line;
next hit_loop;
} elsif( scalar @hspData) {
$hspCount++;
$DEBUG and do{ print STDERR +( $hspCount % 10 ? "+" : "+\n" ); };
my $hspObj = eval { new Bio::Tools::Blast::HSP(-DATA =>\@hspData,
-PARENT =>$self,
-NAME =>$hspCount,
-PROGRAM =>$_prog,
);
};
if($@) {
push @errs, $@;
push @bad_names, "#$hspCount";
$hspObj->destroy if ref $hspObj;
undef $hspObj;
} else {
push @hspList, $hspObj;
if (!defined($self->{'_expect'}) || $hspObj->expect() < $self->{'_expect'}) {
$self->{'_expect'} = $hspObj->expect();
}
if (!defined($self->{'_p'}) || $hspObj->p() < $self->{'_p'}) {
$self->{'_p'} = $hspObj->p();
}
}
@hspData = ();
push @hspData, $line;
next;
} else {
push @hspData, $line;
}
} elsif( $start ) {
if( $line =~ /^(end|>|Parameters|CPU|Database:)/ ) {
$hspCount++;
$DEBUG and do{ print STDERR +( $hspCount % 10 ? "+" : "+\n" ); };
my $hspObj = eval { new Bio::Tools::Blast::HSP(-DATA =>\@hspData,
-PARENT =>$self,
-NAME =>$hspCount,
-PROGRAM =>$_prog,
);
};
if($@) {
push @errs, $@;
push @bad_names, "#$hspCount";
$hspObj->destroy if ref $hspObj;
undef $hspObj;
} else {
push @hspList, $hspObj;
if (!defined($self->{'_expect'}) || $hspObj->expect() < $self->{'_expect'}) {
$self->{'_expect'} = $hspObj->expect();
}
if (!defined($self->{'_p'}) || $hspObj->p() < $self->{'_p'}) {
$self->{'_p'} = $hspObj->p();
}
}
} else {
push @hspData, $line;
}
}
}
$self->{'_length'} or $self->throw( "Can't determine hit sequence length.");
if($_prog =~ /TBLAST[NX]/) {
$self->{'_logical_length'} = $self->{'_length'} / 3;
}
if(@errs) {
my ($str);
if(@errs > 2) {
$str = "SHOWING FIRST EXCEPTION ONLY:\n$errs[0]";
$self->clear_err();
} else {
$str = join("\n",@errs);
}
if( not scalar @hspList) {
$self->throw("Failed to create any HSP objects for $hspCount potential HSP(s).",
"\n\nTRAPPED EXCEPTION(S):\n$str\nEND TRAPPED EXCEPTION(S)\n"
);
} else {
$self->warn(sprintf("Could not create HSP objects for %d HSP(s): %s", scalar(@errs), join(', ',@bad_names)),
"\n\nTRAPPED EXCEPTION(S):\n$str\nEND TRAPPED EXCEPTION(S)\n"
);
}
} else {
$self->{'_hsps'} =\@ hspList;
} } | sub _set_id
{
my( $self, $desc ) = @_;
my ($seqID1, $seqID2, $dbID, @seqDat);
local $_ = $desc;
my @linedat = split();
my $data = $linedat[0];
if($data =~ /^(\S+)\s*/) {
$self->name($1);
}
else {
$self->warn("Can't locate sequence identifier in summary line.", "Line = $data");
$data = 'Unknown sequence ID' if not $data;
$self->name($data);
}
$self->{'_db'} = '-';
} | sub _tile_hsps
{
my $self = shift;
$self->{'_tile_hsps'} = 1;
$self->{'_gaps_query'} = 0;
$self->{'_gaps_sbjct'} = 0;
if((defined($self->{'_n'}) and $self->{'_n'} == 1) or $self->num_hsps == 1) {
my $hsp = $self->hsp;
$self->{'_length_aln_query'} = $hsp->length('query');
$self->{'_length_aln_sbjct'} = $hsp->length('sbjct');
$self->{'_length_aln_total'} = $hsp->length('total');
($self->{'_totalIdentical'},$self->{'_totalConserved'}) = $hsp->matches();
$self->{'_gaps_query'} = $hsp->gaps('query');
$self->{'_gaps_sbjct'} = $hsp->gaps('sbjct');
return;
} else {
$self->{'_length_aln_query'} = 0;
$self->{'_length_aln_sbjct'} = 0;
$self->{'_length_aln_total'} = 0;
$self->{'_totalIdentical'} = 0;
$self->{'_totalConserved'} = 0;
}
my($hsp, $qstart, $qstop, $sstart, $sstop);
my(@qcontigs, @scontigs);
my $qoverlap = 0;
my $soverlap = 0;
my $max_overlap = $self->{'_overlap'};
foreach $hsp ($self->hsps()) {
($qstart, $qstop) = $hsp->range('query');
($sstart, $sstop) = $hsp->range('sbjct');
my ($qgaps, $sgaps) = $hsp->gaps();
$self->{'_gaps_query'} += $qgaps;
$self->{'_gaps_sbjct'} += $sgaps;
$self->{'_length_aln_total'} += $hsp->length;
$qoverlap = &_adjust_contigs('query', $hsp, $qstart, $qstop,\@ qcontigs, $max_overlap);
$soverlap = &_adjust_contigs('sbjct', $hsp, $sstart, $sstop,\@ scontigs, $max_overlap);
if(not defined $self->{'_queryStart'}) {
$self->{'_queryStart'} = $qstart;
$self->{'_queryStop'} = $qstop;
$self->{'_sbjctStart'} = $sstart;
$self->{'_sbjctStop'} = $sstop;
} else {
$self->{'_queryStart'} = ($qstart < $self->{'_queryStart'} ? $qstart : $self->{'_queryStart'});
$self->{'_queryStop'} = ($qstop > $self->{'_queryStop'} ? $qstop : $self->{'_queryStop'});
$self->{'_sbjctStart'} = ($sstart < $self->{'_sbjctStart'} ? $sstart : $self->{'_sbjctStart'});
$self->{'_sbjctStop'} = ($sstop > $self->{'_sbjctStop'} ? $sstop : $self->{'_sbjctStop'});
}
}
foreach(@qcontigs) {
$self->{'_length_aln_query'} += $_->{'stop'} - $_->{'start'} + 1;
$self->{'_totalIdentical'} += $_->{'iden'};
$self->{'_totalConserved'} += $_->{'cons'};
}
foreach(@scontigs) {
$self->{'_length_aln_sbjct'} += $_->{'stop'} - $_->{'start'} + 1;
}
if($qoverlap) {
if($soverlap) { $self->ambiguous_aln('qs');
}
else { $self->ambiguous_aln('q');
}
} elsif($soverlap) {
$self->ambiguous_aln('s');
}
my $prog = $self->parent->program;
if($prog eq 'TBLASTN') {
$self->{'_length_aln_sbjct'} /= 3;
} elsif($prog eq 'BLASTX' ) {
$self->{'_length_aln_query'} /= 3;
} elsif($prog eq 'TBLASTX') {
$self->{'_length_aln_query'} /= 3;
$self->{'_length_aln_sbjct'} /= 3;
}
} | sub ambiguous_aln
{
my $self = shift;
if(@_) { $self->{'_ambiguous_aln'} = shift; }
$self->{'_ambiguous_aln'} || '-'; } | sub bits
{
my $self = shift;
my ($bits);
if(not defined($self->{'_bits'})) {
$bits = $self->hsp->bits;
} else {
$bits = $self->{'_bits'};
}
return $bits; } | sub database
{ my $self = shift;
$self->warn("Bio::Tools::Sbjct::database() is deprecated.\nNo useful information is provided by this method.\nUse Bio::Tools::Blast::database().\n");
return $self->{'_db'};
}
} | sub desc
{
my( $self, $len ) = @_;
$len = (defined $len) ? $len : (CORE::length $self->{'_desc'});
substr( $self->{'_desc'}, 0 ,$len ); } | sub destroy
{
my $self=shift;
if($self->{'_hsps'}) {
foreach($self->hsps) {
$_->destroy;
undef $_;
}
undef $self->{'_hsps'};
}
$DEBUG==2 && print STDERR "DESTROYING $self ${\$self->name}";
$self->SUPER::destroy;
}
} | sub display
{
my( $self, %param) = @_;
$param{-HEADER} = 0;
$self->SUPER::display(%param);
$self->show =~ /hsp/i and $self->_display_hsps( %param);
} | sub end
{
my ($self, $seqType) = @_;
$seqType ||= (wantarray ? 'list' : 'query');
if($self->num_hsps == 1) {
return $self->hsp->end($seqType);
} else {
$self->_tile_hsps() if not $self->{'_tile_hsps'};
if($seqType =~ /list|array/i) {
return ($self->{'_queryStop'}, $self->{'_sbjctStop'});
} else {
$seqType = "_\L$seqType\E";
return $self->{$seqType.'Stop'};
}
}
} | sub expect
{
my ($self, $fmt) = @_;
my $val = $self->{'_expect'};
$fmt ||= $_signif_fmt;
defined($val) or $self->throw("Can't get Expect value: HSPs may not have been set.");
return $val if not $fmt or $fmt =~ /^raw/i;
return &get_exponent($val) if $fmt =~ /^exp/i;
return (split (/eE/, $val)) if $fmt =~ /^parts/i;
return $val; } | sub frac_aligned_hit
{
my $self = shift;
$self->_tile_hsps() if not $self->{'_tile_hsps'};
sprintf( "%.2f", $self->{'_length_aln_sbjct'}/$self->logical_length('sbjct'));
}
} | sub frac_aligned_query
{
my $self = shift;
$self->_tile_hsps() if not $self->{'_tile_hsps'};
sprintf( "%.2f", $self->{'_length_aln_query'}/$self->logical_length('query'));
} | sub frac_conserved
{
my ($self, $seqType) = @_;
$seqType ||= 'total';
$seqType = lc($seqType);
$self->_tile_hsps() if not $self->{'_tile_hsps'};
sprintf( "%.2f", $self->{'_totalConserved'}/$self->{'_length_aln_'.$seqType});
} | sub frac_identical
{
my ($self, $seqType) = @_;
$seqType ||= 'total';
$seqType = lc($seqType);
$self->_tile_hsps() if not $self->{'_tile_hsps'};
sprintf( "%.2f", $self->{'_totalIdentical'}/$self->{'_length_aln_'.$seqType});
} | sub frame
{
my $self = shift;
my ($frame);
if(not defined($self->{'_frame'})) {
$frame = $self->hsp->frame;
} else {
$frame = $self->{'_frame'};
}
return $frame; } | sub gaps
{
my( $self, $seqType ) = @_;
$seqType ||= (wantarray ? 'list' : 'total');
$self->_tile_hsps() if not $self->{'_tile_hsps'};
$seqType = lc($seqType);
if($seqType =~ /list|array/i) {
return ($self->{'_gaps_query'}, $self->{'_gaps_sbjct'});
}
if($seqType eq 'total') {
return ($self->{'_gaps_query'} + $self->{'_gaps_sbjct'}) || 0;
} else {
return $self->{'_gaps_'.$seqType} || 0;
}
} | sub get_exponent
{
my $data = shift;
my($num, $exp) = split /[eE]/, $data;
if( defined $exp) {
$num = 1 if not $num;
$num >= 5 and $exp++;
$num <= -5 and $exp--;
} elsif( $num == 0) {
$exp = -999;
} elsif( not $num =~ /\./) {
$exp = CORE::length($num) -1;
} else {
$exp = 0;
$num .= '0' if $num =~ /\.$/;
my ($c);
my $rev = 0;
if($num !~ /^0/) {
$num = reverse($num);
$rev = 1;
}
do { $c = chop($num);
$c == 0 && $exp++;
} while( $c ne '.');
$exp = -$exp if $num == 0 and not $rev;
$exp -= 1 if $rev;
}
return $exp;
}
1;
__END__
} | sub homol_data
{
my ($self,%param) = @_;
my $data = $self->name();
foreach ($self->hsps) {
$data .="\t".$_->homol_data(%param);
}
$data .= "\t".$self->ambiguous_aln();
$data;
} | sub hsp
{
my( $self, $option ) = @_;
$option ||= 'best';
if (not ref $self->{'_hsps'}) {
$self->throw("Can't get HSPs: data not collected.");
}
my @hsps = @{$self->{'_hsps'}};
return $hsps[0] if $option =~ /best|first|1/i;
return $hsps[$#hsps] if $option =~ /worst|last/i;
$self->throw("Can't get HSP for: $option",
"Valid arguments: 'best', 'worst'");
} | sub hsps
{
my $self = shift;
if (not ref $self->{'_hsps'}) {
$self->throw("Can't get HSPs: data not collected.");
}
return wantarray
? @{$self->{'_hsps'}}
: scalar(@{$self->{'_hsps'}});
} | sub is_signif
{
my $self = shift;
return ($self->{'_significance'} <= $self->parent->signif and
$self->length > $self->parent->signif_len);
}
} | sub length
{
my $self = shift;
$self->{'_length'};
} | sub length_aln
{
my( $self, $type ) = @_;
$type ||= 'query';
$self->_tile_hsps() if not $self->{'_tile_hsps'};
my $data = $self->{'_length_aln_'.$type};
if(!$data) {
$self->throw("Can't get length aln for sequence type\" $type\"",
"Valid types are 'query', 'sbjct'");
}
$data;
} | sub logical_length
{
my $self = shift;
my $seqType = shift || 'query';
$seqType eq 'sbjct' and return
$self->{'_logical_length'} || $self->{'_length'};
my $qlen = $self->parent->length;
my $prog = $self->parent->program;
if($prog =~ /T?BLASTX/ ) {
$qlen /= 3;
}
return $qlen;
} | sub matches
{
my( $self, $arg) = @_;
my(@data,$data);
if(!$arg) {
@data = ($self->{'_totalIdentical'}, $self->{'_totalConserved'});
return @data if @data;
} else {
if($arg =~ /^id/i) {
$data = $self->{'_totalIdentical'};
} else {
$data = $self->{'_totalConserved'};
}
return $data if $data;
}
$self->throw("Can't get identical or conserved data: no data.");
} | sub n
{
my $self = shift;
my ($n);
if(not defined($self->{'_n'})) {
$n = $self->hsp->n;
} else {
$n = $self->{'_n'};
}
$n ||= $self->num_hsps;
return $n; } | sub num_hsps
{
my $self = shift;
if (not defined $self->{'_hsps'}) {
$self->throw("Can't get HSPs: data not collected.");
}
return scalar(@{$self->{'_hsps'}});
} | sub num_unaligned_hit
{
my $self = shift;
$self->_tile_hsps() if not $self->{'_tile_hsps'};
my $num = $self->logical_length('sbjct') - $self->{'_length_aln_sbjct'};
($num < 0 ? 0 : $num );
} | sub num_unaligned_query
{
my $self = shift;
$self->_tile_hsps() if not $self->{'_tile_hsps'};
my $num = $self->logical_length('query') - $self->{'_length_aln_query'};
($num < 0 ? 0 : $num );
} | sub num_unaligned_sbjct
{ my $self=shift; $self->num_unaligned_hit(@_); } | sub overlap
{
my $self = shift;
if(@_) { $self->{'_overlap'} = shift; }
defined $self->{'_overlap'} ? $self->{'_overlap'} : 0; } | sub p
{
my ($self, $fmt) = @_;
my ($val);
$fmt ||= $_signif_fmt;
if(not defined($val = $self->{'_p'})) {
my $note = '';
if($self->parent->_layout() == 2) {
$note = "Blast2 does not report P-values. Use expect() instead.";
}
$self->throw("Can't get P-value: undefined.", $note);
}
return $val if not $fmt or $fmt =~ /^raw/i;
return &get_exponent($val) if $fmt =~ /^exp/i;
return (split (/eE/, $val)) if $fmt =~ /^parts/i;
return $val; } | sub range
{
my ($self, $seqType) = @_;
$seqType ||= 'query';
return ($self->start($seqType), $self->end($seqType));
} | sub rank
{
my $self = shift;
if(@_) {$self->{'_rank'} = shift; }
$self->{'_rank'};
} | sub score
{
my $self = shift;
my ($score);
if(not defined($self->{'_score'})) {
$score = $self->hsp->score;
} else {
$score = $self->{'_score'};
}
return $score; } | sub seq_inds
{
my ($self, $seq, $class, $collapse) = @_;
$seq ||= 'query';
$class ||= 'identical';
$collapse ||= 0;
my (@inds, $hsp);
foreach $hsp ($self->hsps) {
push @inds, $hsp->seq_inds($seq, $class);
}
if(@inds) {
my %tmp = map { $_, 1 } @inds;
@inds = sort {$a <=> $b} keys %tmp;
}
require Bio::Tools::Blast::HSP;
$collapse ? &Bio::Tools::Blast::HSP::collapse_nums(@inds) : @inds;
}
} | sub signif
{
my ($self, $fmt) = @_;
my $val = defined($self->{'_p'}) ? $self->{'_p'} : $self->{'_expect'};
$fmt ||= $_signif_fmt;
defined($val) or $self->throw("Can't get P- or Expect value: HSPs may not have been set.");
return $val if not $fmt or $fmt =~ /^raw/i;
return &get_exponent($val) if $fmt =~ /^exp/i;
return (split (/eE/, $val)) if $fmt =~ /^parts/i;
return $val;
} | sub start
{
my ($self, $seqType) = @_;
$seqType ||= (wantarray ? 'list' : 'query');
if($self->num_hsps == 1) {
return $self->hsp->start($seqType);
} else {
$self->_tile_hsps() if not $self->{'_tile_hsps'};
if($seqType =~ /list|array/i) {
return ($self->{'_queryStart'}, $self->{'_sbjctStart'});
} else {
$seqType = "_\L$seqType\E";
return $self->{$seqType.'Start'};
}
}
} | General documentation
Bio::Tools::Blast::Sbjct.pm is a concrete class that inherits from
Bio::Root::Objectand relies on two other modules:
Bio::Tools::Blast::HSP
Encapsulates a single high-scoring segment pair within a hit.
Bio::Tools::Blast
Provides a container for Sbjct.pm objects.
Bio::Tools::Blast::Sbjct.pm does not currently inherit from
Bio::Root::Vector.pm since Bio::Root::Vector.pm may be re-designed to
make it usable via delegation. Thus, a Blast.pm object would manage a
vector of Sbjct.pm objects. Stay tuned.
One consequence of the HSP tiling is that methods that rely on HSP
tiling such as
frac_identical(),
frac_conserved(),
gaps()etc. may report misleading numbers when -OVERLAP is set to a large
number. For example, say we have two HSPs and the query sequence tile
as follows:
1 8 22 30 40 60
Full seq: ------------------------------------------------------------
* ** * **
HSP1: --------------- (6 identical matches)
** ** **
HSP2: ------------- (6 identical matches)
If -OVERLAP is set to some number over 4, HSP1 and HSP2 will not be
tiled into a single contig and their numbers of identical matches will
be added, giving a total of 12, not 10 if they had be combined into
one contig. This can lead to number greater than 1.0 for methods
frac_identical() and
frac_conserved(). This is less of an issue
with gapped Blast since it tends to combine HSPs that would be listed
separately without gapping. (Fractions >1.0 can be viewed as a
signal for an interesting alignment that warrants further inspection,
thus turning this bug into a feature).
Using large values for -OVERLAP can lead to incorrect numbers
reported by methods that rely on HSP tiling but can be useful if you
care more about detecting ambiguous alignments. Setting -OVERLAP
to zero will lead to the most accurate numbers for the
tiling-dependent methods but will be useless for detecting overlapping
HSPs since all HSPs will appear to overlap.
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to one
of the Bioperl mailing lists. Your participation is much appreciated.
bioperl-l@bioperl.org - General discussion
http://bio.perl.org/MailList.html - About the mailing lists
Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution. Bug reports can be submitted via email
or the web:
bioperl-bugs@bio.perl.org
http://bugzilla.bioperl.org/
Copyright (c) 1996-2000 Steve Chervitz. All Rights Reserved.
This module is free software; you can redistribute it and/or
modify it under the same terms as Perl itself.
Methods beginning with a leading underscore are considered private
and are intended for internal use by this module. They are
not considered part of the public interface and are described here
for documentation purposes only.
Information about the various data members of this module is provided for those
wishing to modify or understand the code. Two things to bear in mind:
1 Do NOT rely on these in any code outside of this module.
All data members are prefixed with an underscore to signify that they are private.
Always use accessor methods. If the accessor doesn't exist or is inadequate,
create or modify an accessor (and let me know, too!). (An exception to this might
be for HSP.pm which is more tightly coupled to Sbjct.pm and
may access Sbjct data members directly for efficiency purposes, but probably
should not).
2 This documentation may be incomplete and out of date.
It is easy for these data member descriptions to become obsolete as
this module is still evolving. Always double check this info and search
for members not described here.
An instance of Bio::Tools::Blast::Sbjct.pm is a blessed reference to a hash containing
all or some of the following fields:
FIELD VALUE
--------------------------------------------------------------
_hsps : Array ref for a list of Bio::Tools::Blast::HSP.pm objects.
:
_db : Database identifier from the summary line.
:
_desc : Description data for the hit from the summary line.
:
_length : Total length of the hit sequence.
:
_score : BLAST score.
:
_bits : BLAST score (in bits). Matrix-independent.
:
_p : BLAST P value. Obtained from summary section. (Blast1/WU-Blast only)
:
_expect : BLAST Expect value. Obtained from summary section.
:
_n : BLAST N value (number of HSPs) (Blast1/WU-Blast2 only)
:
_frame : Reading frame for TBLASTN and TBLASTX analyses.
:
_totalIdentical: Total number of identical aligned monomers.
:
_totalConserved: Total number of conserved aligned monomers (a.k.a. "positives").
:
_overlap : Maximum number of overlapping residues between adjacent HSPs
: before considering the alignment to be ambiguous.
:
_ambiguous_aln : Boolean. True if the alignment of all HSPs is ambiguous.
:
_length_aln_query : Length of the aligned region of the query sequence.
:
_length_aln_sbjct : Length of the aligned region of the sbjct sequence.
INHERITED DATA MEMBERS
----------------------
_name : From Bio::Root::Object.pm. String representing the name of the
: sbjct sequence obtained from the BLAST report.
:
_parent : From Bio::Root::Object.pm. This member contains a reference to the
: Bio::Tools::Blast.pm object to which this hit belongs.