Changeset 310

Show
Ignore:
Timestamp:
05/05/08 10:08:06 (7 months ago)
Author:
t
Message:

added motif searching, better zoom/translation

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/vertex/lib/vertex/data/__init__.py

    r309 r310  
    1414from vertex import coords 
    1515 
    16 genome = maps = colors = None 
    17 refseq_exondb = xeno_refseq_exondb = None 
     16genome = maps = None 
     17refseq_map = conserv_map = None 
     18 
    1819_initialized = False 
    1920 
    2021def load(): 
    21     global genome, maps, colors, _initialized 
    22     global refseq_exondb 
    23     global xeno_refseq_exondb 
     22    global genome, maps, _initialized, refseq_map, conserv_map 
    2423 
    2524    if not _initialized: 
    …… 
    2928        slicedb = shelve.BsdDbShelf(_db) 
    3029 
    31         _db = bsddb.btopen('/mnt/ctb/data/refseq_exons.db', 'r') 
    32         refseq_exondb = shelve.BsdDbShelf(_db) 
     30#        _db = bsddb.btopen('/mnt/ctb/data/refseq_exons.db', 'r') 
     31#        refseq_exondb = shelve.BsdDbShelf(_db) 
    3332 
    3433        _db = bsddb.btopen('/mnt/ctb/data/xeno_refseq.db', 'r') 
    3534        xeno_refseqdb = shelve.BsdDbShelf(_db) 
    3635 
    37         _db = bsddb.btopen('/mnt/ctb/data/xeno_refseq_exons.db', 'r') 
    38         xeno_refseq_exondb = shelve.BsdDbShelf(_db) 
     36#        _db = bsddb.btopen('/mnt/ctb/data/xeno_refseq_exons.db', 'r') 
     37#        xeno_refseq_exondb = shelve.BsdDbShelf(_db) 
    3938 
    4039        _db = bsddb.btopen('/mnt/ctb/data/phastConsMost.db', 'r') 
    …… 
    8180#       maps.append(xeno_refseq_map) 
    8281 
    83         colors = ('red', 'black', 'blue') 
    8482        _initialized = True 

  • trunk/vertex/lib/vertex/web/__init__.py

    r309 r310  
    55    return Publisher(Top()) 
    66 
    7 from . import draw 
    87from . import user 
    98 
    …… 
    1413    _q_exports = ['', 'draw'] 
    1514 
    16     draw = draw.Draw(data.genome, data.maps, data.colors) 
    17      
    1815    def _q_index(self): 
    19         return "hello, world" 
     16        return "Welcome to vertex!<p>Users: <ul><li> <a href='./titus/'>Titus Brown</a></ul>" 
    2017 
    2118    def _q_lookup(self, username): 

  • trunk/vertex/lib/vertex/web/draw/__init__.py

    r309 r310  
    55from vertex import coords 
    66 
     7_cache = {} 
     8 
    79### 
    810 
    …… 
    1012    _q_exports = [''] 
    1113 
    12     def __init__(self, genome, maps, map_colors): 
     14    def __init__(self, cachekey, genome, maps, map_colors): 
     15        self.cachekey = cachekey 
    1316        self.genome = genome 
    1417        self.maps = maps 
    1518        self.map_colors = map_colors 
     19        self.addl_maps = [] 
    1620 
    1721    def _q_index(self): 
    …… 
    2125        (chr, start, stop) = coords.translate(info) 
    2226 
    23         interval = self.genome[chr][start:stop] 
     27#        image = _cache.get((self.cachekey, chr, start, stop)) 
     28        image = None 
     29         
     30        if not image: 
     31            interval = self.genome[chr][start:stop] 
    2432 
    25         picture_class = pygr_draw.BitmapSequencePicture 
    26         p = pygr_draw.draw_annotation_maps(interval, self.maps, 
    27                                            default_colors=self.map_colors, 
    28                                            picture_class=picture_class) 
    29         image = p.finalize() 
     33            picture_class = pygr_draw.BitmapSequencePicture 
     34            p = pygr_draw.draw_annotation_maps(interval, self.maps, 
     35                                               default_colors=self.map_colors, 
     36                                               picture_class=picture_class) 
     37            image = p.finalize() 
     38            _cache[(self.cachekey, chr, start, stop)] = image 
     39        else: 
     40            print '(cache hit: %s)' % (str((chr, start, stop)),) 
    3041 
    3142        response = get_response() 

  • trunk/vertex/lib/vertex/web/user.py

    r309 r310  
     1import quixote 
    12from quixote.directory import Directory 
     3 
     4import motility 
    25 
    36from vertex import data, coords 
    …… 
    1619               'ret-e' : 'chr6:5,894,800-5,895,100' } 
    1720 
     21motifs = { 'WGATAR' : 'purple', 
     22           'AGCGGA' : 'green', 
     23           'GGGTTGC' : 'red' } 
     24 
    1825### 
    1926 
    …… 
    2936 
    3037    def _q_index(self): 
    31         s = "Welcome %s" % (self.username,) 
    32  
     38        s = "Welcome, %s!" % (self.username,) 
     39 
     40        ### 
     41 
     42        keys = views.keys() 
     43        keys.sort() 
     44         
    3345        x = [] 
    34         for k in views: 
     46        for k in keys: 
    3547            x.append("<a href='v/%s'>%s</a>" % (k, k,)) 
    36         s += "<ul><li>" + "<li>".join(x) + "</ul>" 
     48        s += "<p>Bookmarked regions: <ul><li>" + "<li>".join(x) + "</ul>" 
     49 
     50        ### 
     51 
     52        keys = constructs.keys() 
     53        keys.sort() 
     54         
     55        x = [] 
     56        for k in keys: 
     57            x.append("<a href='v/%s'>%s</a>" % (k, k,)) 
     58        s += "Constructs: <ul><li>" + "<li>".join(x) + "</ul>" 
    3759 
    3860        return s 
    …… 
    4466    Translate user viewname into coordinates, and display. 
    4567    """ 
    46     _q_exports = ['draw'] 
     68    _q_exports = ['draw', 'sequence'] 
    4769     
    4870    def __init__(self, username): 
    4971        self.username = username 
    5072 
    51         maps = list(data.maps) 
    52  
     73        maps = [ data.refseq_map, data.conserv_map ] 
     74             
    5375        constructs_nlmsa = make_nlmsa(data.genome, constructs) 
    5476        maps.append(constructs_nlmsa) 
    55          
    56         self.draw = Draw(data.genome, maps, data.colors) 
     77 
     78        maps.append(motif_nlmsa) 
     79 
     80        colors = [ 'red', 'black', 'blue', 'black' ] 
     81         
     82        self.draw = Draw(username, data.genome, maps, colors) 
     83 
     84    def sequence(self): 
     85        request = quixote.get_request() 
     86        where = request.form['i'] 
     87        (chr, start, stop) = coords.translate(where) 
     88 
     89        seq = data.genome[chr][start:stop] 
     90        seq = str(seq) 
     91 
     92        response = quixote.get_response() 
     93        response.set_content_type('text/plain') 
     94        s = ">%s\n" % (where,) 
     95 
     96        for i in range(0, len(seq), 60): 
     97            s += seq[i:i+60] + "\n" 
     98             
     99        return s 
    57100 
    58101    def _q_lookup(self, view): 
    59102        if view in views: 
    60             coords = views[view] 
     103            where = views[view] 
     104        elif view in constructs: 
     105            where = constructs[view] 
    61106        else: 
    62             coords = view 
     107            where = view 
     108 
     109        ### 
     110 
     111        (chr, start, stop) = coords.translate(where) 
     112 
     113        center = (stop + start) / 2 
     114        size = stop - start 
     115 
     116        zoom_out_size = float(size) * 3 
     117        zoom_out = '%s:%d-%d' % (chr, 
     118                                 center - zoom_out_size/2, 
     119                                 center + zoom_out_size/2) 
     120         
     121        zoom_in_size = float(size) / 3 
     122        zoom_in = '%s:%d-%d' % (chr, 
     123                                center - zoom_in_size/2, 
     124                                center + zoom_in_size/2) 
     125 
     126        shift_size = size / 3 
     127        shift_left = '%s:%d-%d' % (chr, start - shift_size, stop - shift_size) 
     128        shift_right = '%s:%d-%d' % (chr, start + shift_size, stop + shift_size) 
     129         
     130        ### 
     131 
     132        x = [] 
     133        for motif, color in motifs.items(): 
     134            x.append('%s (%s)' % (motif, color,)) 
     135        motiflist = ", ".join(x) 
     136 
     137        ### 
    63138 
    64139        return """\ 
    65140<h2>%s</h2> 
    66141%s 
     142<p> 
     143interval size: %s (%s kb) 
     144<p> 
     145<a href='./sequence?i=%s'>get genomic sequence</a> 
     146<p> 
     147<a href='./%s'>zoom out</a>, <a href='./%s'>zoom in</a>, 
     148<a href='./%s'>shift left</a>, <a href='./%s'>shift right</a> 
     149<p> 
     150<a href='..'>return</a> 
    67151<hr> 
    68 <img src='draw/%s'> 
    69 """ % (view, coords, coords) 
     152Display features: <font color='red'>RefSeq mRNA</font>, conservation, 
     153<font color='blue'>constructs</font> 
     154<p> 
     155Motifs: %s 
     156<img src='draw/%s' width='1000'> 
     157""" % (view, where, size, size / 1000, where, 
     158       zoom_out, zoom_in, shift_left, shift_right, motiflist, where) 
    70159 
    71160### 
     161 
     162def make_motif_search_nlmsa(genome, subseqs, motifs): 
     163    from pygr import cnestedlist, seqdb 
     164    from pygr_draw import Annotation 
     165 
     166    d = {} 
     167    m = 0 
     168    for subseq in subseqs: 
     169        assert subseq.start >= 0 
     170 
     171        for motif, color in motifs.items(): 
     172            matches = motility.find_iupac(str(subseq), motif, 
     173                                          offset=subseq.start) 
     174 
     175            for n, (start, stop, orientation, _) in enumerate(matches): 
     176                z = n + m 
     177                d[str(z)] = Annotation(str(z), subseq.id, start, stop, color) 
     178                 
     179            m += len(matches) 
     180 
     181    annodb = seqdb.AnnotationDB(d, genome) 
     182    nlmsa = cnestedlist.NLMSA('simple', mode='memory', use_virtual_lpo=True) 
     183 
     184    for k in annodb: 
     185        nlmsa.addAnnotation(annodb[k]) 
     186 
     187    nlmsa.build() 
     188 
     189    return nlmsa 
    72190 
    73191def make_nlmsa(genome, coords_dict): 
    …… 
    80198        d[str(n)] = Annotation(str(n), chr, start, stop, 'blue') 
    81199 
    82         import sys 
    83         print >>sys.stderr, chr, start, stop 
    84  
    85200    annodb = seqdb.AnnotationDB(d, genome) 
    86201    nlmsa = cnestedlist.NLMSA('simple', mode='memory', use_virtual_lpo=True) 
    …… 
    92207 
    93208    return nlmsa 
     209 
     210data.load() 
     211 
     212x = [] 
     213for v in views.values(): 
     214    (chr, start, stop) = coords.translate(v) 
     215    subseq = data.genome[chr][start:stop] 
     216    x.append(subseq) 
     217 
     218motif_nlmsa = make_motif_search_nlmsa(data.genome, x, motifs) 

  • trunk/vertex/scripts/make-png.py

    r309 r310  
    5656conserv_map = cnestedlist.NLMSA('/mnt/ctb/data/conserv_map', 'r', genomeUnion, 
    5757                              pairwiseMode=True) 
    58 xeno_refseq_map = cnestedlist.NLMSA('/mnt/ctb/data/xeno_refseq_map', 
    59                                     'r', genomeUnion, pairwiseMode=True) 
     58#xeno_refseq_map = cnestedlist.NLMSA('/mnt/ctb/data/xeno_refseq_map', 
     59#                                    'r', genomeUnion, pairwiseMode=True) 
    6060 
    6161sd = refseq_map.seqDict 
    6262 
    63 seq = sd['galGal3.chr6'] 
     63seq = sd['galGal3.chr1'] 
    6464 
    65 subseq = seq[5788957:6042376] 
     65subseq = seq[52958231:53024630] 
    6666 
    6767pprint.pprint(refseq_map[subseq].keys()) 
    …… 
    7272print repr(k[0]) 
    7373 
    74 pprint.pprint([ k.human_name for k in xeno_refseq_map[subseq].keys() ]) 
    75 k = xeno_refseq_map[subseq].keys() 
    76 print k[0].human_name 
     74#pprint.pprint([ k.human_name for k in xeno_refseq_map[subseq].keys() ]) 
     75#k = xeno_refseq_map[subseq].keys() 
     76#print k[0].human_name 
    7777 
    7878k = conserv_map[subseq].keys() 
    …… 
    8181picture_class = pygr_draw.BitmapSequencePicture 
    8282p = pygr_draw.draw_annotation_maps(subseq, 
    83                                    (refseq_map, conserv_map, 
    84                                     xeno_refseq_map), 
     83                                   (refseq_map, conserv_map,), 
    8584                                   picture_class=picture_class) 
     85#                                    xeno_refseq_map), 
    8686 
    8787image = p.finalize()