Pygr: Docs Overview

Previous: 1 Introduction Up: Pygr: Docs Overview Next: About this document ...

---

2 Module Documentation

The following subsections provide details about how to use specific modules of Pygr functionality.

2.1 Installation

tar -xzvf pygr-0.3.tar.gz 
cd pygr
python setup.py install

Once the test framework has completed successfully, the setup script will install pygr into python's respective site-packages directory.

If you don't want to install pygr into your system-wide site-packages, replace the "python setup.py install" command with

python setup.py build

Pygr contains several modules imported as follows:

from pygr import seqdb # IMPORT SEQUENCE DATABASE MODULE

If you did not install pygr in your system-wide site-packages, you must set your PYTHONPATH to the location of your pygr build. For example, if your top-level pygr source directory is PYGRDIR then you'd type something like:

setenv PYTHONPATH PYGRDIR/build/lib.linux-i686-2.3

2.2 graphs and graph query module

The basic idea of Pygr is that all Python data can be viewed as a graph whose nodes are objects and whose edges are object relations (in Python, references from one object to another). This has a number of advantages.

1. All data in a Python program become a database that can be queried through simple but general graph query tools. In many cases the need to write new code for some task can be replaced by a database query.

2. Graph databases are more general and flexible in terms of what they can represent and query than relational databases, which is very important for complex bioinformatics data.

3. Indeed, in Pygr, a query is itself just a graph that can be stored and queried in a database, opening paths to automated query construction.

4. Pygr graphs are fully indexed, making queries about edge relationships (which are often unacceptably slow in relational databases) fast.

5. The interface can be very simple and pythonic: it's just a Mapping. In Python "everything is a dictionary", also known as "the Mapping protocol": a dictionary maps some set of inputs to some set of outputs. e.g. m[a]=b maps a onto b, as a unique relation. In Pygr, if we want to be able to map a node to multiple target nodes (i.e. allow it to have multiple edges), we simply add another layer of mapping: m[a][b]=edgeInfo (where edgeInfo is optional edge info.)

Examples of the Pygr syntax:

graph += node1 # ADD node1 TO graph
graph[node1] += node2 # ADD AN EDGE FROM node1 TO node2
graph[node1][node2]=edge_info # ADD AN EDGE WITH ASSOCIATED edge_info
# ADD SCHEMA BINDING WITH graph[node] BOUND AS node.attr
setschema(node,attr,graph) 
# SEARCH graph FOR SUBGRAPH {1->2; 1->3; 2->3}, 
# I.E. EXONSKIP, WHERE THE SPLICE FROM 2 -> 3 HAS ATTRIBUTE type 'U11/U12' 
for m in GraphQuery(graph,{1:{2:None,3:None},\
                   2:{3:dict(filter=lambda edge,**kwargs:edge.type=='U11/U12')},\
                   3:{}}):
    print m[1].id,m[2].id,m[1,2].id

Let's examine these examples one by one:

• adding a node to a graph is distinct from creating edges between it and other nodes. The graph+=node notation simply adds node to the graph, initially with no edges to other nodes.

• A similar syntax (graph[node1]+=node2) can be used to add an edge between two nodes, but with no edge information. In this case the edge information stored for this relation is simply the Python None value. Note that in Pygr the default type of graph has directed edges; that is a->b does not imply b->a. In the default dictGraph graph class, these are two distinct edges that would have to be added separately if you truly want to have an edge going both from a to b and from b to a.

• To add an edge between two nodes with edge information, use the graph[node1][node2]=edge_info syntax.

• You can bind an object attribute to a graph, using setschema(obj,attr,graph). This acts like Python's built-in setattr(obj,attr,value), but instead of obj.attr simply storing the specified value, it is bound to the graph so that obj.attr is equivalent to graph[obj]. Both syntaxes are interchangeable and can be mixed in different pieces of code accessing the same object.

• Since Pygr adopts the Mapping protocol as its model for storing graphs, you can create graphs simply by creating Python dict objects e.g. foo:bar. In this example we construct a query graph whose "nodes" are just the integers 1, 2, and 3. Since any kind of object is a valid key in Python mappings, they can therefore also be used as "nodes" in a Pygr graph. This query graph illustrates a few simple principles:

• a Pygr graph is just a two-level Python mapping. For example, 1:2,None is a graph with a single edge from 1 to 2, with no edge information. Pygr graphs can have multiple edges from or to a given node.

• edge information in a query graph can be used to specify extra query arguments, again in the form of a Python dictionary. This dictionary is interpreted as a set of "named arguments" to be used by the GraphQuery search method. For example, a filter argument is interpreted as a callable function that is passed a set of named arguments describing the current edge / node matching being tested, and whose return value (True or False) will determine whether this edge "matches" our query graph. In this example, we used it to check whether the edge.type attribute is "U11/U12" (an unusual type of splicing in gene structure graphs).

• Graph query in Pygr simply means finding a subgraph of the datagraph that has node-to-node match to the edge structure given in the query graph. In this example it is a simple exon-skip structure (3 exons, one of which can either be included or skipped). The GraphQuery class provides a general mechanism for performing graph queries on any Python data (see below for full details). It can be used as an iterator that will return all matches to the query (if any).

• Matches are themselves returned as a mapping of nodes and edges of the query graph (in this example, its nodes are the integers 1, 2 and 3) onto nodes and edges of the data graph. In this example the match is returned as m, so m[1] is the node in the data graph corresponding to node 1 in the query graph. This example assumes that object has an id attribute, which is printed out. To refer to an edge, just use a tuple corresponding to a pair of nodes in the query graph. In this example, 1,2 refers to the edge from node 1 to node 2 in the query graph, so m[1,2] is the edge in data graph between nodes m[1] and m[2]. This example also attempts to print an id attribute from that edge object.

• Note on current behavior: currently, GraphQuery will throw a KeyError exception if it tries to search for a query node in the query graph and does not find it. That's why we have to add the "node with no edges" entry 3: for node 3. This will probably be addressed in the future, since this seems like a potential source of many annoying little bugs.

2.2.1 dictGraph

dictGraph is Pygr's main graph class. It provides all the standard behaviors described above. The current reference implementation uses standard Python dict objects to store the graph. All the usual Mapping protocol methods can be used on dictGraph objects (top-level interface, in the examples above graph) and dictEdge objects (second-level interface; in the examples above graph[node]). e.g.

• for node in graph: iterator method returns all nodes in the graph; you could also use graph.items() to get node,dictEdge pairs, etc.

• for node in graph[node]: iterator method returns all nodes that are targets of edges originating at node. Again, you could use graph[node].items() to get node,edgeInfo pairs. Note: if node is not in graph, this will throw a KeyError exception just like any regular Python dict.

• if node in graph: contains method checks whether node is present in the graph, using dict indexing.

• if node2 in graph[node1]: test whether node1 has an edge to node2. Again, if node1 isn't in graph, this will throw a KeyError exception.

2.2.2 Directionality and Reverse Traversal

Note that dictGraph stores directed edges, that is, a->b does not imply b->a; those are two distinct edges that would have to be added separately if you want an edge going both directions. Moreover, the current implementation of dictGraph does not provide a mechanism for traveling an edge backwards. To do so with algorithmic efficiency requires storing each edge twice: once in a forward index and once in a reverse index. Since that doubles the memory requirements for storing a graph, the default dictGraph class does not do this. If you want such a "forward-backwards" graph, use the dictGraphFB subclass that stores both forwad and reverse indexes, and supports the inverse operator (). graph gets the reverse mapping, e.g. ( graph)[node2] corresponds to the set of nodes that have edges to node2. This area of the code hasn't been tested much yet.

2.2.3 Schema: binding object attributes to graphs

The goal of Pygr is to provide a single consistent model for working with data explicitly modeled as graphs (i.e. dictGraph-like objects) and standard Python objects that were not originally designed to be queried (or thought of) as a "graph". Since Python uses the Mapping concept throughout the language and object model, and provides introspection, there is no reason why Pygr can't work with both kinds of data transparently. One mechanism for making this idea explicit is the idea of binding an object attribute to a graph, via the new method we've called setschema(obj,attr,graph). The idea here is that once you bind an object attribute to a graph, the two different data models obj.attr (object model) or graph[obj] (graph model) are made equivalent and interchangeable. Operating on one affects the other and vice versa; they are two ways of referring to the same relation. This concept can be applied at several different levels

• individual objects: just like getattr() and setattr(), you can apply schema methods to individual objects: getschema(obj,attr) (returns the bound graph) or setschema(obj,attr,graph) (binds the object attribute to the graph).

• all instances of a class: you can bind specific attributes of a given class to a graph using the following class attribute syntax:

class ExonForm(object): # ADD ATTRIBUTES STORING SCHEMA INFO
    __class_schema__=SchemaDict(((spliceGraph,'next'),(alt5Graph,'alt5'),(alt3Graph,'alt3')))

In this class we bound the next attribute to spliceGraph, alt5 attribute to alt5Graph, and alt3 attribute to alt3Graph. That means, every instance obj of this class will have an attribute obj.next that is equivalent to spliceGraph[obj], etc. Note that this is schema, not the actual operation of adding the object as a node to the graph. Indeed, when obj is first created, it is not automatically added to spliceGraph; that is up to the user. Unless your code has added the node to the graph (e.g. spliceGraph+=obj), obj.next should throw a KeyError exception.

The general method getschema(obj,attr) works regardless of whether the schema was stored on an individual object or at the class level.

2.2.4 GraphQuery

The GraphQuery class implements simple node-to-node matching, in which each new node-set is generated by an iterator associated with a specific node in the query graph. This iterator model is general: since indexes (mappings) support the iterator protocol, a given iterator may actually be an index lookup (or other clever search algorithm). The GraphQuery constructor takes two arguments: the default data graph being queried, and the query graph. The query graph is just a graph; its nodes can be any object that can be a graph node (i.e. any object that is indexible, e.g. by adding a __hash__() method). Its node objects will not be modified in any way by the GraphQuery. Its edges are expected to be dictionaries that can be checked for specific keyword arguments:

• filter: must be a callable function that accepts keyword arguments and returns True (accept this edge as a match to the queryGraph) or False (do not accept this edge as a match). This function will be called with the following keyword arguments:
  • toNode: the target node of this edge, in the data graph
  • fromNode: the origin node of this edge, in the data graph
  • edge: the edge information for this edge in the data graph
  • queryMatch: a mapping of the query graph to the data graph, based on the partial matchings made so far
  • gqi: the GraphQueryIterator instance associated with this matching operation. Much more data is available from specific attributes of this object.

• dataGraph: graph in which the current edge should be search for. This allows a query to traverse multiple graphs. In other words, when searching for edges from the current node, look up dataGraph[node] instead of defaultGraph[node].

• attr: object attribute name to use as the iterator, instead of the defaultGraph.In other words, generate edges from the current node via getattr(node,attr) instead of defaultGraph[node]. The object obtained from this attribute must act like a mapping; specifically, it must provide an items() method that returns zero or more pairs of targetNode,edgeInfo, just like a standard Pygr dictEdge object.

• attrN: object attribute name to use as the iterator, instead of the defaultGraph. In other words, generate edges from the current node via getattr(node,attr) instead of defaultGraph[node]. The object obtained from this attribute must act like a sequence; specifically, it must provide an iterator that returns zero or more targetNode. The edgeInfo for any edges generated this way will be None.

• f: a callable function that must return an iterator producing zero or more pairs of targetNode,edgeInfo. Typically f is a Python generator function containing a statement like yield targetNode,edgeInfo.

• fN: a callable function that must return an iterator producing zero or more targetNode. Typically fN is a Python generator function containing a statement like yield targetNode. The edgeInfo for any edges generated this way will be None.

• subqueries: a tuple of query graphs to be performed. Since GraphQuery traversalcorresponds to logical AND (i.e. all the query graph nodes must be successfully matched to return a match), the subqueries are currently treated as a union (logical OR), by simply returning every match from each subquery as a match (at least for this node). Each subquery is itself just another query graph. Moreover, since query graphs can share nodes (i.e. the same object can appear as a node in multiple query graphs), subqueries can make reference to nodes that are already matched by the higher query. This is an area that has not been explored much yet, but provides a pretty general model for powerful queries.

Note: an easy way to pass keyword dictionaries (e.g. as edge information) is simply using the dict() constructor, e.g. dict(dataGraph=myGraph,filter=my_filter). I think this is a little more readable than 'dataGraph':myGraph, 'filter':my_filter.

Note on current behavior: currently, the GraphQuery iterator returns the same mapping object for each iteration (simply changing its contents). So to save these multiple values safely in a list comprehension we have to copy each one into a new dict object via dict(m).

2.2.5 What is GraphQuery actually doing?

A GraphQuery is basically an iterator that returns all possible mappings of the query graph onto the datagraph that match all of the nodes and edges of the query graph onto nodes and edges of the data graph. As an iterator, it does not instantiate a list of the matches, but simply returns the matches one by one. The current design is very simple. The GraphQuery constructor builds an "iterator stack" of GraphQueryIterators, each representing one node in the query graph; they are enumerated in order by a breadth-first-search of the query graph. The GraphQuery iterator processes the stack of GraphQueryIterators: any match simply pushes the stack to the next level; any match at the deepest level of the stack is a complete match (yield the queryMatch mapping); the end of any GraphQueryIterator simply pops the stack. One obvious idea for improving all this is to replace this "interpreter" with a "compiler" that compiles Python for loops that are equivalent to this stack, and run that... likely to be many fold faster.

2.3 sequence Module

Base classes for representing sequences and sequence intervals.

2.3.1 Overview

2.3.2 SeqPath

• Python Sequence: of course, SeqPath behaves like a Python sequence. i.e. the length of a SeqPath s is just len(s), and you iterate over the ``letters'' in it using for l in s: (Note, the individual letters produced by this iterator will themselves be SeqPath objects (by default, of length 1)). And all the slicing operations defined for Python Sequences also apply to SeqPath (see below).

• Slicing: SeqPath is designed to represent a slice (subinterval) of a sequence. Like the Python builtin slice class, it has start, stop, and step attributes that indicate the interval beginning, end, and ``stride''. Moreover, it is itself sliceable in the usual pythonic way, i.e. s[start:stop], where start and stop are in the local coordinate system of s (i.e. s[0] is the first letter of the interval represented by s). Note that SeqPath follows the Python slicing coordinate conventions of positive integers as forward coordinates (i.e. counting from the interval start) and negative integers as reverse coordinates (i.e. counting from the interval end).

• String value: to obtain the actual sequence string representation of a SeqPath, just use the Python builtin str(s). Note that in most cases a SeqPath object does not itself store the sequence string associated with it but obtains it from somewhere else when the user requests it.

• comparison and containment: SeqPath implements the interval-ordering and interval-containment relations using the standard Python order operators and containment operators. i.e. s<t iff s.start<t.start, and s in t iff t.start<=s.start and s.stop<=t.stop.

• orientation: SeqPath carefully represents relationships between intervals on opposite strands of a double-stranded nucleotide sequence. A SeqPath object knows whether it is an interval on the forward or reverse strand. Pygr provides a number of operations for manipulating and comparing intervals of different orientations. For example, -s yields the interval of the opposite strand that is base-paired to interval s (i.e. this is not just the reverse-complement of s in the string 'atgc' 'gcat' sense, but is specifically the SeqPath object representing the coordinate interval on the opposite strand that is base-paired to s).

• schema: a SeqPath object knows ``what sequence'' it is an interval of; it is not just a (start,stop) coordinate pair, but is actually bound to a specific parent sequence object. Specifically, s.path is the parent sequence object of which s is a subinterval; s.path will itself be an instance of SeqPath, and its path attribute will simply be itself. All SeqPath objects are descended from such ``top-level'' SeqPath objects. Note that when you have sequence intervals from both forward and reverse strands of a sequence, all of the forward strand intervals will share the same path attribute (your original top-level sequence object representing the whole sequence in forward orientation), while all the reverse strand intervals will reference another top-level SeqPath created automatically to represent the reverse strand.

• graph structure: a SeqPath object itself acts as a graph, whose nodes are the individual letters of the sequence, and whose edges represent the link from each letter to the next (if any). Thus standard graph query works on SeqPath objects, through the usual interfaces:

  for l in s: # GET EACH LETTER OF THE SEQUENCE
      c=str(l)
  
  edge=s[l1][l2] # GET EDGE INFORMATION FOR l1 --> l2
  
  for l1,l2,edge in s.edges(): # GET ALL l1 --> l2 EDGES
      do_something(l1,l2,e)
  
  # DUMB GRAPH QUERY TO FIND 'AG' SUBSTRINGS IN SEQUENCE s
  for d in GraphQuery(s,{0:{1:dict(filter=lambda fromNode,toNode:
                                   str(fromNode)=='A' and str(toNode)=='G')},
                         1:{}})
      l1,l2,edge=d[0],d[1],d[0,1]
  

  For more information about edges, see the LetterEdge class.

• Mutable Sequences: Just as the Python builtin list class implements ``mutable sequence'' objects that can be resized, SeqPath objects can be resized and changed, without breaking existing subinterval objects that are ``part of'' the resized SeqPath object. In particular, just as a list can be resized by extending its ``stop'' coordinate to a higher value, a SeqPath can be resized by extending its stop coordinate to a higher value. Indeed, you can even create a SeqPath for a particular sequence without knowing that sequence's length (computing the length of a genome sequence might take a long time, if all you want to do is create a sequence object to represent that sequence). You can do this by passing None as the stop (or start) coordinate. In that case, SeqPath will automatically determine its own length at a later time iff a specific user operation makes it absolutely necessary to know this length.

• intersection, union, difference: SeqPath uses the Python *, + and - operators to implement interval intersection, union, and difference operations respectively.

2.3.3 Sequence class

sequence.Sequence(

s, id)

The Sequence class provides a SeqPath flavor that stores a sequence string s and identifier id for this sequence.

from pygr import sequence
seq=sequence.Sequence('GPTPCDLMETQ','FOOG_HUMAN')

update(

s)

You can change the actual string sequence to a new string s using the update method:

seq.update('TKRRPLEDKMNEPS')

seqtype(

)

returns DNA_SEQTYPE for DNA sequences, RNA_SEQTYPE for RNA, and PROTEIN_SEQTYPE for protein.

reverse_complement(

s)

returns the reverse complement of the sequence string s.

• id: the id attribute stores the sequence's identifier.

2.3.4 Coordinate System

Each SeqPath object has a number of attributes giving information about its ``location'':

• orientation: +1 if on the forward strand, or -1 if on the reverse strand.

• path: the top-level sequence object that this interval is part of, or self if this object is its top-level (i.e. not a slice of a larger sequence). Note that all forward intervals share the same path attribute, but reverse strand intervals all have a path attribute that represents the entire reverse strand.

• pathForward: same as path, but always the forward strand sequence.

• start: start coordinate of the interval. NB: SeqPath stores coordinates relative to the start of the forward strand. This is necessary for allowing resizing of the top-level SeqPath; if coordinates were relative to the end of the sequence, they would have to be recomputed every time the length of the sequence changed. The main consequence of this is that coordinates for forward intervals are always positive, whereas coordinates for reverse intervals are always negative (i.e. following the Python convention that negative coordinates count backwards from the end, and the fact that the end of the reverse strand corresponds to the start of the forward strand). NB2: if the SeqPath was originally created with start=None, requesting its start attribute will force it to compute its start coordinate, typically requiring a computation of the sequence length. In this case, the start attribute will computed automatically by SeqPath.__getattr__().

• stop: end coordinate of the interval. The above comments for start apply to stop. Note that for reverse intervals, a stop value of 0 means the end of the reverse strand (i.e. -1 is the last nucleotide of the reverse strand, and 0 is one beyond the last nucleotide of the reverse strand).

• _abs_interval: a tuple giving the (start,stop) coordinates of the interval on the forward strand corresponding to this interval (i.e. for a forward interval, itself, or for a reverse interval, the interval that base-pairs to it).

2.3.5 Extending and Customizing

strslice(

start, stop)

called to get the string sequence of the interval (start, stop). You can provide your own strslice() method to customize how sequence is stored and accessed. For example, SQLSequence.strslice() gets the sequence via a SQL query, and BlastSequence.strslice() obtains it using the fastacmd -L start,stop UNIX shell command from the NCBI toolkit.

__len__(

)

called to compute the length of the sequence. You can customize this to provide an efficient length method for your particular sequence storage. e.g. SQLSequence obtains it via a SQL query; BlastSequence obtains it from a precomputed length index. The default Sequence.__len__() method computes it from len(self.seq), assuming that the sequence can be accessed from the seq attribute.

__getitem__(

slice_obj)

if you want to monitor or intercept slicing requests on your sequence, you can do so by providing your own getitem method. See seqdb.BlastSequenceCache class for an example.

__getattr__(

attr)

if you subclass a SeqPath-derived class and supply a __getattr__ method for your subclass, it must call the parent class's __getattr__. This is essential for ``delayed evaluation'' of start and stop attributes, which are generated automatically by SeqPath's __getattr__. If your subclass inherits from more than one parent class, check whether both parents supply a __getattr__, in which case your subclass must supply a __getattr__ that explicitly calls both of them. Failing to do so will lead to strange bugs.

• seq: the Sequence.strslice() method assumes that the actual sequence is stored on the seq attribute. You could customize this behavior by making the seq attribute a property that is computed on the fly by some method of your own.

2.3.6 IntervalTransform

xform=IntervalTransform(srcPath,destPath)
d2=xform(s2) # MAPS s2 FROM srcPath coords to destPath coord system
d3=xform[s2] # CLIPS s2 TO NOT EXTEND OUTSIDE srcPath, THEN XFORMS
s3=xform.reverse(d3) # MAP BACK TO srcPath COORD SYSTEM

2.3.7 Seq2SeqEdge

__init__(

msaSlice, targetPath, sourcePath=None)

Create a Seq2SeqEdge for the targetPath, on the specified alignment slice. If sourcePath is None, it will be calculated automatically by calling the slice's methods.

__iter__(

sourceOnly=True, **kwargs)

iterate over source intervals within this segment of alignment. kwargs will be passed on to the msaSlice's groupByIntervals and groupBySequences methods.

items(

**kwargs)

same as __iter__, but gets tuples of (source_interval,target_interval).

pIdentity(

mode=max,trapOverflow=True)

Compute the percent identity between the source and target sequence intervals in this segment of the alignment. mode controls the method used for determining the denominator based on the lengths of the two aligned sequence intervals. trapOverflow controls whether overflow (due to multiple mappings of the query sequence to different regions of the alignment) is trapped as an error. To turn off such error trapping, set trapOverflow=False.

pAligned(

mode=max,trapOverflow=True)

Compute the percent alignment between the source and target sequence intervals in this segment of the alignment, i.e. the fraction of residues that are actually aligned as opposed to gaps / insertions, in the two intervals.

Warning: if your query sequence has multiple mappings in the alignment (i.e. it is aligned to two or more different regions in the alignment), pIdentity() and pAligned() may return fractions larger than 1.0. This is because the query sequence may align to a given target sequence via more than one region in the alignment. If you encounter this problem, you can iterate through the individual mappings yourself (by calling the iter(), items() or edges() iterator methods for your alignment slice object), and calculating the percentage identity or alignment (via your own algorithm) individually for each specific mapping. For more background on this problem, see ``Multiple Mappings'', below.

Note that the presence of multiple mappings is not a Pygr bug, but simply reflects the alignment data loaded into Pygr; if the input alignment contains multiple mappings, there is nothing that Pygr can do to make this problem magically vanish.

2.3.8 SeqFilterDict

d=SeqFilterDict(seqIntervalList)
overlap=d[ival] # RETURNS INTERSECTION OF ival AND STORED IVAL, OR KeyError

You can pass a list of intervals to store to the class constructor (as shown above). You can also add a single interval using the syntax d[saveInterval]=saveInterval. (This syntax reflects the actual mapping that the dictionary will perform if later called with the same interval).

2.3.9 LetterEdge

__iter__(

)

iterate over seqpos for sequences that traverse this edge.

iteritems(

)

generate origin, target seqpos for sequences that traverse this edge.

__getitem__(

seq)

return origin,target seqpos for sequence seq; raise KeyError if not in this edge

• seqs: returns its sequences that traverse this edge

2.3.10 Functions

guess_seqtype(

s)

based on the letter composition of the string s, returns DNA_SEQTYPE for DNA sequences, RNA_SEQTYPE for RNA, and PROTEIN_SEQTYPE for protein.

absoluteSlice(

seq, start, stop)

returns the sequence interval of top-level sequence object associated with seq, interpreting start and stop according to the Pygr convention: a pair of positive values represents an interval on the forward strand; a pair of negative values represents an interval on the reverse strand (see Coordinate System, above).

2.4 Alignment Module

Pygr interface to sequence alignment, and scalable storage for multigenome alignments

2.4.1 Overview

Pygr provides a general model for interfacing with any kind of sequence alignment, and also a uniquely scalable storage system for working with huge multiple sequence alignments such as multigenome alignments. Specifically, it lets you work with an alignment both in the traditional Row-Column model (each row is a sequence, each column is a set of individual letters from different sequences, that are aligned; we will refer to this as the RC-MSA model), and also as a graph structure (known as a Partial Order Alignment, which we will refer to as the PO-MSA model). This supports ``traditional'' alignment analysis, as well as graph-algorithms, and even graph query of alignments.

This model has a few basic classes:

• MSA: this class represents an entire alignment. It acts as a graph whose nodes are sequences (or sequence intervals) that are aligned, and whose edges represent specific alignment relationships between specific pairs of sequences (or intervals). Specifically, it acts as a dictionary whose keys are SeqPath objects, and whose values are MSASlice objects (representing an alignment segment associated with a specific SeqPath, see below for details). For example, to find out what's aligned to some sequence interval s1:

  for s2 in msa[s1]: # GET ALL INTERVALS s2 ALIGNED TO s1 IN msa
      do_something(s1,s2)
  

  In addition, its letters attribute acts as a graph interface to the Partial Order alignment (PO-MSA) representation of the alignment. I.e. it is a graph whose nodes each represent a set of individual letters from different sequences, that are aligned to each other, and whose edges connect pairs of nodes that are ``adjacent'' to each other in at least one sequence. Specifically, it acts as a dictionary whose keys are MSANode objects (see below), and whose edges are LetterEdge objects (see previous section).

  for node in msa.letters: # GET ALL ALIGNMENT ``COLUMNS'' IN msa
      for l in node: # GET ALL INDIVIDUAL SEQ LETTERS ALIGNED HERE
          say_something(node,l)
  

• MSASlice: this class represents a segment of alignment associated with a specific sequence interval (s1). It acts as dictionary whose keys are sequence intervals s2 aligned to s1, and whose values are MSASeqEdge objects that represent the alignment relationship between s1 s2. It also has a letters attribute, that represents the subgraph of nodes associated specifically with s1, and the edges that interconnect them.

  myslice=msa[s1]: # GET SLICE ALIGNED TO s1 IN msa
      for node in myslice.letters:  # GET ALL ALIGNMENT ``COLUMNS'' FOR s1
          for l1,l2,e in node.edges(): # GET INDIVIDUAL LETTERS ALIGNED TO l1 OF s1
  	    whatever(l1,l2,e)
  

  This class also has a regions method that generates all the alignment interval relationships in this slice according to ``grouping'' criteria such as maximum permissible gap length, etc. (i.e. any region of alignment containing no gaps larger than a specified size would be returned as a single region, whereas any gap larger than the specified size would split it into two separate regions). This provides a general interface for group-by operations in alignment query.

• MSASeqEdge: this class represents a relationship between a pair of sequence intervals s1 and s2 (SeqPath objects). It provides a mapping between subintervals of s1 s2. I.e. it acts as a dictionary that accepts subintervals of s1 as keys, and maps them to aligned subintervals of s2. It also has a letters attribute, that represents the subgraph of nodes associated specifically with them, and the edges that interconnect the nodes.

• MSANode: this class represents a specific ``column'' in the alignment that aligns a set of individual letters from different sequences. This corresponds to a node in the PO-MSA representation of the alignment. It acts as a dictionary whose keys are sequence intervals (typically only one letter long) aligned in this column, and whose values are MSASeqEdges representing the alignment of that letter to the column (see above).

2.4.2 NestedList Storage

The nested list data structure solves this problem, by moving any interval in the database that is contained in another interval out of the top-level interval list, into the sublist of the parent interval. Based on this, one can prove that one can stop the scanning operation at the first non-overlapping interval (i.e. the overlapping intervals in any list form a single contiguous block). Overall, this reduces the query time to , where is the number of intervals in the database that actually overlap the query (i.e. results to return). Moreover, the nested list data structure can be implemented very well both in computer memory (RAM) or as indexed disk files. Pygr's disk-based cnestedlist database can complete a typical interval query of the 26GB UCSC 8 genome alignment in about 60 microseconds, compared with 10-30 seconds per query using MySQL.

2.4.3 Multiple Mappings: a Warning

• It introduces ambiguity in the alignment: you don't know which of the multiple hits is considered to be the ``right'' alignment; the UCSC alignment file does not tell you.

• There is no scoring information to resolve this ambiguity. In a way, this situation is even worse than the common situation we previously faced in search for alignment mappings using BLAST, because (unlike BLAST) the MAF alignment does not give a score that indicates which mapping is best. (We haven't seen such scoring information; if it can be recovered for these alignment files, we'd be love to know about that...).

• It can cause ``buggy'' results in calculations based on the alignment. For example, Pygr's pIdentity() and pAligned() computations can give values larger than one when a query region has multiple hits. This is not, strictly speaking, a Pygr bug: the query region is mapped by the MAF file to the same target region multiple times, resulting in multiple overlaps.

If you encounter multiple mappings, you can always iterate over them one by one, and perform your own computations for each one. However, to avoid them altogether, you can restrict your queries to the reference genome for this specific alignment (UCSC offers different versions of each alignment set, each based on a different reference genome).

2.4.4 NLMSA

• PO-MSA: Conceptually, the alignment is represented as a partial order alignment (PO-MSA), in which aligned sequence intervals are fused together as a single ``node'' in the alignment graph; two nodes are connected by an edge if and only if they are adjacent in at least one of the sequences aligned to them (i.e. if residue i of that sequence is in the first node, and residue i+1 is in the second node, then there is a directed edge from the first PO-MSA node to the second node).

• LPO: This alignment graph is partially ordered. Let's define an ordering relation ``i<j'' to mean ``there exists a path of directed edges from i to j''. For two letters i and j in a sequence, i<j XOR j<i (i.e. all nodes have an ordering relationship). By contrast, if two nodes in the LPO represent insertions in different sequences, then NOT i<j AND NOT j<i. Thus there can be some nodes in the LPO that have no ordering relationship with respect to each other. It is still possible to map the PO-MSA onto a linear coordinate system (i.e. to ``linearize'' the partial order): as long as the graph contains no cycles, we can map the nodes i,j,k,... of the graph onto a linear coordinate system x,y,z,... such that for any pair of nodes i,j mapped to coordinates x<y, we assert NOT j<i. This is called the linearized partial order (LPO). This maps the PO-MSA onto a standard Row-Column MSA format, where the LPO coordinate (just an integer sequence 0,1,2...) can be considered the index value of each alignment column.

• nested list: The actual alignment data are stored in the form of (start,stop) pairs representing aligned intervals. Since this representation uses intervals, not individual letters, it takes no more memory to store an alignment of two 100 kb regions than it does to align two individual letters. This is important for scalable storage (and query) of large multi-genome alignments. (Each alignment interval takes 24 bytes: five int for the (start,stop) pairs and target sequence ID, plus one int for the sublist ID). These interval databases are stored using nested lists. Specifically, the alignment is stored as 1) a mapping of each aligned sequence interval onto an LPO coordinate interval; 2) a reverse mapping of each LPO interval onto all the sequence intervals that are aligned there. To find the alignment of a sequence interval onto the other sequences in the alignment, that interval is first mapped onto the LPO, and from there mapped back to intervals in the other sequences. A nested list database is stored for each of these mappings (i.e. for an alignment of N sequences, there will be N+1 nested list databases to store the MSA). Furthermore, if the size of the LPO coordinate system (i.e. number of columns in its RC-MSA format) grows larger than the range representable by int (typically = 2 GB), the LPO will have to be split into separate nested list databases of a size smaller than the maximum range representable by int. This is necessary for handling alignments of large genomes (e.g. the human genome is approximately 3 GB). Pygr takes care of all this for you automatically. Note, as an entirely separate issue, that Pygr's cnestedlist module uses the long long data type for file offsets and the fseeko() POSIX interface for large file support (i.e. 64-bit file sizes), which is supported by current versions of Linux, Mac OS X, etc; otherwise, check if your filesystem supports this.

This functionality is encapsulated in the NLMSA class, which has a number of methods and attributes.

Construction Methods:

NLMSA(

pathstem='', mode='r', seqDict=None, mafFiles=None, maxOpenFiles=1024, maxlen=None, nPad=1000000, maxint=41666666)

Constructor for the class. pathstem specifies a path and filename prefix for the NLMSA files (since multiple files are used to store one NLMSA, it will automatically add a number of suffixes automatically to open the necessary set of files for the NLMSA). mode is either ``r'' to open an existing NLMSA (from the pathstem disk files), or ``w'' to create a new one (which will be saved to the pathstem disk files). seqDict specifies a dictionary which maps sequence names to actual sequence objects representing those sequences. If seqDict is None, the constructor will look for a file pathstem.seqDict that is a prefixUnionDict header file for opening all the sequence database files for you automatically.

mafFiles can be used to specify a list of filenames storing a multiple sequence alignment in the UCSC MAF format. maxOpenFiles limits the open file descriptors the NLMSA will use. Since each sequence has a separate nested list database file, a large multi-genome alignment (with each genome containing 20 different chromosomes, say) can rapidly open a large number of file descriptors. Note: NLMSA only opens a given sequence's nested list database when one of your queries actually requires access to that sequence; it then keeps that file descriptor open to make subsequent queries to it fast. If the number of open file descriptors would exceed maxOpenFiles, it will close other open database files, which may slow down query performance (due to having to open and close databases repeatedly to process queries). maxlen specifies the maximum coordinate value for an LPO database. Its default value is 2GB, to prevent int overflow. Using a smaller value can be useful, to 1) limit the size of the LPO in memory during initial construction, and 2) to limit the size of LPO database files on disk (if for example, your file system does not support files above some maximum size). During initial construction of the NLMSA (from MAF files or user-specified interval alignments), the algorithm performs a one-pass sort of the LPO intervals. Thus, this set of intervals is briefly held in RAM for this sort. If you have insufficient RAM, the construction step may raise a MemoryError. If so, you can avoid this problem by using a smaller maxlen value. The maxint option provides another way of limiting the size of LPO databases. It specifies the maximum number of intervals to store per database. Since each interval takes 24 bytes, the default setting limits each LPO to a total size of 1 GB. Note that the current NLMSA construction algorithm requires loading each database index into memory as one-time operation during construction. If your NLMSA build fails due to running out of memory, simply reduce this value. The nPad option sets the maximum number of LPO coordinate systems (specifically, the offset for the start of real sequence IDs in the NLMSA sequence index). You are unlikely to need to change this default value.

__setitem__(

i1, s1)

store an alignment interval; specifically, align the sequence interval s1 to the LPO coordinate interval i1. This provides a trivial user interface for you to build any desired alignment. If the sequence containing interval s1 is not already in the NLMSA, it will be added for you automatically (i.e. creating the necessary indexing, nested list database files, etc.). In this case, the sequence must supply a unique string identifier, which will be used on subsequent attempts to open the NLMSA database, to match the individual sequence nested-list databases against corresponding sequence objects (using seqDict, see above).

build(

)

to construct the final nested list databases, after all the desired alignment intervals have been saved (using either the mafFiles constructor argument, or using setitem above). This method simply calls the build() method on all the constituent NLMSASequence objects in this alignment.

Alignment Usage Methods:

__getitem__(

s1)

get the alignment slice for the sequence interval s1, i.e. get an NLMSASlice object representing the set of intervals aligned to s1. You can also use a regular Python slice object using integer indices ie. nlmsa[1:45], in which case, it gets the NLMSA slice corresponding to that region of the LPO coordinate system.

• seqs: This attribute provides a dictionary of the sequences in the NLMSA, whose keys are top-level sequence objects, and whose values are the associated NLMSASequence object for each sequence. Ordinarily you will have no need to access the NLMSASequence object directly; only do so if you know what you're doing (details below). This dictionary is of type NLMSASeqDict (see below).

2.4.5 NLMSASlice

__iter__(

)

iterates over all sequence intervals that have a 1:1 mapping (i.e. a block of alignment containing no indels) to all or part of the source interval.

keys(

maxgap=0, maxinsert=0, mininsert= 0, filterSeqs=None, mergeAll=False, ivalMethod=None, sourceOnly=False, indelCut=False, seqGroups=None, minAligned=1, pMinAligned=0., seqMethod=None, **kwargs)

Provides a more general interface than iter(), with two types of group-by capabilities, ``group-by'' operations on the alignment intervals contained within this slice (``horizontal'' grouping), and on the sets of sequences aligned to this slice (``vertical'' grouping).

1. ``group-by'' operations on the alignment intervals contained within this slice. It allows the user to supply various parameters for controlling when alignment intervals will be merged or split in the results that it returns. mergeAll forces it to combine intervals of a given sequence irrespective of the size of gaps or inserts separating them. maxgap sets the maximum gap size for merging two adjacent intervals. If the target sequence for the two alignment intervals has a gap longer than maxgap letters between the two alignment intervals, they will be returned as separate intervals; otherwise they will be merged as a single alignment region. maxinsert sets the maximum length of insert in the target sequence that allows to adjacent intervals to be merged as a single alignment region in the results. mininsert is specifically for handling alignments that may have small ``cycles'' (due to slight inconsistencies in the reported alignment intervals, for example, if a portion of sequence can align at both the end of one interval or at the beginning of another, and the intervals are actually added to the NLMSA that way, then the start of the second interval will actually be before the stop of the first interval; this corresponds to a negative insert value). A mininsert value of zero (the default), prevents any such interval pairs from being merged. Giving a negative mininsert value will allow interval pairs whose insert value is greater than or equal to this value, to be merged. filterSeqs, if not None, should be a dict of sequences used to filter the group-by analysis; i.e. only alignment intervals containing these sequences are considered in the analysis. More specifically, filterSeqs can be used to mask the group-by analysis to a specific interval of a sequence, by having filterSeqs return only the intersection between the interval it is passed as a key, and the masking interval that it stores. If there is no overlap, it must raise KeyError. The sequence.SeqFilterDict class provides exactly this masking capability, i.e.

d=sequence.SeqFilterDict(someIntervals)
overlap=d[ival] # RETURNS INTERSECTION BETWEEN ival AND someIntervals, OR KeyError

ivalMethod, if not None, allows the user to provide a Python function that performs interval grouping. Specifically it is called as ivalMethod(l, ns,msaSlice=self, **kwargs), where l is the list of intervals for NLMSASequence ns within the current slice msaSlice; all other args are passed as a dict in kwargs.

2. merge groups of sequences using "vertical" group-by rules. seqGroups: a list of one or more lists of sequences to group. If None, the whole set of sequences will be treated as a single group. Each group will be analyzed separately, as follows: sourceOnly: output intervals will be reported giving only the corresponding interval on the source sequence; redundant output intervals (mapping to the same source interval) are culled. Has the effect of giving a single interval traversal of each group. indelCut: for sourceOnly mode, do not merge separate intervals that the groupByIntervals analysis separated due to an indel). minAligned: the minimum number of sequences that must be aligned to the source sequence for masking the output. Regions below this threshold are masked out; no intervals will be reported in these regions. pMinAligned: the minimum fraction of sequences (out of the total in the group) that must be aligned to the source sequence for masking the output. seqMethod: you may supply your own function for grouping. Called as seqMethod(bounds,seqs,**kwargs), where bounds is a sorted list of (ipos,isStart,i,ns,isIndel,(start,end,targetStart,targetEnd)) and seqs is a list of sequences in the group. Must return a list of (sourceIval,targetIval). See the docs.

iteritems(

**kwargs)

same keys as iter, but for each provides the source interval to target interval mapping (Seq2SeqEdge). Uses same group-by arguments as keys().

edges(

**kwargs)

same interval mappings as iteritems, but for each provides a tuple of three objects: the source interval, the corresponding target interval, and the Seq2SeqEdge providing detailed information about the alignment between the source and target intervals (such as percent identity, etc.). Uses same group-by arguments as keys().

__getitem__(

s1)

treats s1 as a key (target sequence interval), and returns an Seq2SeqEdge object providing detailed information about the alignment between this target interval and the source interval.

__len__(

)

returns the number of distinct sequences that are aligned to the source interval. Note: this is NOT necessarily equal to the number of items that will be returned by the above iterators, since a single target sequence might have multiple 1:1 intervals of alignment to the source interval, due to indels.

In addition to these standard dictionary methods, NLMSASlice provides several additional methods and attributes:

• letters: this attribute provides an interface to the individual alignment columns (NLMSANode objects) containing the source interval, in order from start to stop. This provides an easy way to obtain detailed information about the letter-to-letter alignment of different sequences within this region of the alignment. For details on the kinds of information you can obtain for each alignment column, see NLMSANode, below.

  It also provides a graph interface to subset of the partial order alignment graph corresponding to this slice. For details, see NLMSASliceLetters, below.

split(

**kwargs)

this method provides a way to perform group-by operations on the slice; the output of split() is one or more NLMSASlice objects; if the group-by analysis results in no splitting of the current slice, then it is returned unchanged (i.e. the method just returns self). Uses same group-by arguments as keys(). For further details on group-by operations, see keys() above.

regions(

**kwags)

performs the same group-by analysis as split(), but replaces the source interval by the corresponding interval in the LPO. The main practical consequence of this is that target sequence inserts are included in the resulting slice (because they are present in the LPO interval corresponding to the original source interval), whereas they were NOT included in the original slice (because they are not aligned to the source interval). The main place where this matters is in graph traversal of the slice's letters attribute: whereas the nodes and edges corresponding to these inserts are not considered to be part of the letters graph for the original slice, they are part of the LPO slice. Also, the ``source interval'' in any subsequent operations with the LPO slice will be LPO coordinates instead of subintervals of the original source sequence interval. Uses same group-by arguments as keys().

groupByIntervals(

maxgap=0, maxinsert=0, mininsert= 0, filterSeqs=None, mergeAll=True, ivalMethod=None, **kwargs)

This method performs the interval grouping analysis for all the iterators described above. Users will not need to call it directly. Its arguments are described above (see keys()). It returns a dictionary whose keys are sequences aligned to this slice, and whose values are the list of intervals produced by the group-by analysis for the corresponding sequence. The values are tuples of the form (source_start, source_stop, target_start, target_stop), showing the mapping of a source sequence interval onto a target sequence interval. This dictionary is the primary input to the groupBySequences() method below.

groupBySequences(

seqIntervals, sourceOnly=False, indelCut=False, seqGroups=None, minAligned=1, pMinAligned=0., seqMethod=None, **kwargs)

This method performs the sequence grouping analysis for all the iterators described above. seqIntervals must be a dictionary of sequences and their associated list of intervals (produced by groupByIntervals() above). It returns a list of output sequence intevals, which is either a list of source sequence intervals (sourceOnly mode), or a list of tuples of the form (source_interval, target_interval).

matchIntervals(

seq=None)

this method returns the set of 1:1 match intervals for the target sequence seq (or all aligned sequences, if seq is None), as a dictionary whose keys are target sequence intervals, and whose values are the corresponding source sequence intervals to which they are aligned.

findSeqEnds(

seq)

returns the largest possible interval of seq that is aligned to this slice, i.e. it merges all alignment intervals in this slice containing seq, and returns the merged sequence interval based on the minimum start value and maximum stop value found.

2.4.6 NLMSASliceLetters

This graph has the following methods:

__iter__(

)

generates all the nodes in the slice, in order from left to right.

items(

)

also iteritems(). Generate the same set of nodes as above, as keys, but for each also returns a value representing its outgoing directed edges (see getitem, below).

__getitem__(

node)

gets a dictionary indicating all the outgoing directed edges from node to subsequence nodes, whose keys are the target nodes, and whose edges are the sequence.LetterEdge objects representing each edge.

2.4.7 NLMSANode

__iter__(

)

generates all the individual sequence letters (as SeqPath intervals, presumably of length 1) that are aligned to the source sequence, in this column of the alignment.

edges(

)

generates the same list of of target sequence letters as the iterator, but as a tuple of (target letter, source letter, edge). Currently, edge is just None.

__len__(

)

returns the number of distinct sequences aligned to the source interval, in this column.

Other, internal methods that regular users are unlikely to need:

getSeqPos(

seq)

returns the sequence interval of seq that is aligned to this column, or raises KeyError if it is not aligned here.

getEdgeSeqs(

node2)

returns a dictionary of sequences that traverse the edge directly from this node to node2, i.e. if letter i of seq is aligned to this node, then letter letter i+1 is aligned to node2. The dictionary's keys are top-level sequence objects, and its value for each is the letter position index i as defined above.

nodeEdges(

)

returns a dictionary of the outgoing edges from this node, whose keys are target nodes, and whose values are the corresponding edge objects (of type sequence.LetterEdge).

2.4.8 NLMSASequence

However, one thing you should know is that for a sequence to be stored in a NLMSA, it needs to have a unique string identifier. NLMSASequence obtains a string identifier for the sequence in one of the following ways (in decreasing order of precedence): 1) the sequence ``object'' can itself just be a Python string, in which case that string is used as the identifier. 2) otherwise, the object should be a SeqPath instance. If it has a name attribute, that will be used as the identifier. 3) Otherwise, if it has a id attribute (which is present by default on sequence.Sequence objects), that will be used.

2.5 Seqdb Module

Pygr interface to sequence databases stored in FASTA, BLAST or relational databases.

2.5.1 Overview

The seqdb module provides a simple, consistent interface to sequence databases from a variety of different storage sources such as FASTA, BLAST and relational databases. Sequence databases are modeled (like other Pygr container classes) as dictionaries, whose keys are sequence IDs and whose values are sequence objects. Pygr sequence objects use the Python sequence protocol in all the ways you'd expect: a subinterval of a sequence object is just a Python slice (s[0:10]), which just returns a sequence object representing that interval; the reverse complement is just -s; the length of a sequence is just len(s); to obtain the actual string sequence of a sequence object is just str(s). Pygr sequence objects work intelligently with different types of back-end storage (e.g. relational databases or BLAST databases) to efficiently access just the parts of sequence that are requested, only when an actual sequence string is needed.

2.5.2 External Requirements

• NCBI toolkit: The BLAST database functionality in this module requires that the NCBI toolkit be installed on your system. Specifically, some functions will call the command line programs formatdb, fastacmd, blastall, and megablast.

• RepeatMasker: the BlastDB.megablast() method calls the command line program RepeatMasker to mask out repetitive sequences from seeding alignments, but to allow extension of alignments into masked regions.

• Python DB-API 2.0: the SQLTable class, and dependent classes such as SQLSequence and StoredPathMapping, conform to the Python DB-API 2.0. Typically you must supply a DB-API 2.0-compliant database cursor to the SQLTable constructor. To do so, you must have some DB-API 2.0-compliant module (such as MySQLdb) installed for connecting to a database server.

If you are lacking one or more of these requirements, you can still install Pygr and use all Pygr functionality that does not depend on the missing requirements. If you try to use a function for which a requirement is missing, Pygr will raise an appropriate exception (e.g. unable to run blastall).

2.5.3 BlastDB

from pygr.sequence import *
db=BlastDB('sp') # OPEN SWISSPROT BLAST DB
s=Sequence(str(db['CYGB_HUMAN'][40:-40]),'boo')
m=db.blast(s) # DO BLAST SEARCH
myg=db['MYG_CHICK']
for i in m[s][myg]:
    print repr(i.srcPath),repr(i.destPath),i.blast_score,i.percent_id

Let's go through this example line by line:

• construction of a BlastDB object: This looks for either a FASTA file with the path 'sp' or BLAST database formatted files based on this path (e.g. 'sp.psd' for protein sequences, or 'sp.nsd' for nucleotide sequences).

• db['CYGB_HUMAN'] obtains a sequence object representing the SwissProt sequence whose ID is CYGB_HUMAN. The slice operation [40:-40] behaves just like normal Python slicing: it obtains a sequence object representing the subinterval omitting the first 40 letters and last 40 letters of the sequence. The str() operation obtains the actual string representation of this subinterval.

• Sequence(letter_string, name) creates a new sequence object whose sequence is letter_string, and whose ID is name.

• Running the db.blast(s) method searches the BLAST database for homologies to s, using NCBI BLAST. It chooses reasonable parameters based upon the sequence types of the database and supplied query. However, you can specify extra parameter options if you wish. It returns a Pygr sequence mapping (multiple alignment) that represents a standard Pygr graph of alignment relationships between s and the homologies that were found.

• The expression m[s][myg] obtains the "edge information" for the graph relationship between the two sequence nodes s and myg. (if there was no edge in the m graph representing a relationship between these two sequences, this would produce a KeyError). This edge information consists of a set of interval alignment relationships (described in detail below), which are printed out in this example.

Options for constructing a BlastDB:

BlastDB(

filepath=None,skipSeqLenDict=False,ifile=None,idFilter=None, blastReady=False)

Open a sequence file as a ``database'' object, giving the user access to its sequences, easy searching via blast or megablast, etc. skipSeqLenDict prevents construction of a sequence length index file (which will be named ``filepath.seqlen'') and a fast sequence access file (which will be named ``filepath.pureseq''). Setting this option to True can be useful if you either wish to speed up initial opening of the BlastDB (note: construction of these indexes is only a one-time event) or avoid the extra disk space required by these indexes. Explanation: To facilitate the rapid creation of sequence objects (which requires the length of the sequence), it creates a sequence length index (as a Python shelve). This enables it to avoid actually loading the sequence string into memory each time a sequence object is created; instead it just looks up the sequence length. While this speeds up access to genomic sequence databases (where each sequence tends to be extremely long), this initial step may be slow for databases of short sequences. Setting skipSeqLenDict option to True, will prevent construction of this sequence length index. ifile lets you open the database directly from a file object rather than a filename. If you have a file object, you can pass it directly to BlastDB instead of a filepath. NB: the BlastDB() constructor will close ifile when it is done reading from the file object. idFilter allows you to provide a function for re-mapping the FASTA sequence identifiers read from the sequence file. This can be useful in the case of NCBI FASTA files, since NCBI often treats the sequence ID as a ``blob'' into which any number of database fields can be stuffed, rather than a true ID. blastReady option specifies whether BLAST index files should be automatically constructed (using formatdb). Note, if you attempt to run the blast() method, it will automatically create the index files for you if they are missing.

Useful methods:

iter(

)

iterate over all IDs in the BLAST database.

len(

)

returns number of sequences in the BLAST database.

blast(

seq,al=None,blastpath='blastall',blastprog=None,expmax=0.001,maxseq=None)

run a BLAST search on sequence object seq. Maxseq will limit the number of returned hits to the best maxseq hits.

megablast(

seq,al=None,blastpath='megablast',expmax=1e-20,maxseq=None,minIdentity =None,maskOpts='-U T -F m',rmOpts='')

first performs repeat masking on the sequence by converting repeats to lowercase, then runs megablast with command line options to prevent seeding new alignments within repeats, but allowing extension of alignments into repeats. minIdentity should be a number (maximum value, 100) indicating the minimum percent identity for hits to be returned.

__invert__(

)

The invert operator (~, the ``tilde'' character) enables reverse-mapping of sequence objects to their string ID.

id=(~db)[seq] # GET IDENTIFIER FOR THIS SEQUENCE FROM ITS DATABASE

Useful attributes:

• seqClass: the object class to use for instantiating new sequence objects from this BLAST database. You can set this to create customized sequence behaviors.

2.5.4 FileDBSequence

2.5.5 BlastSequence

This was previously the default class for sequence objects returned from BlastDB, but has been deprecated because we found that NCBI fastacmd was much too slow and consumed enormous amounts of memory. BlastSequence relies on fastacmd for ``fast'' access to individual sequence slices. The advantage is that it only requires BLAST database files (produced by Pygr using formatdb), whereas the new FileDBSequence requires a specially indexed sequence file (constructed by default by BlastDB), which may be a disadvantage if you are low on disk space.

BlastSequence has several optimizations for working with BLAST databases:

• it uses the NCBI tool fastacmd to retrieve sequence efficiently from a BLAST database, when your program requests an actual string of sequence text. Moreover, for subintervals (slices) of the sequence, it uses fastacmd's -L option to request just the desired subinterval of the sequence, rather than the whole sequence. This makes it efficient for requesting specific intervals of large genomic contigs. Basically, just use Python slicing and str() methods on sequence objects, and subsequences will be obtained in an efficient manner.

• the len() method is implemented using the seqLenDict, a precalculated index of the sequence lengths. So again no sequence has to be read by Python.

2.5.6 BlastSequenceCache

Implements a variant of BlastSequence designed to merge and cache requests for local intervals of sequence so that repeated accesses to these regions are bundled and cached for efficiency. You work with sequence objects of this type normally, using Python slicing to obtain subintervals, and str() to get the sequence string for a subinterval. But behind the scenes, it does two things:

• all slicing operations are recorded, in the form of a cache of superintervals. Overlapping or adjacent intervals are merged into a superinterval up to a maximum superinterval size (default 20000). It will automatically create as many superintervals as needed to cover the requested subinterval slices. Each superinterval is represented by an object of the FastacmdIntervalCache class.

• when the sequence string of a subinterval is requested, the cache actually retrieves (and caches) the entire superinterval containing that subinterval. Fastacmd only needs to be called once for this superinterval. Subsequent subinterval string requests that fall within this cached superinterval are simply returned directly from the cache, without calling fastacmd.

2.5.7 SQLSequence

Implements a subclass inheriting from SQLRow and SequenceBase, to use a relational database table to obtain the actual sequence. There are three minor variants DNASQLSequence, RNASQLSequence, ProteinSQLSequence (so that the sequence does not have to analyze itself to determine what kind of sequence it is). Its constructor takes the same arguments as SQLRow(table, id), where table is the SQLTable object representing the table in which the sequence is stored, and id is the primary key of the row representing this sequence. However, normally this class is simply passed to the Table object itself so that it will use it to instantiate new row objects whenever they are requested via its dictionary interface. Here's a simple example:

class YiProteinSequence(ProteinSQLSequence): # CREATE A NEW SQL SEQUENCE CLASS
    def __len__(self): return self.protein_length  # USE LENGTH STORED IN DATABASE
protein=jun03[protein_seq_t] # protein IS OUR SQLTable OBJECT REPRESENTING PROTEIN SEQUENCE TABLE
protein.objclass(YiProteinSequence) # FORCE PROTEIN SEQ TABLE TO USE THIS TO INSTANTIATE ROW OBJECTS
pseq=protein['Hs.1162'] # GET PROTEIN SEQUENCE OBJECT FOR A SPECIFIC CLUSTER

Let's go through this line by line:

• we create a subclass of ProteinSQLSequence to show how Python makes it easy to create customized behaviors that can make database access more efficient. Here we've simply added a __len__ method that uses the protein_length attribute obtained directly from the database, courtesy of SQLRow.__getattr__, which knows what columns exist in the database, and provides them transparently as object attributes. (The ordinary SequenceBase __len__ method calculates it by obtaining the whole sequence string and calculating its length. Clearly it's more efficient for the database to retrieve this number (stored as a column called protein_length) and return it, rather than making it send us the whole sequence).

• next we call the protein.objclass() method to inform the table object that it should use our new class for instantiating any row objects for this table. It will call this class with the usual SQLRow contructor arguments (table, id).

2.5.8 StoredPathMapping

A second major area in Pygr is representation and query of multiple sequence alignment databases in a way that is scalable to whole genomes. We have previously showed (in our work on Partial Order Alignment) that graphs provide both a compact and algorithmically powerful way to store alignments. Combining this with "interval alignment" makes it scalable and gives a simple interface. In Pygr, alignments are just another kind of graph, whose nodes are sequence intervals, and edges are alignment relations. This provides a general-purpose facility for working with sets of sequence intervals, sequence annotation databases, and multiple sequence alignments, all queryable via Pygr graph queries. We have implemented different container subclasses to work with these data in memory or to work transparently with data stored in relational databases. The consistency and simplicity of the Pygr framework makes it a good interface both to run external tools like BLAST, and to store or query the results in persistent storage like a MySQL database.

hg17=BlastDB('/data/ucsc/hg17') # GET CONTAINER FOR HUMAN GENOME DATABASE
bcl2m=hg17['chr22'][16544303:16588541] # GET INTERVAL WITH BCL2L13 GENE
al=hg17.megablast(mouse_bcl2,maxseq=1) # GET REPEAT-MASKED MEGABLAST ALIGNMENT, ONLY TOP HIT
al[bcl2m[1000:1100] ]+=mrna[210:310] #ADD ALIGNMENT OF A 100nt SEGMENT TO mrna SEGMENT
al.storeSQL('test.table',db_cursor) # STORE COMPLETE ALIGNMENT IN RELATIONAL DATABASE
for e in MAFStoredPathMapping(bcl2m,'ucsc_maf8',u).edges(): #GET ITS MULTIGENOMEALIGNMENTS
   print str(e.srcPath),str(e.destPath) # PRINT THE ACTUAL ALIGNED SEQUENCE INTERVALS

Note: We will be unifying all sequence alignment functionality under the NLMSA interface design sometime in the near future. Specifically, the PathMapping and related classes, while similar to NLMSA, will be replaced with interfaces that are identical to the NLMSA interface.

2.5.9 SliceDB

• slice database: must accept a string identifier as a key, and return a slice information object as a value.

• slice information: a slice information object must have the following attributes: name gives the identifier of the sequence containing the slice; start and stop give the coordinates of the sequence interval (which should be positive integers following standard Python slice coordinate conventions); ori gives the sequence orientation as an integer (1 for positive orientation, -1 for negative orientation).

• sequence database: must accept a string identifier as a key, and return a sliceable sequence object as a value.

Both databases should raise KeyError for bad identifiers. The current SliceDB implementation caches sequence objects so that subsequent calls for the same identifier will not require repeating the database queries to the two databases. To remove a sequence object from the cache, just use del db[id] as usual.

SliceDB inherits from the builtin Python dict class, so all standard methods can be used.

db=SliceDB(sliceDB,seqDB) # CREATE OUR DATABASE
gene=db[cluster_id] # USE IT TO GET A GENE SEQUENCE...

2.5.10 VirtualSeq

s=VirtualSeq('FOOG_HUMAN')
len(s) # ONLY 1 LETTER LONG BY DEFAULT
s1=s[100:215] # GET A SLICE OF THIS SEQUENCE
len(s) # NOW IT'S 215

The associated VirtualSeqDB class provides a ``sequence database'' that returns a VirtualSeq object for every identifier requested of it. It acts like a Python dictionary:

db=VirtualSeqDB()
s=db['FOOG_HUMAN'] # ASK FOR A SEQUENCE BY ITS IDENTIFIER
s1=s[100:215] # GET A SLICE OF THIS SEQUENCE

2.5.11 PrefixUnionDict

hg17=BlastDB('/usr/tmp/ucsc_msa/hg17')
mm5=BlastDB('/usr/tmp/ucsc_msa/mm5')
... # LOAD A BUNCH OF OTHER GENOMES TOO...
genomes={'hg17':hg17,'mm5':mm5, 'rn3':rn3, 'canFam1':cf1, 'danRer1':dr1,
'fr1':fr1, 'galGal2':gg2, 'panTro1':pt1} # PREFIX DICTIONARY FOR THE UNION 
					 # OF ALL OUR GENOMES
genomeUnion=PrefixUnionDict(genomes)
ptChr7=genomeUnion['panTro1.chr7'] # GET CHIMP CHROMOSOME 7

if 'panTro1.chr5' in genomeUnion: # CHECK IF THIS ID IN OUR UNION
    pass # DO SOMETHING...

s= -(ptChr7[1000:2000]) # GET A BIT OF THIS SEQUENCE
if s in genomeUnion: # THIS IS HOW TO CHECK IF s DERIVED FROM OUR UNION
    pass # DO SOMETHING...

It provides a __contains__ method that tests whether a given sequence object is derived from the PrefixUnionDict (see example above). Here are some additional methods:

__init__(

prefixDict=None,separator='.',filename=None,dbClass=BlastDB)

You can create a PrefixUnionDict either using a prefixDict (whose keys are string prefixes, and whose values are sequence databases), or using a previously created header file filename. Using the header file, the constructor will automatically open all the sequence databases for you. When opening from a header file, you can also specify a dbClass to be used for opening individual sequence databases listed in the header file; the default is BlastDB. The database class constructor must take a single argument, which is the filepath for opening the database. The separator character is used to form ``prefix.suffix'' identifiers.

__invert__(

)

The invert operator (~, the ``tilde'' character) enables reverse-mapping of sequence objects to their string ID. This is the recommended way to get the ``fully qualified sequence ID'', i.e. with the appropriate prefix prepended.

id=(~db)[seq] # GET PROPERLY PREFIXED-IDENTIFIER FOR THIS SEQUENCE

For a given sequence object seq derived from the union (or a slice of a sequence from the union), return a string identifier in the form of ``foo.bar''.

getName(

path)

This method is deprecated; instead use the __invert__ operator above.

writeHeaderFile(

filename)

Save a header file for this union, to reopen later. It saves the separator character, and a list of prefixes and filepaths to the various sequence databases (which must have a filepath attribute). This header file can be used for later reopening the prefix-union in a single step.

2.5.12 Functions

read_fasta(

ifile, onlyReadOneLine=False)

a generator function that yields tuples of id,title,seq from ifile. The onlyReadOneLine option is useful if you only want to determine the sequence type (e.g. DNA vs protein) and the whole sequence might be extremely long (e.g. a genome!).

write_fasta(

ofile, s, chunk=60, id=None)

writes the sequence s to the output file ofile, using chunk as the line width. id can provide an identifier to use instead of the default s.id.

2.6 sqlgraph Module

2.6.1 SQLTable

This class is derived from the Python builtin dict class, so all standard methods can be used.

objclass(

oclass)

Specify a object class to use for creating new ``row'' objects. oclass must accept a single argument, a tuple object representing a row in the database.

Otherwise, the default oclass for SQLTable is the TupleO class, which provides a named attribute interface to the tuple values representing the row.

load(

oclass=None)

Load all data from the table, using oclass as the row object class if specified (otherwise use the oclass for this table). All rows are loaded from the database and saved as row objects in the Python dictionary of this class.

select(

whereClause,params=None,oclass=None)

Generate the list of objects that satisfy the whereClause via a SQL SELECT query. This function is a generator, so you use it as an iterator. params is passed to the cursor execute statement to allow additional control over the query.

There are several variants of this class:

2.6.2 SQLTableNoCache

2.6.3 SQLMultiTableNoCache

2.6.4 TupleO

2.6.5 SQLRow

2.7 Coordinator Module

Framework for running subtasks distributed over many computers, in a pythonic way, using SSH for secure process invocation and XMLRPC for message passing. Also provides simple interface for queuing and managing any number of such "batch jobs".

2.7.1 Overview

The coordinator module provides a simple system for running a large collection of tasks on a set of cluster nodes. It assumes:

• authentication is handled using ssh-agent. The coordinator module does no authentication itself; it simply tries to spawn jobs to remote nodes using ssh, assuming that you have previously authenticated yourself to ssh-agent.

• the client nodes can access your scripts using the same path as on the initiating system. In other words, if you launch a coordinator job /home/bob/mydir/myscript.py, your client nodes must also be able to access /home/bob/mydir/myscript.py (e.g. via NFS).

• your job consists of a large set of task IDs, and a computation to be performed on each ID. To run this job, you provide an iterator that generates the list of task IDs for the Coordinator to distribute to your client nodes. You start your script to run a Coordinator that serves your list of task IDs to the client nodes. You also provide a function that performs your desired computation on each task ID it receives from the Coordinator. Typically, you provide both the server function (i.e. the iterator that generates the list of task IDs) and the client function (that runs your desired computation for each ID) within a single Python script file. Running this script without extra flags starts the Coordinator, which in turn launches your script as a Processor on one or more client nodes. The Processors andCoordinator work together to complete all the task IDs.

• a ResourceController performs load balancing and resource allocation functions, including: dividing up loads from one or more Coordinators over a set of hosts (each with one or more CPUs); serving a Resource database to Processors requesting specific resources; resource-locking on a per node basis for preventing Processors from using a Resource that is under construction by another Processor. For very large files that are used repeatedly by your computation, it is preferable to first copy them to local disk on each cluster node (fast), rather than reading them over and over again from NFS (slow). Resources provide a simple mechanism for doing this.

from pygr.apps.leelabdb import *
from pygr import coordinator

def map_clusters(server,genome_rsrc='hg17',dbname='HUMAN_SPLICE_03',
                 result_table='GENOME_ALIGNMENT.hg17_cluster_JUN03_all',
                 rmOpts='',**kwargs):
    "map clusters one by one"
    # CONSTRUCT RESOURCE FOR US IF NEEDED
    genome=BlastDB(ifile=server.open_resource(genome_rsrc,'r'))
    # LOAD DB SCHEMA
    (clusters,exons,splices,genomic_seq,spliceGraph,alt5Graph,alt3Graph,mrna, \
     protein,clusterExons,clusterSplices)=getSpliceGraphFromDB(spliceCalcs[dbname])

    for cluster_id in server:
        g=genomic_seq[cluster_id]
        m=genome.megablast(g,maxseq=1,minIdentity=98,rmOpts=rmOpts) # MASK, BLAST, READ INTO m
        # SAVE ALIGNMENT m TO DATABASE TABLE test.mytable USING cursor
        createTableFromRepr(m.repr_dict(),result_table,clusters.cursor,
                            {'src_id':'varchar(12)','dest_id':'varchar(12)'})
        yield cluster_id # WE MUST FUNCTION AS GENERATOR

def serve_clusters(dbname='HUMAN_SPLICE_03',
                   source_table='HUMAN_SPLICE_03.genomic_cluster_JUN03',**kwargs):
    "serve up cluster_id one by one"
    cursor=getUserCursor(dbname)
    t=SQLTable(source_table,cursor)
    for id in t:
        yield id

if __name__=='__main__':
    coordinator.start_client_or_server(map_clusters,serve_clusters,['hg17'],__file__)

Let's analyze the script line by line:

• mapclusters() is a client generator function to be run in a Processor on a client node. It takes one argument representing its connection to the server (a Processor object), and optional keyword arguments read from the command line. It first does some initial setup (opens a BLAST database and loads a schema from a MySQL database), then iterates over task IDs returned to it from the server. A few key points:

• server.open_resource(genome_rsrc,'r') requests a resource given by the genome_rsrc argument from the ResourceController, does whatever is necessary to copy this resource to local disk, and then opens it for reading, returning a file-like object. This can then be used however you like, but you MUST call its close() method (just as you should always do for any file object) to indicate that you're done using it. Failure to close() the file object will leave the Resource "hg17" permanently locked on this specific node. (You would then have to unlock it by hand using the ResourceController.release_rule() method).

• yield cluster_id: the client function must be a Python generator function (i.e. it must use the yield statement), and it must yield the list of IDs that it has processed. Python's generator construct is extremely convenient for many purposes: here it lets us perform both our initializations and iteration over IDs within a single function, while at the same time wrapping each iteration within the Processor's error trapping code (to prevent a single error in your code from causing the entire Processor to shut down). The Processor will trap any errors in your code and and send tracebacks to your Coordinator, which will report them in its logfile. The Processor will tolerate occasional errors and continue processing more IDs. However, if more than a certain number of IDs in a row fail with errors (controlled by the Processor.max_errors_in_a_row attribute), the Processor will exit, on the assumption that either your code or this specific client node don't work correctly.

• serve_clusters() is the server generating function to be run in the Coordinator. It returns an iterator that generates all the task IDs that we want to run. Again, the Python generator construct provides a very clean way of doing this: we simply yield each ID that we want to process in our client Processors.

• if __name__=="__main__": this final clause automatically launches our script as either a Coordinator or Processor depending on the command line options (which are automatically parsed by start_client_or_server()). All we have to do is pass the client generator function, the server generator function, a list of the resources this job will use, and the name of the script file to be run on client nodes. Since that is just this script itself, we use the Python builtin symbol __file__ (which just evaluates to the name of the current script).

• Command-line arguments are parsed (GNU-style, ie. -foo=bar) by start_client_or_server() and passed to your client and server functions as Python named parameters. Because the same list of arguments is passed to your client and server functions, and each of these functions won't necessarily want to get all the named arguments, you should include the **kwargs at the end of the argument list. Any unmatched arguments will be stored in kwargs as a Python mapping (dictionary). If you fail to do this, your client or server function will crash if called with any named parameters other than the ones it expects.

2.7.2 Log and Error Information

Process logging and error information go to three different types of logs:

• Processor logfile(s): every individual Processor (and all subprocesses run by it) send stdout and stderr to a logfile on local disk of the host on which it is running. Currently the filename is /usr/tmp/NAME_N.log, where NAME is the name you assigned to the job when you started the Coordinator, and N is the numeric ID of the Processor assigned by the coordinator (just an auto-increment integer beginning at 0, and increasing by one for each Processor the Coordinator starts). This logfile is the place to look if your job is failing mysteriously-look in the logfile and see its last words before its demise. You can get a complete list of the logfiles for all the Coordinator's Processors by inspecting the logfile attribute of the CoordinatorMonitor (see below).

• Coordinator logfile: all XMLRPC requests from client Processors, as well as error messages from them, are logged here. All Python errors (tracebacks) in your client (Processor) code are reported here. Also, the actual SSH commands used to invoke your Processors on cluster nodes, are logged here. This is usually the place to start, to see whether things are going well (you should see a long stream of next requests as Processors finish a task and request the next one), or failing with errors.

• ResourceController logfile: all XMLRPC requests from Processors and Coordinatorsare logged here, including register() and unregister(), resource requests, and load reporting from cluster nodes. If things are working well, you should see a stream of regular report_load() messages showing steady, full utilization of all the host processors. Excessive register/unregister churning (jobs that start and immediately exit) is a common sign of trouble with your jobs.

2.7.3 Coordinator

To start a job coordinator (which in turn will the run the whole job by starting Processors on cluster nodes using SSH):

python mapclusters5.py mm5_jan02 --errlog=/usr/tmp/leec/mm5_jan02.log \ 
  --dbname=MOUSE_SPLICE --source_table=genomic_cluster_jan02 \
  --genome_rsrc=mm5 --result_table=GENOME_ALIGNMENT.mm5_cluster_jan02_all \ 
  --rmOpts=-rodent \

Here we have told the Coordinator to name itself "mm5_jan02" in all its communications with the ResourceController. Since we gave no command-line flags, the Coordinator will assume that a ResourceController is already running on port 5000 of the current host. You must have an ssh-agent running BEFORE you start the Coordinator, since the Coordinator will attempt to spawn jobs using SSH. The Coordinator will exit with an error message if it is unable to connect to ssh-agent. A few notes:

• The Coordinator will run as a demon process (i.e. in the background, and detached from your terminal session), and redirect its output into a file (here, given by the -errlog option). If you don't specify an -errlog filename, it will create a filename determined by the name we told it to run as, in this case "mm_jan02.log".

• You must ensure that SSH can launch processes on your client nodes "unattended" i.e. without a connection to a controlling terminal. If SSH has to ask for userconfirmations when connecting to a given host (e.g. if it asks whether you want to accept the host key), the Coordinator will not be able to use that host.

• Python errors (tracebacks) in your will be GNU-style command-line options (e.g. -port=8889) are automatically parsed by start_client_or_server() and passed to the Coordinator.__init__() as keyword arguments. This constructor takes the following optional arguments:
  • port: the port number on which this Coordinator should run

  • priority: a floating point number specifying the priority level at which this Coordinator should be run by the ResourceController. The default value is 1.0. A value of 2.0 will give it twice as many Processors as a competing Coordinator of priority 1.0.

  • rc_url: the URL for the ResourceController. Defaults to http://THISHOST:5000
  • errlog: logfile path for saving all output to. Defaults to NAME.log, where NAME is the name you assigned to this Coordinator. Can be an absolute path.

  • immediate: if True, make the job run immediately, without waiting for previous jobs to finish. Default: False.

  • demand_ncpu: if set to a non-zero value, specifies the exact number of Processors you want to run your job.

  • NB: command line arguments are also passed to your server function, and to your client function, as Python named parameters. See the mapclusters5.py example above.

2.7.4 ResourceController

Whereas you start a separate Coordinator for each set of jobs you want to run, you only need a single ResourceController running. To start the ResourceController, run:

python coordinator.py --rc=bi