BioPerl
view release on metacpan or search on metacpan
Bio/Assembly/Tools/ContigSpectrum.pm view on Meta::CPAN
# Score a contig spectrum (the more abundant the contigs and the larger their
# size, the larger the score)
my $csp_score = $csp->score( $csp->nof_seq );
=head1 DESCRIPTION
The Bio::Assembly::Tools::ContigSpectrum Perl module enables to
manually create contig spectra, import them from assemblies,
manipulate them, transform between different types of contig spectra
and output them.
Bio::Assembly::Tools::ContigSpectrum is a module to create, manipulate
and output contig spectra, assembly-derived data used in metagenomics
(community genomics) for diversity estimation.
=head2 Background
A contig spectrum is the count of the number of contigs of different
size in an assembly. For example, the contig spectrum [100 5 1 0 0
...] means that there were 100 singlets (1-contigs), 5 contigs of 2
sequences (2-contigs), 1 contig of 3 sequences (3-contig) and no
larger contigs.
An assembly can be produced from a mixture of sequences from different
metagenomes. The contig obtained from this assembly is a mixed contig
spectrum. The contribution of each metagenome in this mixed contig
spectrum can be obtained by determining a dissolved contig spectrum.
Finally, based on a mixed contig spectrum, a cross contig spectrum can
be determined. In a cross contig spectrum, only contigs containing
sequences from different metagenomes are kept; "pure" contigs are
excluded. Additionally, the total number of singletons (1-contigs)
from each region that assembles with any fragments from other regions
is the number of 1-contigs in the cross contig spectrum.
=head2 Implementation
The simplest representation of a contig spectrum is as a hash
representation where the key is the contig size (number of sequences
making up the contig) and the value the number of contigs of this
size.
In fact, it is useful to have more information associated with the
contig spectrum, hence the Bio::Assembly::Tools::ContigSpectrum module
implements an object containing a contig spectrum hash and additional
information. The get/set methods to access them are:
id contig spectrum ID
nof_rep number of repetitions (assemblies) used
max_size size of (number of sequences in) the largest contig
spectrum hash representation of a contig spectrum
nof_seq number of sequences
avg_seq_len average sequence length
eff_asm_params reports effective assembly parameters
nof_overlaps number of overlaps (needs eff_asm_params)
min_overlap minimum overlap length in a contig (needs eff_asm_params)
min_identity minimum sequence identity percentage (needs eff_asm_params)
avg_overlap average overlap length (needs eff_asm_params)
avg_identity average overlap identity percentage (needs eff_asm_params)
Operations on the contig spectra:
to_string create a string representation of the spectrum
spectrum import a hash contig spectrum
assembly determine a contig spectrum from an assembly, contig or singlet
dissolve calculate a dissolved contig spectrum (depends on assembly)
cross produce a cross contig spectrum (depends on assembly)
add add a contig spectrum to an existing one
average make an average of several contig spectra
score score a contig spectrum: the higher the number of contigs
and the larger their size, the higher the score.
When using operations that rely on knowing "where" (from what
metagenomes) a sequence came from (i.e. when creating a dissolved or
cross contig spectrum), make sure that the sequences used for the
assembly have a name header, e.g. E<gt>metagenome1|seq1,
E<gt>metagenome2|seq1, ...
Note: The following operations require the C<Graph::Undirected> module:
eff_asm_params, cross, dissolve
=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 lists 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
the bugs and their resolution. Bug reports can be submitted via email
or the web:
bioperl-bugs@bio.perl.org
https://github.com/bioperl/bioperl-live/issues
=head1 AUTHOR - Florent E Angly
Email florent_dot_angly_at_gmail_dot_com
=head1 APPENDIX
Bio/Assembly/Tools/ContigSpectrum.pm view on Meta::CPAN
sub nof_overlaps {
my ($self, $nof_overlaps) = @_;
if (defined $nof_overlaps) {
$self->throw("The number of overlaps must be strictly positive. Got ".
"'$nof_overlaps'") if $nof_overlaps < 1;
$self->{'_nof_overlaps'} = $nof_overlaps;
}
$nof_overlaps = $self->{'_nof_overlaps'};
return $nof_overlaps;
}
=head2 min_overlap
Title : min_overlap
Usage : $csp->min_overlap
Function: get/set the assembly minimum overlap length
Returns : integer
Args : integer [optional]
=cut
sub min_overlap {
my ($self, $min_overlap) = @_;
if (defined $min_overlap) {
$self->throw("The minimum of overlap length must be strictly positive. Got".
" '$min_overlap'") if $min_overlap < 1;
$self->{'_min_overlap'} = $min_overlap;
}
$min_overlap = $self->{'_min_overlap'};
return $min_overlap;
}
=head2 avg_overlap
Title : avg_overlap
Usage : $csp->avg_overlap
Function: get/set the assembly average overlap length
Returns : decimal
Args : decimal [optional]
=cut
sub avg_overlap {
my ($self, $avg_overlap) = @_;
if (defined $avg_overlap) {
$self->throw("The average overlap length must be strictly positive. Got ".
"'$avg_overlap'") if $avg_overlap < 1;
$self->{'_avg_overlap'} = $avg_overlap;
}
$avg_overlap = $self->{'_avg_overlap'};
return $avg_overlap;
}
=head2 min_identity
Title : min_identity
Usage : $csp->min_identity
Function: get/set the assembly minimum overlap identity percent
Returns : 0 < decimal < 100
Args : 0 < decimal < 100 [optional]
=cut
sub min_identity {
my ($self, $min_identity) = @_;
if (defined $min_identity) {
$self->throw("The minimum overlap percent identity must be strictly ".
"positive. Got '$min_identity'") if $min_identity < 1;
$self->{'_min_identity'} = $min_identity;
}
$min_identity = $self->{'_min_identity'};
return $min_identity;
}
=head2 avg_identity
Title : avg_identity
Usage : $csp->avg_identity
Function: get/set the assembly average overlap identity percent
Returns : 0 < decimal < 100
Args : 0 < decimal < 100 [optional]
=cut
sub avg_identity {
my ($self, $avg_identity) = @_;
if (defined $avg_identity) {
$self->throw("The average overlap percent identity must be strictly ".
"positive. Got '$avg_identity'") if $avg_identity < 1;
$self->{'_avg_identity'} = $avg_identity;
}
$avg_identity = $self->{'_avg_identity'};
return $avg_identity;
}
=head2 avg_seq_len
Title : avg_seq_len
Usage : $csp->avg_seq_len
Function: get/set the assembly average sequence length
Returns : avg_seq_len
Args : real [optional]
=cut
sub avg_seq_len {
my ($self, $avg_seq_len) = @_;
if (defined $avg_seq_len) {
$self->throw("The average sequence length must be strictly positive. Got ".
"'$avg_seq_len'") if $avg_seq_len < 1;
$self->{'_avg_seq_len'} = $avg_seq_len;
}
$avg_seq_len = $self->{'_avg_seq_len'};
return $avg_seq_len;
}
=head2 eff_asm_params
Title : eff_asm_params
Usage : $csp->eff_asm_params(1)
Function: Get/set the effective assembly parameters option. It defines if the
effective assembly parameters should be determined when a contig
spectrum based or derived from an assembly is calculated. The
effective assembly parameters include avg_seq_length, nof_overlaps,
min_overlap, avg_overlap, min_identity and avg_identity.
1 = get them, 0 = don't.
Returns : integer
Args : integer [optional]
=cut
sub eff_asm_params {
my ($self, $eff_asm_params) = @_;
if (defined $eff_asm_params) {
$self->throw("eff_asm_params can only take values 0 or 1. Input value was ".
"'$eff_asm_params'") unless $eff_asm_params == 0 || $eff_asm_params == 1;
$self->{'_eff_asm_params'} = $eff_asm_params;
}
$eff_asm_params = $self->{'_eff_asm_params'};
return $eff_asm_params;
}
=head2 spectrum
Title : spectrum
Bio/Assembly/Tools/ContigSpectrum.pm view on Meta::CPAN
# Average number of sequences
$avg->{'_nof_seq'} /= $avg->{'_nof_rep'};
# Average number of overlaps
$avg->{'_nof_overlaps'} /= $avg->{'_nof_rep'};
return $avg;
}
=head2 score
Title : score
Usage : my $score = $csp->score();
Function: Score a contig spectrum (or cross-contig spectrum) such that the
higher the number of contigs (or cross-contigs) and the larger their
size, the higher the score.
Let n : total number of sequences
c_q : number of contigs of size q
q : number of sequence in a contig
We define: score = n/(n-1) * (X - 1/n)
where X = sum ( c_q * q^2 ) / n**2
The score ranges from 0 (singlets only) to 1 (a single large contig)
It is possible to specify a value for the number of sequences to
assume in the contig spectrum.
Returns : contig score, or undef if there were no sequences in the contig spectrum
Args : number of total sequences to assume [optional]
=cut
sub score {
my ($self, $nof_seqs) = @_;
# Sanity check
my $n = $self->nof_seq;
return undef if ($n <= 0);
# Calculate X
my $score = 0;
my $q_max = $self->max_size;
my $spec = $self->spectrum;
for my $q ( 1 .. $q_max ) {
my $c_q = $spec->{$q};
if ( $q == 1 && $nof_seqs ) {
$c_q += $nof_seqs - $n;
$n = $nof_seqs;
}
next if not $c_q;
$score += $c_q * $q ** 2;
}
$score /= $n ** 2;
# Rescale X to obtain the score
$score = $n/($n-1) * ($score - 1/$n);
return $score;
}
=head2 _naive_assembler
Title : _naive_assembler
Usage :
Function: Reassemble the specified sequences only based on their position in
the contig. This naive assembly only verifies that the minimum
overlap length and percentage identity are respected. No actual
alignment is done
Returns : arrayref of contigs and singlets
Args : Bio::Assembly::Contig
array reference of sequence IDs to use [optional]
minimum overlap length (integer) [optional]
minimum percentage identity (integer) [optional]
=cut
sub _naive_assembler {
my ($self, $contig, $seqlist, $min_overlap, $min_identity) = @_;
# Use all reads if none was specified:
if (not defined $seqlist) {
for my $seq ($contig->each_seq) {
push @$seqlist, $seq->id;
}
}
# Sanity checks
if ( ! ref $seqlist || ! ref($seqlist) eq 'ARRAY') {
$self->throw('Expecting an array reference. Got ['.ref($seqlist)."] \n");
}
my $max = scalar @$seqlist;
$self->throw("Expecting at least 2 sequences as input for _naive_assembler")
if ($max < 2);
# Build contig graph
my %seq_hash = map { $_ => undef } (@$seqlist) if (scalar @$seqlist > 0);
my ($g, $overlaps) = $self->_contig_graph($contig, \%seq_hash, $min_overlap, $min_identity);
# Construct sub-contigs
my @contig_objs;
my $num = 1;
if (defined $g) {
for my $connected_reads ($g->connected_components) { # reads that belong in contigs
my $sub_id = $contig->id.'_'.$num;
my $sub_contig = $self->_create_subcontig($contig, $connected_reads, $sub_id);
push @contig_objs, $sub_contig;
$num++;
for my $read_id ( @$connected_reads ) {
delete $seq_hash{$read_id};
}
}
}
# Construct sub-singlets
my @singlet_objs;
for my $read_id ( keys %seq_hash ) {
my $read = $contig->get_seq_by_name($read_id);
my $sub_singlet = Bio::Assembly::Singlet->new(
-id => $contig->id.'_'.$num,
-seqref => $self->_obj_copy($read)
);
$num++;
push @singlet_objs, $sub_singlet;
}
return [@contig_objs, @singlet_objs];
}
=head2 _create_subcontig
Title : _create_subcontig
Usage :
Bio/Assembly/Tools/ContigSpectrum.pm view on Meta::CPAN
my @contig_stats = (0, 0);
# contig_stats = (avg_length, nof_seq)
for my $seqobj ($contigobj->each_seq) {
next if defined $seq_hash && !defined $$seq_hash{$seqobj->id};
my $seq_string;
if ($contigobj->isa('Bio::Assembly::Singlet')) { # a singlet
$seq_string = $contigobj->seqref->seq;
} else { # a contig
$seq_string = $seqobj->seq;
}
# Number of non-gap characters in the sequence
my $seq_len = $seqobj->_ungapped_len;
my @seq_stats = ($seq_len);
@contig_stats = $self->_update_seq_stats(@contig_stats, @seq_stats);
}
return @contig_stats;
}
=head2 _update_seq_stats
Title : _update_seq_stats
Usage :
Function: Update the number of sequences and their average length 1
average identity 1
minimum length 1
minimum identity 1
number of overlaps 1 average sequence length
Returns : average sequence length
number of sequences
Args : average sequence length 1
number of sequences 1
average sequence length 2
number of sequences 2
=cut
sub _update_seq_stats {
my ($self, $p_avg_length, $p_nof_seq, $n_avg_length, $n_nof_seq) = @_;
# Defaults
if (not defined $n_nof_seq) {
$n_nof_seq = 1;
}
# Update overlap statistics
my $avg_length = 0;
my $nof_seq = $p_nof_seq + $n_nof_seq;
if ($nof_seq != 0) {
$avg_length = ($p_avg_length * $p_nof_seq + $n_avg_length * $n_nof_seq) / $nof_seq;
}
return $avg_length, $nof_seq;
}
=head2 _get_assembly_overlap_stats
Title : _get_assembly_overlap_stats
Usage : my ($avglength, $avgidentity, $minlength, $min_identity, $nof_overlaps)
= $csp->_get_assembly_overlap_stats($assemblyobj);
Function: Get statistics about pairwise overlaps in contigs of an assembly
Returns : average overlap length
average identity percent
minimum overlap length
minimum identity percent
number of overlaps
Args : Bio::Assembly::Scaffold, Contig or Singlet object
hash reference with the IDs of the sequences to consider [optional]
=cut
sub _get_assembly_overlap_stats {
my ($self, $assembly_obj, $seq_hash) = @_;
# Sanity check
if ( !defined $assembly_obj ||
( !$assembly_obj->isa('Bio::Assembly::ScaffoldI') &&
!$assembly_obj->isa('Bio::Assembly::Contig') ) ) {
$self->throw("Must provide a Bio::Assembly::ScaffoldI, Contig or Singlet object");
}
$self->throw("Expecting a hash reference. Got [".ref($seq_hash)."]")
if (defined $seq_hash && ! ref($seq_hash) eq 'HASH');
# Look at all the contigs (no singlets!)
my @asm_stats = (0, 0, undef, undef, 0);
# asm_stats = (avg_length, avg_identity, min_length, min_identity, nof_overlaps)
for my $contig_obj ( $self->_get_contig_like($assembly_obj) ) {
@asm_stats = $self->_update_overlap_stats( @asm_stats,
$self->_get_contig_overlap_stats($contig_obj, $seq_hash) );
}
return @asm_stats;
}
=head2 _get_contig_overlap_stats
Title : _get_contig_overlap_stats
Usage : my ($avglength, $avgidentity, $minlength, $min_identity, $nof_overlaps)
= $csp->_get_contig_overlap_stats($contigobj);
Function: Get statistics about pairwise overlaps in a contig or singlet. The
statistics are obtained using graph theory: each read is a node
and the edges between 2 reads are weighted by minus the number of
conserved residues in the alignment between the 2 reads. The
minimum spanning tree of this graph represents the overlaps that
form the contig. Overlaps that do not satisfy the minimum overlap
length and similarity get a malus on their score.
Note: This function requires the optional BioPerl dependency
module called 'Graph'
Returns : average overlap length
average identity percent
minimum overlap length
minimum identity percent
number of overlaps
Args : Bio::Assembly::Contig or Singlet object
hash reference with the IDs of the sequences to consider [optional]
=cut
sub _get_contig_overlap_stats {
my ($self, $contig_obj, $seq_hash) = @_;
# Sanity check
$self->throw("Must provide a Bio::Assembly::Contig object")
if (!defined $contig_obj || !$contig_obj->isa("Bio::Assembly::Contig"));
$self->throw("Expecting a hash reference. Got [".ref($seq_hash)."]")
if (defined $seq_hash && ! ref($seq_hash) eq 'HASH');
my @contig_stats = (0, 0, undef, undef, 0);
# contig_stats = (avg_length, avg_identity, min_length, min_identity, nof_overlaps)
# Build contig graph
### consider providing the minima to _contig_graph here too?
my ($g, $overlaps) = $self->_contig_graph($contig_obj, $seq_hash);
if ( defined $g ) {
# Graph minimum spanning tree (tree that goes through strongest overlaps)
$g = $g->MST_Kruskal();
# Calculate minimum overlap length and identity for this contig
for my $edge ( $g->edges ) {
# Retrieve overlap information
my ($id1, $id2) = @$edge;
if (not exists $$overlaps{$id1}{$id2}) {
($id2, $id1) = @$edge;
}
my ($score, $length, $identity) = @{$$overlaps{$id1}{$id2}};
# Update contig stats
my @overlap_stats = ($length, $identity);
@contig_stats = $self->_update_overlap_stats(@contig_stats, @overlap_stats);
}
}
return @contig_stats;
}
=head2 _update_overlap_stats
Title : _update_overlap_stats
Usage :
Function: update the number of overlaps and their minimum and average length
and identity
Returns :
Args : average length 1
average identity 1
minimum length 1
minimum identity 1
number of overlaps 1
average length 2
average identity 2
minimum length 2
minimum identity 2
number of overlaps 2
=cut
sub _update_overlap_stats {
my ($self,
$p_avg_length, $p_avg_identity, $p_min_length, $p_min_identity, $p_nof_overlaps,
$n_avg_length, $n_avg_identity, $n_min_length, $n_min_identity, $n_nof_overlaps)
= @_;
# Defaults
if (not defined $n_nof_overlaps) { $n_nof_overlaps = 1 };
if ((not defined $n_min_length) && ($n_avg_length != 0)) { $n_min_length = $n_avg_length };
if ((not defined $n_min_identity) && ($n_avg_identity != 0)) { $n_min_identity = $n_avg_identity };
# Update overlap statistics
my ($avg_length, $avg_identity, $min_length, $min_identity, $nof_overlaps)
= (0, 0, undef, undef, 0);
$nof_overlaps = $p_nof_overlaps + $n_nof_overlaps;
if ($nof_overlaps > 0) {
$avg_length = ($p_avg_length * $p_nof_overlaps + $n_avg_length * $n_nof_overlaps) / $nof_overlaps;
$avg_identity = ($p_avg_identity * $p_nof_overlaps + $n_avg_identity * $n_nof_overlaps) / $nof_overlaps;
}
if ( not defined $p_min_length ) {
$min_length = $n_min_length;
} elsif ( not defined $n_min_length ) {
$min_length = $p_min_length;
} else { # both values are defined
if ($n_min_length < $p_min_length) {
$min_length = $n_min_length;
} else {
$min_length = $p_min_length;
}
}
if ( not defined $p_min_identity ) {
$min_identity = $n_min_identity;
} elsif ( not defined $n_min_identity ) {
$min_identity = $p_min_identity;
} else { # both values are defined
if ($n_min_identity < $p_min_identity) {
$min_identity = $n_min_identity;
} else {
$min_identity = $p_min_identity;
}
}
return $avg_length, $avg_identity, $min_length, $min_identity, $nof_overlaps;
}
=head2 _overlap_alignment
Title : _overlap_alignment
Usage :
Function: Produce an alignment of the overlapping section of two sequences of
a contig. Minimum overlap length and percentage identity can be
specified. Return undef if the sequences do not overlap or do not
meet the minimum overlap criteria.
Return : Bio::SimpleAlign object reference
alignment overlap length
alignment overlap identity
Args : Bio::Assembly::Contig object reference
Bio::LocatableSeq contig sequence 1
Bio::LocatableSeq contig sequence 2
minimum overlap length [optional]
minimum overlap identity percentage[optional]
=cut
sub _overlap_alignment {
my ($self, $contig, $qseq, $tseq, $min_overlap, $min_identity) = @_;
# get query and target sequence position
my $qpos = $contig->get_seq_coord($qseq);
my $tpos = $contig->get_seq_coord($tseq);
# check that there is an overlap
my $qend = $qpos->end;
my $tstart = $tpos->start;
return if $qend < $tstart;
my $qstart = $qpos->start;
my $tend = $tpos->end;
return if $qstart > $tend;
# get overlap boundaries and check overlap length
my $left;
if ($qstart >= $tstart) {
$left = $qstart
} else {
$left = $tstart;
}
my $right;
if ($qend > $tend) {
$right = $tend;
} else {
$right = $qend;
}
my $overlap = $right - $left + 1;
return if defined $min_overlap && $overlap < $min_overlap;
# slice query and target sequence to overlap boundaries
my $qleft =
$contig->change_coord('gapped consensus', "aligned ".$qseq->id, $left);
my $qstring = substr($qseq->seq, $qleft - 1, $overlap);
my $tleft =
$contig->change_coord('gapped consensus', "aligned ".$tseq->id, $left);
my $tstring = substr($tseq->seq, $tleft - 1, $overlap);
# remove gaps present in both sequences at the same position
for (my $pos = 0 ; $pos < $overlap ; $pos++) {
my $qnt = substr($qstring, $pos, 1);
if ($qnt eq '-') {
my $tnt = substr($tstring, $pos, 1);
if ($tnt eq '-') {
substr($qstring, $pos, 1, '');
substr($tstring, $pos, 1, '');
$pos--;
$overlap--;
}
}
}
return if defined $min_overlap && $overlap < $min_overlap;
# count the number of gaps remaining in each sequence
my $qgaps = ($qstring =~ tr/-//);
my $tgaps = ($tstring =~ tr/-//);
# make an alignment object with the query and target sequences
my $aln = Bio::SimpleAlign->new;
my $alseq = Bio::LocatableSeq->new(
-id => 1,
-seq => $qstring,
-start => 1,
-end => $overlap - $qgaps,
-alphabet => 'dna',
);
$aln->add_seq($alseq);
$alseq = Bio::LocatableSeq->new(
-id => 2,
-seq => $tstring,
-start => 1,
-end => $overlap - $tgaps,
-alphabet => 'dna',
);
$aln->add_seq($alseq);
# check overlap percentage identity
my $identity = $aln->overall_percentage_identity;
return if defined $min_identity && $identity < $min_identity;
# all checks passed, return alignment
return $aln, $overlap, $identity;
}
=head2 _contig_graph
Title : _contig_graph
Usage :
Function: Creates a graph data structure of the contig.The graph is undirected.
The vertices are the reads of the contig and edges are the overlap
between the reads. The edges are weighted by the opposite of the
overlap, so it is negative and the better the overlap, the lower the
weight.
Return : Graph object or undef
hashref of overlaps (score, length, identity) for each read pair
Args : Bio::Assembly::Contig object reference
hash reference with the IDs of the sequences to consider [optional]
minimum overlap length (integer) [optional]
minimum percentage identity (integer) [optional]
=cut
sub _contig_graph {
my ($self, $contig_obj, $seq_hash, $min_overlap, $min_identity) = @_;
# Sanity checks
if( !ref $contig_obj || ! $contig_obj->isa('Bio::Assembly::Contig') ) {
$self->throw("Unable to process non Bio::Assembly::Contig ".
"object [".ref($contig_obj)."]");
}
if (not eval { require Graph::Undirected }) {
$self->throw("Error: the module 'Graph' is needed by the method ".
"_contig_graph but could not be found\n$@");
}
# Skip contigs of 1 sequence (they have no overlap)
my @seq_objs = $contig_obj->each_seq;
my $nof_seqs = scalar @seq_objs;
return if ($nof_seqs <= 1);
# Calculate alignment between all pairs of reads
my %overlaps;
for my $i (0 .. $nof_seqs-1) {
my $seq_obj = $seq_objs[$i];
my $seq_id = $seq_obj->id;
# Skip this read if not in list of wanted sequences
next if defined $seq_hash && !exists $$seq_hash{$seq_id};
# What is the best sequence to align to?
my ($best_score, $best_length, $best_identity);
for my $j ($i+1 .. $nof_seqs-1) {
# Skip this sequence if not in list of wanted sequences
my $target_obj = $seq_objs[$j];
my $target_id = $target_obj->id;
next if defined $seq_hash && !exists $$seq_hash{$target_id};
# How much overlap with this sequence?
my ($aln_obj, $length, $identity)
= $self->_overlap_alignment($contig_obj, $seq_obj, $target_obj, $min_overlap, $min_identity);
next if ! defined $aln_obj; # there was no sequence overlap or overlap not good enough
# Score the overlap as the number of conserved residues. In practice, it
# seems to work better than giving +1 for match and -3 for errors
# (mismatch or indels)
my $score = $length * $identity / 100;
# Apply a malus (square root) for scores that do not satisfy the minimum
# overlap length similarity. It is necessary for overlaps that get a high
# score without satisfying both the minimum values.
if ( ( $min_overlap && ($length < $min_overlap ) ) ||
( $min_identity && ($identity < $min_identity) ) ) {
$score = sqrt($score);
}
$overlaps{$seq_id}{$target_id} = [$score, $length, $identity];
( run in 2.481 seconds using v1.01-cache-2.11-cpan-39bf76dae61 )