ZeroDivisionError

I m loading ~ 150 aligned sequences and want to generate a Phylotree,
yet at some point, a ZeroDivisionError is raised.

The stacktrace:

Traceback (most recent call last):
  File "/Users/marcel/dev/prototype/wrapper.py", line 56, in tree
    self._tree = DistanceTreeConstructor().nj(self.distance_matrix)
  File "/Users/marcel/dev/prototype/wrapper.py", line 49, in distance_matrix
    self._distance_matrix = DistanceCalculator('blosum62').get_distance(self.msa)
  File "/Users/marcel/dev/venv/prototype/lib/python2.7/site-packages/Bio/Phylo/TreeConstruction.py", line 465, in get_distance
    dm[seq1.id, seq2.id] = self._pairwise(seq1, seq2)
  File "/Users/marcel/dev/venv/prototype/lib/python2.7/site-packages/Bio/Phylo/TreeConstruction.py", line 448, in _pairwise
    return 1 - (score * 1.0 / max_score)
ZeroDivisionError: float division by zero

The sequences which seem to fail are the following:
seq1:

-------------------------------------------------------GCTTGCTGCTTGGCTGACGAGTGGCGGACGGGTGAGTAATGT-CTGGGAAACTG-----CCTGATGGAGGGGGATAAC-TACTGGAAACGGT-AGCTAATA--CCGCATAA----TGT----CGCAAGA--CCAAAGAGGGGGAC--CTTCGGGCCTCTTGCCAT-CGGATGTGCCC-AGATGGGATTAGC-TAGTTGGTGGGGT---AACGGCTC---ACCTAGGCGACG--ATCCCTAGCTGGTCTGAGA-GGATGAC--CAGCCACA--CTGGAACTGAGACAC---GGTCCAGACTCCTACGGGAGGCAGCAG-TGGGGAATAT--------TGCACAA--TGGGCG---CAAGCCTG-ATGCAGC-------CATGCCGCGTGTATGAA-GAAGGCCTTCGGGTTGTAAAGTACTTTCAGCGGGGAGGAAGGG-AGTAAAGTT-AATACCTTT----GCTCATTGACGTTACCCGCAGAAGAAGCACC-GGCTAAC-TCCGTGCCAGCAGCCGCGGTAATACGG-A-GGGTGC-AAGCGTTAATCGG-AATTACT-GGGCGTAAAGCGC-ACGC-AGGC-GGT--TTGTT-AAGTC-AGATGTGAAATC-----CCCG-GGCT----CAA-CCTGGGAACTGCATCTGATAC--TGG-CAA--------GCTTGAGTCT---CGTAGAGGGGGGTAGAA-TTCCAGGTGT-AGCGGTGAAATGCGTAGAGATCTGG-AGGAATA---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

seq2:

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CCGGAGCTAACGCGTT-AAGTCGACCGCC-TGGG-----------GAGTACGGCCGCAAGGTTAAAACTCAAATGAATTGACG-----GGGGCCCGCACAAGCGGT-GGAGCA-TGTGGTTTAATTCGATGCAACGCGAAGAACCTTA-CCTACTC-TTGACATCCAGA-GAACTTAGCAGAG-ATG--CTATGGTGCCTTCGGG-AACTCTGAGACA-GGTGCT-GCAT-GGCTGTCGTCAGCTCGTGTTGTG-AAATGTT-GGGTTAAGTCCCGCAACGA-GCGC-AACCCTTATCCTTTGTTGCCAGCG-GTTCGGTCGGGAACTCAAAGGAGACTGCCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCC-TTACGAGTAGGGCTACACACGTGCTACAATGGC-ATATACAAAGAGAAGC-GACCTCGCGAGAGCAA-GCGGACCTCATAAAGTAT-GTCGTAGTCCGGATCGGAGTCTGCAACTCG-ACTCCGTGAAGTCGGAAT-CGCTAGTAATCGTGGATCA-GAATGCCA-----CGGTGAATACGTTCCCGGGCCTTGTACACACCGCCCGTCACACCATGGGAGTGGGTTGCAAAAGAA------------------------------------------------------------------

What is interesting about that pair of aligned sequences is every position (column in the full alignment) has a gap in either sequence. Nowhere can two bases be compared directly to get a score.

So, also long as gap-vs-gap scores zero, and gap-vs-base scores zero, this pairwise alignment scores zero. I guess we need to special case this to return the maximum possible distance?

@etal

@etal how about we insert this into the _pairwise method?

if len(seq1) == 0:
    return 1  # max possible scaled distance

@etal

Yes,I think that should do it.
On Mar 21, 2015 3:19 AM, "Peter Cock" [email protected] wrote:

@etal https://github.com/etal how about we insert this into the
_pairwise method?

if len(seq1) == 0:
return 1 # max possible scaled distance

—
Reply to this email directly or view it on GitHub
#491 (comment).

@etal

@etal - Can you construct a little test case for this, and verify the fix?

Will this case ever get fulfilled?
In the example, both sequences have a length of 1723.
As a fallback I would propose the following:

def _pairwise(self, seq1, seq2):
        """Calculate pairwise distance from two sequences"""
        score = 0
        max_score = 0
        if ....
            ....
        else ....
            ....
        # everything failed so far
        if max_score == 0:
            return 1

        return 1 - (score * 1.0 / max_score)

or

def _pairwise(self, seq1, seq2):
        """Calculate pairwise distance from two sequences"""
        score = 0
        max_score = 0
        if ....
            ....
        else ....
            ....
        try:
            return 1 - (score * 1.0 / max_score)
        except ZeroDivisionError, e:
            return 1

Sorry, I'm not going to have time to test this myself within a reasonable amount of time.

@etal

@etal @peterjc I can take care of this if it is still open.

@kianho

Yes, please @kianho - another set of eyes on this code would be great.

Avoid a ZeroDivisionError in TreeConstruction.DistanceCalculator (#491)

Add unit test for TreeConstructor.DistanceCalculator edge case (#491)

@lijax

I pushed a fix and a test.

I noticed an interesting behavior that @lijax must have intended (?), but wasn't obvious at first -- if no scoring matrix is specified, the distances are based on character identity, which is like the "identity" scoring scheme but without skipping gap characters. The test includes checking this behavior, too. It would be more clear to have another name for this type of scoring (maybe it's useful for odd encodings of phenotypic data), e.g. character, and for this _distance method to simply fall back to either the character or identity scoring matrix if self.scoring_matrix isn't specified.

@etal

Sorry, I'm a little too late to notice this issue.

Yes, @etal . That was Intended for some general usage. I agree with your suggestion. It'll be more obvious to give it a name. How about making skip_letters changeable to make this class more generalized?

@lijax

@lijax @etal is this scoring question still relevant, or can we close this issue?

TreeConstruction: Add 'skip_letters' parameter to DistanceCalculator (#…

@lijax

I made the skip_letters a parameter as @lijax suggested. I'll leave it to your discretion if/when the default behavior should change to make the "identity" model skip gap and stop characters in calculating the score.

TreeConstruction: Add 'skip_letters' parameter to DistanceCalculator (b…

Navigation Menu

Repository navigation

ZeroDivisionError #491

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

ZeroDivisionError #491

Description

Activity

peterjc commented on Mar 21, 2015

peterjc commented on Mar 21, 2015

etal commented on Mar 25, 2015

peterjc commented on Mar 25, 2015

quiqua commented on Apr 7, 2015

etal commented on Apr 7, 2015

kianho commented on May 25, 2015

peterjc commented on May 25, 2015

etal commented on May 27, 2015

lijax commented on May 28, 2015

peterjc commented on May 16, 2017

etal commented on Jul 3, 2017

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions