← Back to branch summary

~ubuntu-branches/ubuntu/natty/python-cogent/natty

« back to all changes in this revision

Viewing changes to cogent/align/pycompare.py

  • Committer: Bazaar Package Importer
  • Author(s): Steffen Moeller
  • Date: 2010-12-04 22:30:35 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20101204223035-j11kinhcrrdgg2p2
Tags: 1.5-1
* Bumped standard to 3.9.1, no changes required.
* New upstream version.
  - major additions to Cookbook
  - added AlleleFreqs attribute to ensembl Variation objects.
  - added getGeneByStableId method to genome objects.
  - added Introns attribute to Transcript objects and an Intron class.
  - added Mann-Whitney test and a Monte-Carlo version
  - exploratory and confirmatory period estimation techniques (suitable for
    symbolic and continuous data)
  - Information theoretic measures (AIC and BIC) added
  - drawing of trees with collapsed nodes
  - progress display indicator support for terminal and GUI apps
  - added parser for illumina HiSeq2000 and GAiix sequence files as 
    cogent.parse.illumina_sequence.MinimalIlluminaSequenceParser.
  - added parser to FASTQ files, one of the output options for illumina's
    workflow, also added cookbook demo.
  - added functionality for parsing of SFF files without the Roche tools in
    cogent.parse.binary_sff
  - thousand fold performance improvement to nmds
  - >10-fold performance improvements to some Table operations

Show diffs side-by-side

added added

removed removed

Lines of Context:
cogent/align/pycompare.py
1
1
#!/usr/bin/env python
2
2
# Very slow.  See compare.pyx
3
3
 
 
4
from __future__ import division
 
5
import cogent.util.progress_display as UI
 
6
from cogent.util.modules import importVersionedModule, ExpectedImportError
 
7
 
4
8
__author__ = "Peter Maxwell"
5
9
__copyright__ = "Copyright 2007-2009, The Cogent Project"
6
10
__credits__ = ["Peter Maxwell", "Gavin Huttley"]
7
11
__license__ = "GPL"
8
 
__version__ = "1.4.1"
 
12
__version__ = "1.5.0"
9
13
__maintainer__ = "Gavin Huttley"
10
14
__email__ = "gavin.huttley@anu.edu.au"
11
15
__status__ = "Production"
12
16
 
13
 
def dotplot(seq1, seq2, window, threshold, min_gap_length=0, band=None, progress=False):
 
17
def py_segments_from_diagonal(seq1, seq2, window, threshold, min_gap_length,
 
18
        diagonal):
 
19
    d_segments = []
 
20
    was_high = False
 
21
    scores = [0] * window
 
22
    score = 0
 
23
    (i_lo, i_hi) = max(0, -diagonal), min(len(seq1), len(seq2)-diagonal)
 
24
    for i in range(i_lo, i_hi):
 
25
        j = i + diagonal
 
26
        k = i % window
 
27
        score -= scores[k]
 
28
        scores[k] = seq1[i] == seq2[j]
 
29
        score += scores[k]
 
30
        if score >= threshold:
 
31
            if not was_high:
 
32
                start = max(i_lo, i - window)
 
33
                if d_segments and start-d_segments[-1][1] < min_gap_length:
 
34
                    (start, jumped_end) = d_segments.pop()
 
35
                was_high = True
 
36
        else:
 
37
            if was_high:
 
38
                d_segments.append((start, i))
 
39
                was_high = False
 
40
    if was_high:
 
41
        d_segments.append((start, i_hi))
 
42
        
 
43
    return d_segments
 
44
 
 
45
try:
 
46
    _compare = importVersionedModule('_compare', globals(),
 
47
            (1, 3), "slow Python dotplot")
 
48
    segments_from_diagonal = _compare.segments_from_diagonal
 
49
except ExpectedImportError:
 
50
    segments_from_diagonal = py_segments_from_diagonal
 
51
 
 
52
@UI.display_wrap
 
53
def dotplot(seq1, seq2, window, threshold, min_gap_length=0, band=None, ui=None):
14
54
    """A list of line segments covering the window-mers with identical matches > threshold
15
55
    
16
56
    Gaps of size less than min_gap will be hidden, which saves on line segments.
17
57
    if 'band' is not None then it limits the searched area
18
58
    """
 
59
    def one_diagonal(dia):
 
60
        segs = segments_from_diagonal(seq1, seq2, window, threshold, 
 
61
                min_gap_length, dia)
 
62
        return [((start, start+dia), (end, end+dia)) for (start, end) in segs]
19
63
    
20
 
    result = []
21
64
    if band is None:
22
65
        band = max(len(seq1), len(seq2))
23
 
    for diagonal in range(-min(len(seq1), band), min(len(seq2), band)+1):
24
 
        d_segments = []
25
 
        was_high = False
26
 
        scores = [0] * window
27
 
        score = 0
28
 
        (i_lo, i_hi) = max(0, -diagonal), min(len(seq1), len(seq2)-diagonal)
29
 
        if progress:
30
 
            if diagonal %100 == 0:
31
 
                print diagonal, (i_lo, i_lo+diagonal), (i_hi, i_hi+diagonal), len(result)
32
 
        for i in range(i_lo, i_hi):
33
 
            j = i + diagonal
34
 
            k = i % window
35
 
            score -= scores[k]
36
 
            scores[k] = seq1[i] == seq2[j]
37
 
            score += scores[k]
38
 
            if score >= threshold:
39
 
                if not was_high:
40
 
                    start = max(i_lo, i - window)
41
 
                    was_high = True
42
 
            else:
43
 
                if was_high:
44
 
                    end = i - 1
45
 
                    if d_segments and start-d_segments[-1][1] < min_gap_length:
46
 
                        (start, jumped_end) = d_segments.pop()
47
 
                    d_segments.append((start, end))
48
 
                    was_high = False
49
 
        if was_high:
50
 
            end = i - 1
51
 
            d_segments.append((start, end))
52
 
        
53
 
        for (start, end) in d_segments:
54
 
            result.append(((start, start+diagonal), (end, end+diagonal)))
 
66
    diagonals = range(-min(len(seq1), band), min(len(seq2), band)+1)
 
67
    result = []
 
68
    for diag_segments in ui.imap(one_diagonal, diagonals, noun='offset'):
 
69
        result.extend(diag_segments)
55
70
    return result