Showing posts with label tree. Show all posts
Showing posts with label tree. Show all posts

13 July 2012

Parsing the Newick format in C using flex and bison.

The following post is my answer for this question on biostar "Newick 2 Json converter".
The Newick tree format is a simple format used to write out trees (using parentheses and commas) in a text file .
The original question asked for a parser based on perl but here, I've implemented a C parser using flex/bison.


Example:

((Human:0.3, Chimpanzee:0.2):0.1, Gorilla:0.3, (Mouse:0.6, Rat:0.5):0.2);

A formal grammar for the Newick format is available here
Items in { } may appear zero or more times.
   Items in [ ] are optional, they may appear once or not at all.
   All other punctuation marks (colon, semicolon, parentheses, comma and
         single quote) are required parts of the format.


              tree ==> descendant_list [ root_label ] [ : branch_length ] ;

   descendant_list ==> ( subtree { , subtree } )

           subtree ==> descendant_list [internal_node_label] [: branch_length]
                   ==> leaf_label [: branch_length]

            root_label ==> label
   internal_node_label ==> label
            leaf_label ==> label

                 label ==> unquoted_label
                       ==> quoted_label

        unquoted_label ==> string_of_printing_characters
          quoted_label ==> ' string_of_printing_characters '

         branch_length ==> signed_number
                       ==> unsigned_number

The Flex Lexer

The Flex Lexer is used to extract the terminal tokens of the grammar from the input stream.
Those terminals are '(' ')' ',' ';' ':' , strings and numbers. For the simple and double quoted strings, we tell the lexer to enter in a specific state ( 'apos' and 'quot').

The Bison Scanner

The Bison scanner reads the tokens returned by Flex and implements the grammar.
The simple structure holding the tree is defined in 'struct tree_t'. The code also contains some methods to dump the tree as JSON.

Makefile


Testing

compile:
$ make
bison -d newick.y
flex newick.l
gcc -Wall -O3 newick.tab.c lex.yy.c
lex.yy.c:1265:17: warning: ‘yyunput’ defined but not used [-Wunused-function]
lex.yy.c:1306:16: warning: ‘input’ defined but not used [-Wunused-function]
test:
echo "((Human:0.3, Chimpanzee:0.2):0.1, Gorilla:0.3, (Mouse:0.6, Rat:0.5):0.2);" | ./a.out

{
    "children": [
        {
            "length": 0.1,
            "children": [
                {
                    "label": "Human",
                    "length": 0.3
                },
                {
                    "label": "Chimpanzee",
                    "length": 0.2
                }
            ]
        },
        {
            "label": "Gorilla",
            "length": 0.3
        },
        {
            "length": 0.2,
            "children": [
                {
                    "label": "Mouse",
                    "length": 0.6
                },
                {
                    "label": "Rat",
                    "length": 0.5
                }
            ]
        }
    ]
}


That's it,

Pierre

07 January 2011

SVG image in wikipedia + Zoom.it =a Zoomable "Tree of life"

Zoom.it is a free service for viewing and sharing high-resolution imagery. You give us the link to any image ( including SVG, pdfs) on the web along with a nice short URL.. As a test, I used the SVG file "Tree of life with genome size.svg" on wikipedia and here is the awesome result generated by Zoom.it:



That's it,
Piere

21 September 2010

Trees in Mongodb, my notebook with Gene Ontology

In the current post I've loaded the Gene Ontology into MongoDB and played with the tree structure of the database:

Loading GeneOntology into MongoDB

First, download GO as RDF at http://archive.geneontology.org/latest-termdb/go_daily-termdb.rdf-xml.gz and transform it with my XSLT stylesheet go2mongo.xsl (available here):
<?xml version='1.0' encoding="UTF-8" ?>
<xsl:stylesheet
xmlns:xsl='http://www.w3.org/1999/XSL/Transform'
xmlns:go="http://www.geneontology.org/dtds/go.dtd#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
version='1.0'
>
<xsl:output method="text"/>

<xsl:param name="colName">go</xsl:param>

<xsl:template match="/">
<xsl:apply-templates/>
</xsl:template>

<xsl:template match="go:go">
<xsl:apply-templates select="rdf:RDF"/>
</xsl:template>

<xsl:template match="rdf:RDF">

db.<xsl:value-of select="$colName"/>.drop();

<xsl:apply-templates select="go:term"/>


</xsl:template>

<xsl:template match="go:term">
<xsl:text>term={_id:</xsl:text><xsl:apply-templates select="go:accession" mode="text"/>
<xsl:text>,name:</xsl:text><xsl:apply-templates select="go:name" mode="text"/>
<xsl:if test="go:synonym">
<xsl:text>,synonyms:[</xsl:text>
<xsl:for-each select="go:synonym">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="." mode="text"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:definition">
<xsl:text>,definition:</xsl:text>
<xsl:apply-templates select="go:definition" mode="text"/>
</xsl:if>

<xsl:if test="go:comment">
<xsl:text>,comments:[</xsl:text>
<xsl:for-each select="go:comment">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="." mode="text"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:part_of">
<xsl:text>,part_of:[</xsl:text>
<xsl:for-each select="go:part_of">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="@rdf:resource"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:is_a">
<xsl:text>,is_a:[</xsl:text>
<xsl:for-each select="go:is_a">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="@rdf:resource"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:negatively_regulates">
<xsl:text>,negatively_regulates:[</xsl:text>
<xsl:for-each select="go:negatively_regulates">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="@rdf:resource"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:positively_regulates">
<xsl:text>,positively_regulates:[</xsl:text>
<xsl:for-each select="go:positively_regulates">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="@rdf:resource"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:regulates">
<xsl:text>,regulates:[</xsl:text>
<xsl:for-each select="go:regulates">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="@rdf:resource"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:dbxref">
<xsl:text>,dbxrefs:[</xsl:text>
<xsl:for-each select="go:dbxref">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="."/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:association">
<xsl:text>,associations:[</xsl:text>
<xsl:for-each select="go:association">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:text>{evidences:[</xsl:text>
<xsl:for-each select="go:evidence">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="." mode="text"/>
</xsl:for-each>
<xsl:text>],gene_product:{name:</xsl:text>
<xsl:apply-templates select="go:gene_product/go:name" mode="text"/>
<xsl:text>,dbxref:</xsl:text>
<xsl:apply-templates select="go:gene_product/go:dbxref" />
<xsl:text>}}</xsl:text>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>

<xsl:if test="go:is_obsolete">
<xsl:text>,is_obsolete:[</xsl:text>
<xsl:for-each select="go:is_obsolete">
<xsl:if test="position()!=1"><xsl:text>,</xsl:text></xsl:if>
<xsl:apply-templates select="@rdf:resource"/>
</xsl:for-each>
<xsl:text>]</xsl:text>
</xsl:if>
<xsl:text>};
db.</xsl:text>
<xsl:value-of select="$colName"/>
<xsl:text>.save(term);
</xsl:text>
</xsl:template>

<xsl:template match="go:dbxref">
<xsl:text>{database_symbol:</xsl:text>
<xsl:apply-templates select="go:database_symbol" mode="text"/>
<xsl:text>,reference:</xsl:text>
<xsl:apply-templates select="go:reference" mode="text"/>
<xsl:text>}</xsl:text>
</xsl:template>

<xsl:template match="*" mode="text">
<xsl:text>&quot;</xsl:text>
<xsl:call-template name="escape">
<xsl:with-param name="s" select="."/>
</xsl:call-template>
<xsl:text>&quot;</xsl:text>
</xsl:template>

<xsl:template match="@rdf:resource">
<xsl:text>{&apos;$ref&apos;:&apos;</xsl:text>
<xsl:value-of select="$colName"/>
<xsl:text>&apos;,&apos;$id&apos;:&apos;</xsl:text>
<xsl:value-of select="substring-after(.,'#')"/>
<xsl:text>&apos;}</xsl:text>
</xsl:template>


<xsl:template name="escape">
<xsl:param name="s"/>
<xsl:choose>
<xsl:when test="contains($s,'&quot;')">
<xsl:value-of select="substring-before($s,'&quot;')"/>
<xsl:text>\&quot;</xsl:text>
<xsl:call-template name="escape">
<xsl:with-param name="s" select="substring-after($s,'&quot;')"/>
</xsl:call-template>
</xsl:when>
<xsl:otherwise>
<xsl:value-of select="$s"/>
</xsl:otherwise>
</xsl:choose>
</xsl:template>



</xsl:stylesheet>


unzip and transform 'go_daily-termdb.rdf-xml' with the stylesheet to generate the javascript:
xsltproc go2mongo.xsl go_daily-termdb.rdf-xml > input.js
The file input.js looks like this:
term={
_id:"GO:0000001",
name:"mitochondrion inheritance",
synonyms:["mitochondrial inheritance"],
definition:"The distribution of mitochondria, including the mitochondrial genome, into daughter cells after mitosis or meiosis, mediated by interactions between mitochondria and the cytoskeleton.",
is_a:[
{'$ref':'go','$id':'GO:0048308'},
{'$ref':'go','$id':'GO:0048311'}
]
};
db.go.save(term);
term={_id:"GO:0000002",name:"mitochondrial genome maintenance",definition:"The maintenance of the structure and integrity of the mitochondrial genome; includes replication and segregation of the mitochondrial chromosome.",is_a:[{'$ref':'go','$id':'GO:0007005'}],dbxrefs:[{database_symbol:"InterPro",reference:"IPR009446"},{database_symbol:"Pfam",reference:"PF06420"}]};
db.go.save(term);
term={_id:"GO:0000003",name:"reproduction",synonyms:["GO:0019952","GO:0050876","reproductive physiological process"],definition:"The production by an organism of new individuals that contain some portion of their genetic material inherited from that organism.",is_a:[{'$ref':'go','$id':'GO:0008150'}],dbxrefs:[{database_symbol:"Wikipedia",reference:"Reproduction"}]};
db.go.save(term);
term={_id:"GO:0000005",name:"ribosomal chaperone activity",definition:"OBSOLETE. Assists in the correct assembly of ribosomes or ribosomal subunits in vivo, but is not a component of the assembled ribosome when performing its normal biological function.",comments:["This term was made obsolete because it refers to a class of gene products and a biological process rather than a molecular function."],is_a:[{'$ref':'go','$id':'obsolete_molecular_function'}]};
db.go.save(term);
term={_id:"GO:0042254",name:"ribosome biogenesis",synonyms:["GO:0007046","ribosomal chaperone activity","ribosome biogenesis and assembly"],definition:"The process of the formation of the constituents of the ribosome subunits, their assembly, and their transport to the sites of protein synthesis.",is_a:[{'$ref':'go','$id':'GO:0022613'}],dbxrefs:[{database_symbol:"InterPro",reference:"IPR001790"},{database_symbol:"InterPro",reference:"IPR004037"},{database_symbol:"InterPro",reference:"IPR007023"},{database_symbol:"InterPro",reference:"IPR012948"},{database_symbol:"SP_KW",reference:"KW-0690"},{database_symbol:"HAMAP",reference:"MF_00554"},{database_symbol:"HAMAP",reference:"MF_00699"},{database_symbol:"HAMAP",reference:"MF_00803"},{database_symbol:"HAMAP",reference:"MF_01852"},{database_symbol:"Pfam",reference:"PF00466"},{database_symbol:"Pfam",reference:"PF04939"},{database_symbol:"Pfam",reference:"PF08142"},{database_symbol:"PROSITE",reference:"PS01082"},{database_symbol:"Wikipedia",reference:"Ribosome_biogenesis"},{database_symbol:"SMART",reference:"SM00785"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR00436"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR01575"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR02729"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR03594"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR03596"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR03597"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR03598"}]};
db.go.save(term);
term={_id:"GO:0044183",name:"protein binding involved in protein folding",synonyms:["chaperone activity"],definition:"Interacting selectively and non-covalently with any protein or protein complex (a complex of two or more proteins that may include other nonprotein molecules) that contributes to the process of protein folding.",is_a:[{'$ref':'go','$id':'GO:0005515'}]};
db.go.save(term);
term={_id:"GO:0051082",name:"unfolded protein binding",synonyms:["binding unfolded ER proteins","chaperone activity","fimbrium-specific chaperone activity","glycoprotein-specific chaperone activity","histone-specific chaperone activity","ribosomal chaperone activity","tubulin-specific chaperone activity"],definition:"Interacting selectively and non-covalently with an unfolded protein.",is_a:[{'$ref':'go','$id':'GO:0005515'}],dbxrefs:[{database_symbol:"InterPro",reference:"IPR000397"},{database_symbol:"InterPro",reference:"IPR001305"},{database_symbol:"InterPro",reference:"IPR001404"},{database_symbol:"InterPro",reference:"IPR002194"},{database_symbol:"InterPro",reference:"IPR002777"},{database_symbol:"InterPro",reference:"IPR002939"},{database_symbol:"InterPro",reference:"IPR003095"},{database_symbol:"InterPro",reference:"IPR003708"},{database_symbol:"InterPro",reference:"IPR004127"},{database_symbol:"InterPro",reference:"IPR004226"},{database_symbol:"InterPro",reference:"IPR004487"},{database_symbol:"InterPro",reference:"IPR004961"},{database_symbol:"InterPro",reference:"IPR008971"},{database_symbol:"InterPro",reference:"IPR009033"},{database_symbol:"InterPro",reference:"IPR009169"},{database_symbol:"InterPro",reference:"IPR010236"},{database_symbol:"InterPro",reference:"IPR011599"},{database_symbol:"InterPro",reference:"IPR012713"},{database_symbol:"InterPro",reference:"IPR012714"},{database_symbol:"InterPro",reference:"IPR012715"},{database_symbol:"InterPro",reference:"IPR012716"},{database_symbol:"InterPro",reference:"IPR012717"},{database_symbol:"InterPro",reference:"IPR012718"},{database_symbol:"InterPro",reference:"IPR012719"},{database_symbol:"InterPro",reference:"IPR012720"},{database_symbol:"InterPro",reference:"IPR012721"},{database_symbol:"InterPro",reference:"IPR012722"},{database_symbol:"InterPro",reference:"IPR012724"},{database_symbol:"InterPro",reference:"IPR012725"},{database_symbol:"InterPro",reference:"IPR016153"},{database_symbol:"InterPro",reference:"IPR016154"},{database_symbol:"InterPro",reference:"IPR019805"},{database_symbol:"HAMAP",reference:"MF_00117"},{database_symbol:"PROSITE",reference:"MF_00117"},{database_symbol:"HAMAP",reference:"MF_00175"},{database_symbol:"PROSITE",reference:"MF_00175"},{database_symbol:"HAMAP",reference:"MF_00307"},{database_symbol:"PROSITE",reference:"MF_00307"},{database_symbol:"HAMAP",reference:"MF_00308"},{database_symbol:"PROSITE",reference:"MF_00308"},{database_symbol:"PROSITE",reference:"MF_00332"},{database_symbol:"HAMAP",reference:"MF_00505"},{database_symbol:"PROSITE",reference:"MF_00505"},{database_symbol:"HAMAP",reference:"MF_00600"},{database_symbol:"PROSITE",reference:"MF_00679"},{database_symbol:"HAMAP",reference:"MF_00790"},{database_symbol:"PROSITE",reference:"MF_00821"},{database_symbol:"HAMAP",reference:"MF_00822"},{database_symbol:"HAMAP",reference:"MF_01046"},{database_symbol:"HAMAP",reference:"MF_01152"},{database_symbol:"PROSITE",reference:"MF_01152"},{database_symbol:"HAMAP",reference:"MF_01183"},{database_symbol:"ProDom",reference:"PD010430"},{database_symbol:"Pfam",reference:"PF00684"},{database_symbol:"Pfam",reference:"PF01430"},{database_symbol:"Pfam",reference:"PF01556"},{database_symbol:"Pfam",reference:"PF01920"},{database_symbol:"Pfam",reference:"PF02556"},{database_symbol:"Pfam",reference:"PF02970"},{database_symbol:"Pfam",reference:"PF02996"},{database_symbol:"Pfam",reference:"PF03280"},{database_symbol:"PIRSF",reference:"PIRSF002356"},{database_symbol:"PIRSF",reference:"PIRSF002583"},{database_symbol:"PIRSF",reference:"PIRSF005261"},{database_symbol:"PRINTS",reference:"PR00625"},{database_symbol:"PRINTS",reference:"PR01594"},{database_symbol:"PROSITE",reference:"PS00298"},{database_symbol:"PROSITE",reference:"PS00750"},{database_symbol:"PROSITE",reference:"PS00751"},{database_symbol:"PROSITE",reference:"PS00995"},{database_symbol:"PROSITE",reference:"PS51188"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR00074"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR00115"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR00382"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR00809"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR02350"},{database_symbol:"JCVI_TIGRFAMS",reference:"TIGR03142"}]};
db.go.save(term);
term={_id:"GO:0000006",name:"high affinity zinc uptake transmembrane transporter activity",definition:"Catalysis of the transfer of a solute or solutes from one side of a membrane to the other according to the reaction: Zn2+(out) = Zn2+(in), probably powered by proton motive force. In high affinity transport the transporter is able to bind the solute even if it is only present at very low concentrations.",is_a:[{'$ref':'go','$id':'GO:0005385'}]};
db.go.save(term);
(...)
Here the notation {'$ref':'go','$id':'GO:0048308'} is a special object interpreted by mongo as a "Database Reference", a kind of forein-key/link/pointer to another document with a special method named 'fetch' retrieving the linked document.

Load 'input.js' into mongo
mongo mygodatabase input.js

Playing with the GeneOntology Tree

I'm going to look if a go-term is a descendant of one another. First, let's define two useful javascript recursive functions looking for the parent(s) of a given node threw the property is_a.
var goNodeIsA= function (childNode, parentId) {
if (childNode == null) {
return false;
}
if (childNode._id == parentId) {
return true;
}
if (!childNode.is_a) {
return false;
}
for (var i = 0; i < childNode.is_a.length; ++i) {
if (goNodeIsA(childNode.is_a[i].fetch(), parentId))
{
return true;
}
}
return false;
}

var goIsA=function (childId, parentId)
{
return goNodeIsA(db.go.findOne({_id:childId}), parentId);
}

Now, let's find if GO:0003723 (RNA binding) is a descendant of GO:0005488 (binding) ?
> goIsA("GO:0003723","GO:0005488");
true

And is GO:0003723 (RNA binding) is a descendant of GO:0050355 (triphosphatase activity) ?
> goIsA("GO:0003723","GO:0050355");
false

Loop over all the GO terms and find the descendants of GO:0050355 (triphosphatase activity):
> db.go.find({},{name:1,is_a:1}).forEach(function(term) { if(goIsA(term._id,'GO:0005488')) printjson(term); })

(...)
{
"_id" : "GO:0080084",
"name" : "5S rDNA binding",
"is_a" : [
{
"$ref" : "go",
"$id" : "GO:0000182"
}
]
}
{
"_id" : "GO:0080087",
"name" : "callose binding",
"is_a" : [
{
"$ref" : "go",
"$id" : "GO:0030247"
}
]
}
{
"_id" : "GO:0080115",
"name" : "myosin XI tail binding",
"is_a" : [
{
"$ref" : "go",
"$id" : "GO:0032029"
}
]
}
{
"_id" : "GO:0090079",
"name" : "translation regulator activity, nucleic acid binding",
"is_a" : [
{
"$ref" : "go",
"$id" : "GO:0003676"
},
{
"$ref" : "go",
"$id" : "GO:0045182"
}
]
}
(...)




That's it

Pierre

23 February 2010

A Mysql user defined function (UDF) for Gene Ontology (GO)

In this post I'll create a Mysql Defined function (UDF) answering if a defined GO term is a descendant of another term. This post is not much different from two previous posts I wrote here:

Here I built a binary file containing an array of pairs(term,parent_term) of GO identifiers from the XML/RDF file o_daily-termdb.rdf-xml.gz. This file is then used by the mysql UDF to find all the ancestors of a given term. The source code is available here:

Building the database

The libxml API is used to load go_daily-termdb.rdf-xml.gz. For each go:term, the parent terms are searched and each pair(term,parent-term) is saved into the array.
static void scanterm(xmlNode *term)
{
xmlAttrPtr rsrc;
xmlNode *n=NULL;
xmlChar * acn=NULL;
for(n=term->children; n!=NULL; n=n->next)
{
if (n->type != XML_ELEMENT_NODE) continue;
if(strcmp(n->name,"accession")==0 &&
n->ns!=NULL && n->ns->href!=NULL &&
strcmp(n->ns->href,GO_NS)==0)
{
acn= xmlNodeGetContent(n);
break;
}
}
if(acn==NULL) return;

for(n=term->children; n!=NULL; n=n->next)
{
if (n->type != XML_ELEMENT_NODE) continue;
if(strcmp(n->name,"is_a")==0 &&
n->ns!=NULL && n->ns->href!=NULL &&
strcmp(n->ns->href,GO_NS)==0
)
{
xmlChar* is_a=xmlGetNsProp(n,"resource",RDF_NS);
if(is_a!=NULL)
{
char* p=strstr(is_a,"#GO:");
if(p!=NULL)
{
++p;
termdb.terms=(TermPtr)realloc(termdb.terms,sizeof(Term)*(termdb.n_terms+1));
if(termdb.terms==NULL)
{
fprintf(stderr,"out of memory\n");
exit(EXIT_FAILURE);
}
strncpy(termdb.terms[termdb.n_terms].parent,p,MAX_TERM_LENGTH);
strncpy(termdb.terms[termdb.n_terms].child,acn,MAX_TERM_LENGTH);
++termdb.n_terms;
}
xmlFree(is_a);
}
}
}
xmlFree(acn);
}
When the RDF file has been read, the content of the array is sorted by term and saved to a file
qsort(termdb.terms,termdb.n_terms,sizeof(Term),cmp_term);
out= fopen(argv[1],"w");
if(out==NULL)
{
fprintf(stderr,"Cannot open %s\n",argv[1]);
return -1;
}
fwrite(&(termdb.n_terms),sizeof(int),1,out);
fwrite(termdb.terms,sizeof(Term),termdb.n_terms,out);
fflush(out);
fclose(out);

Initializing the mysql UDF

When the mysql UDF is inited the binary file is loaded into memory (error handling ignored)
my_bool go_isa_init(
UDF_INIT *initid,
UDF_ARGS *args,
char *message
)
{
TermDBPtr termdb;
FILE* in=NULL;

initid->maybe_null=1;
initid->ptr= NULL;

termdb=(TermDBPtr)malloc(sizeof(TermDB));


in=fopen(GO_PATH,"r");
fread(&(termdb->n_terms),sizeof(int),1,in);
termdb->terms=malloc(sizeof(Term)*(termdb->n_terms));
fread(termdb->terms,sizeof(Term),termdb->n_terms,in);
fclose(in);
initid->ptr=(void*)termdb;
return 0;
}

Disposing the mysql UDF

When the UDF is disposed, we just free the memory
/* The deinitialization function */
void go_isa_deinit(UDF_INIT *initid)
{
/* free the memory **/
if(initid->ptr!=NULL)
{
TermDBPtr termdb=(TermDBPtr)initid->ptr;
if(termdb->terms!=NULL) free(termdb->terms);
free(termdb);
initid->ptr=NULL;
}
}

invoking the UDF

The UDF itself will scan the GO directed tree and will get all the ancestors of a given term:
long long go_isa(UDF_INIT *initid, UDF_ARGS *args,
char *is_null, char *error)
{
(...)
index=termdb_findIndexByName(termdb,term);
if(index==-1)
{
return 0;
}
return recursive_search(termdb,index,parent);
}

static int recursive_search(const TermDBPtr db,int index, const char* parent)
{
int rez=0;
int start=index;
int parent_idx=0;

if(start<0 || start>=db->n_terms) return 0;
if(strcmp(db->terms[index].child,parent)==0) return 1;
while(index < db->n_terms)
{
if(strcmp(db->terms[index].child,db->terms[start].child)!=0) break;
if(strcmp(db->terms[index].parent,parent)==0) return 1;
parent_idx= termdb_findIndexByName(db,db->terms[index].parent);
rez= recursive_search(db,parent_idx,parent);
if(rez==1) return 1;
++index;
}
return 0;
}
As the terms have been ordered by their names, a binary_search is used to find the index of a given term:
static int lower_bound(const TermDBPtr termsdb, const char* name)
{
int low = 0;
int len= termsdb->n_terms;

while(len>0)
{
int half=len/2;
int mid=low+half;
if( strncmp(termsdb->terms[mid].child,name,MAX_TERM_LENGTH)<0)
{
low=mid;
++low;
len=len-half-1;
}
else
{
len=half;
}
}
return low;
}


static int termdb_findIndexByName(const TermDBPtr termsdb,const char* name)
{
int i=0;
if(name==NULL || termsdb==NULL || termsdb->terms==NULL || termsdb->n_terms==0) return -1;
i= lower_bound(termsdb,name);
if(i<0 || i >= termsdb->n_terms || strcmp(termsdb->terms[i].child,name)!=0) return -1;
return i;
}

Compiling

On my laptop, the UDF was compiled using
gcc -shared -DGO_PATH='"/tmp/terms.bin"' -I /usr/include -I /usr/local/include -I /usr/include/mysql -o /usr/lib/go_udf.so src/go_udf.c

Creating the function

mysql> create function go_isa RETURNS INTEGER SONAME 'go_udf.so';

Invoking the UDF function

Finf all the GO terms that are a descendant of GO:0016859 (cis-trans isomerase activity):
mysql> select LEFT(T.name,20) as name,T.acc from go_latest.term as T where go_isa(T.acc,"GO:0016859");
+----------------------+------------+
| name | acc |
+----------------------+------------+
| peptidyl-prolyl cis- | GO:0003755 |
| retinal isomerase ac | GO:0004744 |
| maleylacetoacetate i | GO:0016034 |
| cis-trans isomerase | GO:0016859 |
| cis-4-[2-(3-hydroxy) | GO:0018839 |
| trans-geranyl-CoA is | GO:0034872 |
| carotenoid isomerase | GO:0046608 |
| 2-chloro-4-carboxyme | GO:0047466 |
| 4-hydroxyphenylaceta | GO:0047467 |
| farnesol 2-isomerase | GO:0047885 |
| furylfuramide isomer | GO:0047907 |
| linoleate isomerase | GO:0050058 |
| maleate isomerase ac | GO:0050076 |
| maleylpyruvate isome | GO:0050077 |
| retinol isomerase ac | GO:0050251 |
+----------------------+------------+
15 rows in set (1.27 sec)

Disposing the UDF function

drop function go_isa;


That's it
Pierre

05 January 2010

My Python notebook: displaying the state of the Genome Projects

This post if about my first Python program. I'm still a newbie so please, don't flame ! :-)

The NCBI hosts a XML file describing the state of the Genome Projects at ftp://ftp.ncbi.nih.gov/genomes/genomeprj/gp.xml. The very first lines of the file look like this:

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE gp:DocumentSet SYSTEM "gp4v.dtd">
<gp:DocumentSet xmlns:SOAP-ENV="http://schemas.xmlsoap.org/soap/envelope/" xmlns:SOAP-ENC="http://s
chemas.xmlsoap.org/soap/encoding/" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:gp="gp">
<gp:Document>
<gp:ProjectID>1</gp:ProjectID>
<gp:IsPrivate>false</gp:IsPrivate>
<gp:eType>eGenome</gp:eType>
<gp:eMethod>eWGS</gp:eMethod>
<gp:eTechnologies>eUnspecified</gp:eTechnologies>
<gp:OriginDB>NCBI</gp:OriginDB>
<gp:CreateDate>03/05/2003 15:03:16</gp:CreateDate>
<gp:ModificationDate>07/18/2006 16:09:02</gp:ModificationDate>
<gp:ProjectName/>
<gp:OrganismName>Acidobacterium capsulatum ATCC 51196</gp:OrganismName>
<gp:StrainName>ATCC 51196</gp:StrainName>
<gp:TaxID>240015</gp:TaxID>
<gp:LocusTagPrefix/>
<gp:DNASource/>
<gp:SequencingDepth/>
<gp:ProjectURL>http://www.tigr.org/tdb/mdb/mdbinprogress.html</gp:ProjectURL>
<gp:DataURL/>
<gp:OrganismDescription/>
<gp:EstimatedGenomeSize>4.15</gp:EstimatedGenomeSize>
(...)
<gp:ChromosomeDefinitions/>
<gp:SubmittedSequences/>
<gp:LegacyPrefixes/>
</gp:Document>
<gp:Document>
<gp:ProjectID>3</gp:ProjectID>
<gp:IsPrivate>false</gp:IsPrivate>
<gp:eType>eGenome</gp:eType>
(...)
The program I wrote reads this XML file and generates the following visualization:
The left part is the tree of the NCBI taxonomy. In front of each taxon, the circles on the right are the Genome Projects. The size of a circle is the log10 of the estimated size of the genome (if available). The position on the x-axis of those circles is the date given by the tag "<p:CreateDate>".
There are two main classes: A Taxon is a node in the NCBI taxonomy, it contains a link to his children and to his parent:
class Taxon:
def __init__(self):
self.id=-1
self.name=None
self.parent=None
self.children=set()
. A Project is an item in the XML file. It contains a link to a Taxon:
class Project:
def __init__(self):
self.id=-1
self.taxon=None
self.genomeSize=None
self.date=None
The XML file is processed with a SAX handler. Each time a new project is found, the corresponding taxon is downloaded :
class GPHandler(handler.ContentHandler):
LIMIT=1E7
SLEEP=0
def __init__(self,owner):
self.owner=owner
self.content=None
self.project=None

(...)

def startElementNS(self,name,qname,attrs):
if name[1]=="Document":
self.project=Project()
elif name[1] in ["ProjectID","TaxID","CreateDate","EstimatedGenomeSize"] :
self.content=""

def endElementNS(self,name,qname):
if name[1]=="Document":
if len(self.owner.id2project)<GPHandler.LIMIT:
self.owner.id2project[self.project.id]= self.project
self.project=None
elif name[1]=="TaxID":
if len(self.owner.id2project)<GPHandler.LIMIT:
self.project.taxon=self.owner.findTaxonById(int(self.content))
time.sleep(GPHandler.SLEEP) #don't be evil with NCBI
elif name[1]=="ProjectID":
self.project.id=int(self.content)
sys.stderr.write("Project id:"+self.content+" n="+str(len(self.owner.id2project))+"\n")
elif name[1]=="EstimatedGenomeSize" and len(self.content)>0:
self.project.genomeSize=float(self.content)
elif name[1]=="CreateDate":
self.project.date=self.sql2date(self.content)
self.content=None
def characters(self,content):
if self.content!=None:
self.content+=content
The XML description of each taxons is downloaded and parsed using the NCBI EFetch service. A XML for a taxon looks like this:
<TaxaSet>
<Taxon>
<TaxId>240015</TaxId>
<ScientificName>Acidobacterium capsulatum ATCC 51196</ScientificName>
<OtherNames>
<EquivalentName>Acidobacterium capsulatum strain ATCC 51196</EquivalentName>
<EquivalentName>Acidobacterium capsulatum str. ATCC 51196</EquivalentName>
</OtherNames>
<ParentTaxId>33075</ParentTaxId>
<Rank>no rank</Rank>
<Division>Bacteria</Division>
<GeneticCode>
<GCId>11</GCId>
<GCName>Bacterial, Archaeal and Plant Plastid</GCName>
</GeneticCode>
<MitoGeneticCode>
<MGCId>0</MGCId>
<MGCName>Unspecified</MGCName>
</MitoGeneticCode>
<Lineage>cellular organisms; Bacteria; Fibrobacteres/Acidobacteria group; Acidobacteria; Acidobacteria (class); Acidobacteriales; Acidobacteriaceae; Acidobacterium; Acidobacterium capsulatum</Lineage>
<LineageEx>
<Taxon>
<TaxId>131567</TaxId>
<ScientificName>cellular organisms</ScientificName>
<Rank>no rank</Rank>
</Taxon>
<Taxon>
<TaxId>2</TaxId>
<ScientificName>Bacteria</ScientificName>
<Rank>superkingdom</Rank>
</Taxon>
<Taxon>
<TaxId>131550</TaxId>
<ScientificName>Fibrobacteres/Acidobacteria group</ScientificName>
<Rank>superphylum</Rank>
</Taxon>
<Taxon>
<TaxId>57723</TaxId>
<ScientificName>Acidobacteria</ScientificName>
<Rank>phylum</Rank>
</Taxon>
<Taxon>
<TaxId>204432</TaxId>
<ScientificName>Acidobacteria (class)</ScientificName>
<Rank>class</Rank>
</Taxon>
<Taxon>
<TaxId>204433</TaxId>
<ScientificName>Acidobacteriales</ScientificName>
<Rank>order</Rank>
</Taxon>
<Taxon>
<TaxId>204434</TaxId>
<ScientificName>Acidobacteriaceae</ScientificName>
<Rank>family</Rank>
</Taxon>
<Taxon>
<TaxId>33973</TaxId>
<ScientificName>Acidobacterium</ScientificName>
<Rank>genus</Rank>
</Taxon>
<Taxon>
<TaxId>33075</TaxId>
<ScientificName>Acidobacterium capsulatum</ScientificName>
<Rank>species</Rank>
</Taxon>
</LineageEx>
<CreateDate>2003/07/25</CreateDate>
<UpdateDate>2005/01/19</UpdateDate>
<PubDate>2005/01/29</PubDate>
</Taxon>
</TaxaSet>
This time I used a DOM implementation for python to analyse the file:
def findTaxonById(self,taxId):
if taxId in self.id2taxon:
return self.id2taxon[taxId]

url="http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&retmode=xml&id="+str(taxId)
new_taxon=Taxon()
new_taxon.id=taxId
new_taxon.parent

dom = xml.dom.minidom.parse(urllib.urlopen(url))
top= dom.documentElement
e1= self.first(top,"Taxon")
new_taxon.name = self.textContent(self.first(e1,"ScientificName"))

lineage=[]
lineage.append(self.taxonRoot)
for e2 in self.elements(self.first(e1,"LineageEx"),"Taxon"):
t2=Taxon()
t2.id= int(self.textContent(self.first(e2,"TaxId")))
t2.name= self.textContent(self.first(e2,"ScientificName"))
if t2.id in self.id2taxon:
t2= self.id2taxon[t2.id]
else:
self.id2taxon[t2.id]=t2
lineage.append(t2)
lineage.append(new_taxon)

i=1
while i < len(lineage):
lineage[i-1].children.add(lineage[i])
lineage[i].parent=lineage[i-1]
i+=1

self.id2taxon[new_taxon.id]=new_taxon
return new_taxon
At the end, the figure is generated using SVG:

Source code

#Author Pierre Lindenbaum
#Mail: plindenbaum@yahoo.fr
#WWW: http://plindenbaum.blogspot.com
#usage:
# wget -O gp.xml "ftp://ftp.ncbi.nih.gov/genomes/genomeprj/gp.xml"
# python gp.py gp.xml
from xml.sax import make_parser, handler,saxutils
import sys,time,math,pickle,os.path
import xml.dom.minidom
import urllib


class SVG:
NS="http://www.w3.org/2000/svg"
class XLINK:
NS="http://www.w3.org/1999/xlink"



class Main:
DB_FILE="gp.db"
def __init__(self):
self.id2project=dict()
self.id2taxon=dict()
self.min_size=1.0E13
self.max_size=0.0
self.min_date=1.0E13
self.max_date=0.0
self.taxonRoot=Taxon()
self.taxonRoot.id=0
self.taxonRoot.name="Tree of Life"
self.id2taxon[self.taxonRoot.id]=self.taxonRoot

def escape(self,s):
x=""
for c in s:
if c=='<':
x+="&lt;"
elif c=='>':
x+="&gt;"
elif c=='&':
x+="&amp;"
elif c=='\'':
x+="&apos;"
elif c=='\'':
x+="&quot;"
else:
x+=c
return x

def first(self,root,name):
for c1 in root.childNodes:
if c1.nodeType!= xml.dom.minidom.Node.ELEMENT_NODE or \
c1.nodeName!=name:
continue
return c1
return None
def elements(self,root,name):
a=[]
for c1 in root.childNodes:
if c1.nodeType!= xml.dom.minidom.Node.ELEMENT_NODE or \
c1.nodeName!=name:
continue
a.append(c1)
return a

def textContent(self,root):
content=""
for c1 in root.childNodes:
if c1.nodeType== xml.dom.minidom.Node.TEXT_NODE:
content+= c1.nodeValue
else:
content+=self.textContent(c1)
return content

def findTaxonById(self,taxId):
if taxId in self.id2taxon:
return self.id2taxon[taxId]

url="http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&retmode=xml&id="+str(taxId)
sys.stderr.write(url+"\n")

new_taxon=Taxon()
new_taxon.id=taxId
new_taxon.parent

dom = xml.dom.minidom.parse(urllib.urlopen(url))
top= dom.documentElement
e1= self.first(top,"Taxon")
new_taxon.name = self.textContent(self.first(e1,"ScientificName"))

lineage=[]
lineage.append(self.taxonRoot)
for e2 in self.elements(self.first(e1,"LineageEx"),"Taxon"):
t2=Taxon()
t2.id= int(self.textContent(self.first(e2,"TaxId")))
t2.name= self.textContent(self.first(e2,"ScientificName"))
if t2.id in self.id2taxon:
t2= self.id2taxon[t2.id]
else:
self.id2taxon[t2.id]=t2
lineage.append(t2)
lineage.append(new_taxon)

i=1
while i < len(lineage):
lineage[i-1].children.add(lineage[i])
lineage[i].parent=lineage[i-1]
i+=1

self.id2taxon[new_taxon.id]=new_taxon
return new_taxon

def svg(self):
width= 1000
height=10000
for p in self.id2project.values():
if p.genomeSize!=None:
self.min_size = min(self.min_size, p.genomeSize)
self.max_size = max(self.max_size, p.genomeSize)
self.min_date = min(self.min_date, p.date)
self.max_date = max(self.max_date, p.date)
shift= (self.max_date-self.min_date)*0.05
self.min_date-=shift
self.max_date+=shift
sys.stderr.write(" size:"+str(self.min_size)+"/"+str(self.max_size))
sys.stderr.write(" date:"+str(self.min_date)+"/"+str(self.max_date))
print "<?xml version='1.0' encoding='UTF-8'?>"
print "<svg xmlns='"+SVG.NS+"' xmlns:xlink='"+XLINK.NS+"' "+\
"width='"+str(width+1)+"' "+\
"height='"+str(height+1)+"' "+\
"style='stroke:black;fill:none;stroke-width:0.5px;' "+\
"><title>Genome Projects</title>"
print "<rect x='0' y='0' width='"+str(width/2.0)+"' height='"+str(height)+"' style='stroke:green;'/>"
print "<rect x='"+str(1+width/2.0)+"' y='0' width='"+str(width/2.0)+"' height='"+str(height)+"' style='stroke:blue;'/>"
self.recurs(self.taxonRoot,0,0,width/2,height,(width/2.0)/self.taxonRoot.depth())
print "</svg>"
def recurs(self,taxon,x,y,width,height,h_size):
total=0.0
midy=y+height/2.0
print "<circle cx='"+str(x)+"' cy='"+str(midy)+"' r='"+str(2)+"' title='"+self.escape(taxon.name)+"' style='fill:black;'/>"

for p in self.id2project.values():
if p.taxon!=taxon:
continue
cx= width+ ((p.date - self.min_date)/(self.max_date - self.min_date))*width
cy= midy
print "<a target='ncbi' title='"+self.escape(taxon.name)+" projectId:"+str(p.id)+" size:"+str(p.genomeSize)+"' xlink:href='http://www.ncbi.nlm.nih.gov/sites/entrez?Db=genomeprj&amp;cmd=ShowDetailView&amp;TermToSearch="+str(p.id)+"'>"
print "<g style='stroke:red;' >"
if p.genomeSize!=None:
radius=max(2.0,((math.log10(p.genomeSize) - math.log10(self.min_size))/(math.log10(self.max_size) - math.log10(self.min_size)))*100.0)
print "<circle cx='"+str(cx)+"' cy='"+str(cy)+"' r='"+str(radius)+"' style='fill:orange;fill-opacity:0.1;'/>"
print "<line x1='"+str(cx-5)+"' y1='"+str(cy)+"' x2='"+str(cx+5)+"' y2='"+str(cy)+"'/>"
print "<line x1='"+str(cx)+"' y1='"+str(cy-5)+"' x2='"+str(cx)+"' y2='"+str(cy+5)+"'/>"
print "</g></a>"

for c in taxon.children:
total+= c.weight()
for c in taxon.children:
h2=height* (c.weight()/total)
y2=y+h2/2.0
w2=h_size
if len(c.children)==0:
w2=width-x
print "<polyline points='"+\
str(x)+","+str(midy)+" "+\
str(x+w2/2.0)+","+str(midy)+" "+\
str(x+w2/2.0)+","+str(y2)+" "+\
str(x+w2)+","+str(y2)+" "+\
"' title='"+self.escape(c.name)+"'/>"
self.recurs(c,x+w2,y,width,h2,w2)
y+=h2


class Taxon:
def __init__(self):
self.id=-1
self.name=None
self.parent=None
self.children=set()
def __hash__(self):
return self.id.__hash__()
def __eq__(self,other):
if other==None:
return False
return self.id==other.id
def weight(self):
w=1
for t in self.children:
w+= t.weight()
return w
def depth(self):
##sys.stderr.write(self.name+"\n")
d=0.0
for t in self.children:
d= max(t.depth(),d)
return 1.0+d

class Project:
def __init__(self):
self.id=-1
self.taxon=None
self.genomeSize=None
self.date=None

class GPHandler(handler.ContentHandler):
LIMIT=1E7
SLEEP=0
def __init__(self,owner):
self.owner=owner
self.content=None
self.project=None
def sql2date(self,s):
s=s[:s.index(' ')]
cal =time.mktime( time.strptime(s,"%m/%d/%Y") )
return cal
def startElementNS(self,name,qname,attrs):
if name[1]=="Document":
self.project=Project()
elif name[1] in ["ProjectID","TaxID","CreateDate","EstimatedGenomeSize"] :
self.content=""
def endElementNS(self,name,qname):
if name[1]=="Document":
if len(self.owner.id2project)<GPHandler.LIMIT:
self.owner.id2project[self.project.id]= self.project
self.project=None
elif name[1]=="TaxID":
if len(self.owner.id2project)<GPHandler.LIMIT:
self.project.taxon=self.owner.findTaxonById(int(self.content))
time.sleep(GPHandler.SLEEP) #don't be evil with NCBI
elif name[1]=="ProjectID":
self.project.id=int(self.content)
sys.stderr.write("Project id:"+self.content+" n="+str(len(self.owner.id2project))+"\n")
elif name[1]=="EstimatedGenomeSize" and len(self.content)>0:
self.project.genomeSize=float(self.content)
elif name[1]=="CreateDate":
self.project.date=self.sql2date(self.content)
self.content=None
def characters(self,content):
if self.content!=None:
self.content+=content




main=Main();
sys.stderr.write(main.escape("<>!")+"\n")
parser = make_parser()
# we want XML namespaces
parser.setFeature(handler.feature_namespaces,True)
# tell parser to disable validation
parser.setFeature(handler.feature_validation,False)
parser.setFeature(handler.feature_external_ges, False)
parser.setContentHandler(GPHandler(main))
parser.parse(sys.argv[1])

main.svg()


That's it,

Pierre

21 March 2007

Geni, Graphiz, Dot & Family Tree.

Geni is a genealogy-related social networking website launched in beta mode in January 2007. Since yesterday, the family tree can now be exported as a gedcom-xml file (alpha version).

The following xslt stylesheet transforms the gedcom file into a Graphiz/DOT input which can be used to generate the family tree.

Family Tree



<?xml version='1.0' ?>
<xsl:stylesheet xmlns:xsl='http://www.w3.org/1999/XSL/Transform' version='1.0'>
<xsl:output method='text' omit-xml-declaration="yes" />

<xsl:template match="GEDCOM">
digraph &quot;G&quot; {
<xsl:apply-templates select="FamilyRec"/>
<xsl:apply-templates select="IndividualRec"/>
}
</xsl:template>

<xsl:template match="IndividualRec">
<xsl:value-of select="@Id"/>[ shape=box, label=&quot;<xsl:value-of select="IndivName/GivenName"/><xsl:text> </xsl:text><xsl:value-of select="IndivName/SurName"/> <xsl:if test="DeathStatus=&apos;dead;&apos;">(d)</xsl:if>&quot;

<xsl:choose>
<xsl:when test="Gender=&apos;M&apos;">
,color=blue
</xsl:when>
<xsl:when test="Gender=&apos;F&apos;">
,color=pink
</xsl:when>
<xsl:otherwise>
,color=black
</xsl:otherwise>
</xsl:choose>


];
</xsl:template>

<xsl:template match="FamilyRec">
<xsl:variable name="famId"><xsl:value-of select="@Id"/></xsl:variable>

<xsl:value-of select="$famId"/>[shape=point];

<xsl:if test="HusbFath">
<xsl:value-of select="HusbFath/Link/@Ref"/>-&gt;<xsl:value-of select="$famId"/>;
</xsl:if>

<xsl:if test="WifeMoth">
<xsl:value-of select="WifeMoth/Link/@Ref"/>-&gt;<xsl:value-of select="$famId"/>;
</xsl:if>

<xsl:for-each select="Child">
<xsl:value-of select="$famId"/>-&gt;<xsl:value-of select="Link/@Ref"/>;
</xsl:for-each>
</xsl:template>
</xsl:stylesheet>



Pierre

28 January 2007

Social Genealogy

From the FAQ: Geni is a unique approach to solving the problem of genealogy, which is the question of how everyone is related. Geni lets you create a family tree through our fun simple interface. When you add a relative's email address, he or she will be invited to join your tree. That relative can then add other relatives, and so on. Your tree will continue to grow as relatives invite other relatives. Each family member has a profile which can be viewed by clicking their name in the tree. This helps family members learn more about each other and stay in touch. Family members can also share information and work together to build profiles for common ancestors. Geni is a private network. Only the people in your tree can see your tree and your profile. Geni will not share your personal information with third parties.



See also, my previous post about how to draw pedigrees using DOT.

Pierre