Showing posts with label sax. Show all posts
Showing posts with label sax. Show all posts

11 October 2010

Playing with the Wordle algorithm: a tag cloud of Mesh Terms

The paper describing Wordle has been recently published (http://www.research.ibm.com/visual/papers/wordle_final2.pdf ). The algorithm was briefly described: “The most distinctive geometric aspect of a Wordle is the layout algorithm, which packs words to make efficient use of space. While many space-filling visualizations exist, they typically work by recursiveLayout proceeds according to this pseudocode:

sort words by weight, decreasing
for each word w:
w.position := makeInitialPosition(w);
while w intersects other words:
updatePosition(w);

The two key procedures here are "makeInitialPosition" and "updatePosition". The makeInitialPosition routine picks a point at random according to a distribution that takes into account the desired overall shape, and, if desired, alphabetical order. The updatePosition routine moves the word on a spiral of increasing radius, radiating from the word's starting position. The updatePosition routine is aware of constraints on the overall shape of the Wordle. Constraining the layout to a rectangular shape causes updatePosition to prefer positions inside of the strict boundaries of the playing field; a blobby overall shape accepts boundary violations. The rectangular constraint is relaxed when the spiral radius exceeds either playing field dimension. ”
Your browser does not support the <CANVAS> element !


Just for fun , I've implemented my own version of the Wordle Alogrithm. The Java code is available on github at http://github.com/lindenb/jsandbox/blob/master/src/sandbox/MyWordle.java. I won't describe the program here, but I'll just say that the code invokes a java.awt.font.TextLayout class to get the shape of the text:
Graphics2D g=(...)
FontRenderContext frc = g.getFontRenderContext();
Font font=new Font("Dialog",Font.BOLD,fontSize);
TextLayout textLayout=new TextLayout(w.getText(), font, frc);
Shape shape=textLayout.getOutline(null);
and a java.awt.geom.Area to test if two shapes intersects.

Ok, let's test this code. FIrst I'm going to dump a pubmed query as XML with another simple tool named PubmedDump. This XML file is then parsed with a javascript program called from the SAX parser saxstream.jar (previously described here):
mesh.js
importPackage(Packages.sandbox);
importPackage(Packages.java.io);
importPackage(Packages.java.awt);
var content=null;
var mesh2count={};

function startElement(uri,localName,name,atts)
{
if(name=="DescriptorName")
{
content="";
}
}
function characters(s)
{
if(content!=null) content+=s;
}
function endElement(uri,localName,name)
{
if(content!=null)
{
var count=mesh2count[content];
if(count===undefined) count=0;
mesh2count[content]=count+1;
}
content=null;
}
function endDocument()
{
var w= new MyWordle();
for(var s in mesh2count)
{
var word= new MyWordle.Word(s,mesh2count[s]);
w.add(word);
}
w.setUseArea(true);/* use shape area instead of bounding boxes */
w.setAllowRotate(true);
w.setSortType(1);/* sort by weight */
w.doLayout();

var f=new File("result.svg");
w.saveAsSVG(f);
}
This javascript code counts the occurence of each MESH term (<DescriptorName>), and when the document has been parsed, a new instance of MyWordle class is created, filled and the result is saved to a SVG file.

Invocation


Here is an example for the query "Rotavirus NSP3 NSP1":
java -jar pubmeddump.jar "Rotavirus NSP3 NSP1" |\
java -cp mywordle.jar:saxscript.jar org.lindenb.tinytools.SAXScript -n -f mesh.js


Result





That's it,
Pierre

12 February 2010

Processing large XML documents with XSLT

I've resurrected an old java program called xsltstream which might be useful for biohackathon2010. This program applies a XSLT stylesheet only for the given node from a large xml document. The DOM is read from a SAX stream, built in memory for each target element , processed with XSLT and then disposed. Now, say you want to transform a XML file from dbSNP with XSLT to make a RDF document. You cannot do that with xsltproc because the XML file is just too big ( e.g. ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/XML/ds_ch1.xml.gz is 1,099,375 KB ).
But an xslt stylesheet can be applied with xsltstream to all the <Rs> elements of 'ds_ch1.xml.gz':

java -jar xsltstream.jar -x 'http://lindenb.googlecode.com/svn/trunk/src/xsl/dbsnp2rdf.xsl' -q Rs \
'ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/XML/ds_ch1.xml.gz' |\
grep -v "rdf:RDF" | grep -v "<?xml version"


(...)
<o:SNP rdf:about="http://www.ncbi.nlm.nih.gov/snp/830">
<dc:title>rs830</dc:title>
<o:taxon rdf:resource="http://www.ncbi.nlm.nih.gov/taxonomy/9606"/>
<o:het rdf:datatype="http://www.w3.org/2001/XMLSchema#float">0.02</o:het>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:WIAF"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:SNP500CANCER"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:SEQUENOM"/>
<o:hasMapping>
<o:Mapping>
<o:build rdf:resource="urn:void:ncbi:build:Celera/36_3"/>
<o:chrom rdf:resource="urn:void:ncbi:chromosome:9606/chr1"/>
<o:start rdf:datatype="http://www.w3.org/2001/XMLSchema#int">66444409</o:start>
<o:end rdf:datatype="http://www.w3.org/2001/XMLSchema#int">66444410</o:end>
<o:orient>+</o:orient>
</o:Mapping>
</o:hasMapping>
<o:hasMapping>
<o:Mapping>
<o:build rdf:resource="urn:void:ncbi:build:HuRef/36_3"/>
<o:chrom rdf:resource="urn:void:ncbi:chromosome:9606/chr1"/>
<o:start rdf:datatype="http://www.w3.org/2001/XMLSchema#int">66263806</o:start>
<o:end rdf:datatype="http://www.w3.org/2001/XMLSchema#int">66263807</o:end>
<o:orient>-</o:orient>
</o:Mapping>
</o:hasMapping>
<o:hasMapping>
<o:Mapping>
<o:build rdf:resource="urn:void:ncbi:build:reference/36_3"/>
<o:chrom rdf:resource="urn:void:ncbi:chromosome:9606/chr1"/>
<o:start rdf:datatype="http://www.w3.org/2001/XMLSchema#int">67926134</o:start>
<o:end rdf:datatype="http://www.w3.org/2001/XMLSchema#int">67926135</o:end>
<o:orient>+</o:orient>
</o:Mapping>
</o:hasMapping>
</o:SNP>


<o:SNP rdf:about="http://www.ncbi.nlm.nih.gov/snp/844">
<dc:title>rs844</dc:title>
<o:taxon rdf:resource="http://www.ncbi.nlm.nih.gov/taxonomy/9606"/>
<o:het rdf:datatype="http://www.w3.org/2001/XMLSchema#float">0.42</o:het>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:WIAF"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:LEE"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:HGBASE"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:SC_JCM"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:TSC-CSHL"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:LEE"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:YUSUKE"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:CGAP-GAI"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:CSHL-HAPMAP"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:PERLEGEN"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:ABI"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:SI_EXO"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:BCMHGSC_JDW"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:HUMANGENOME_JCVI"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:SNP500CANCER"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:1000GENOMES"/>
<o:hasHandle rdf:resource="urn:void:ncbi:snp:handle:ILLUMINA-UK"/>
<o:hasMapping>
<o:Mapping>
<o:build rdf:resource="urn:void:ncbi:build:Celera/36_3"/>
<o:chrom rdf:resource="urn:void:ncbi:chromosome:9606/chr1"/>
<o:start rdf:datatype="http://www.w3.org/2001/XMLSchema#int">134750981</o:start>
<o:end rdf:datatype="http://www.w3.org/2001/XMLSchema#int">134750982</o:end>
<o:orient>+</o:orient>
</o:Mapping>
</o:hasMapping>
<o:hasMapping>
<o:Mapping>
<o:build rdf:resource="urn:void:ncbi:build:HuRef/36_3"/>
<o:chrom rdf:resource="urn:void:ncbi:chromosome:9606/chr1"/>
<o:start rdf:datatype="http://www.w3.org/2001/XMLSchema#int">132892081</o:start>
<o:end rdf:datatype="http://www.w3.org/2001/XMLSchema#int">132892082</o:end>
<o:orient>-</o:orient>
(...)


The java archive for xsltstream is available at http://lindenb.googlecode.com/files/xsltstream.jar

Usage:
  -x <xslt-stylesheet file/url> required
-p <param-name> <param-value> (add parameter to the xslt engine)
-d depth (0 based) default:-1
-q qName target default:null
<file>|stdin



That's it !
Pierre

07 February 2010

I Really need to sleep: inserting the SNPs into MongoDB with C++

In the previous post I showed how to parse dbSNP/XML with libxml. In the current post I'll insert the results into MongoDB using the native Mongo C++ API. Some others have already posted about MongoDB: for example see Jan's, Brad's or Neil's posts The Boost-C++library is required: my code failed to compile with boost 4.* but it compiled fine with boost 3.9.

Starting mongoDB

> mkdir ~/tmp/MONGODB/data
> mongodb-linux-i686-1.2.2/bin/mongod -dbpath ~/tmp/MONGODB/data

Sun Feb 7 19:22:36 Mongo DB : starting : pid = 12123 port = 27017 dbpath = /home/pierre/tmp/MONGODB/data master = 0 slave = 0 32-bit
** NOTE: when using MongoDB 32 bit, you are limited to about 2 gigabytes of data
** see http://blog.mongodb.org/post/137788967/32-bit-limitations for more

Sun Feb 7 19:22:36 db version v1.2.2, pdfile version 4.5
Sun Feb 7 19:22:36 git version: 8a4fb8b1c7cb78648c55368d806ba35054f6be54
Sun Feb 7 19:22:36 sys info: Linux domU-12-31-39-01-70-B4 2.6.21.7-2.fc8xen #1 SMP Fri Feb 15 12:39:36 EST 2008 i686 BOOST_LIB_VERSION=1_37
Sun Feb 7 19:22:36 waiting for connections on port 27017

Refactoring the XML parser for dbSNP


I've added a 'mongo::DBClientConnection connection' to the MondoDB database in the class DBSNPHandler. This connection is simply opened with
state.connection.connect("localhost");
.
Now, in the method "endElement" when a new SNP is found, a new mongo::BSONObjBuilder object is filled with the fields describing a snp. The "rs####" is used as the key of the database.
mongo::BSONObjBuilder b;
b.append("_id", state->rs_id);
b.append("name", state->rs_id);
b.append("seq5", state->seq5);
b.append("observed", state->observed);
b.append("seq3", state->seq3);
This object is then converted to a mongo::BSONObj and inserted into the MondoDB database:
mongo::BSONObj p = b.obj();
state->connection.insert("ncbi.dbsnp", p);
At the end, we can loop over all the items:
std::cout << "count:" << state.connection.count("ncbi.dbsnp") << std::endl;
mongo::BSONObj emptyObj;
std::auto_ptr<mongo::DBClientCursor> cursor = state.connection.query("ncbi.dbsnp", emptyObj);
while( cursor->more() )
{
std::cout << cursor->next().toString() << std::endl;
}

Compilation

Here , the program was compiled using
export LD_LIBRARY_PATH=../boost/boost_1_39_0/stage/lib:../mongodb-linux-i686-1.2.2/lib

g++ `xml2-config --cflags --libs` \
-I ../mongodb-linux-i686-1.2.2/include/mongo\
-I ../boost/boost_1_39_0\
-L ../mongodb-linux-i686-1.2.2/lib \
-L ../boost/boost_1_39_0/stage/lib \
dbsnp.c \
-lmongoclient \
-lboost_thread-gcc43-mt \
-lboost_filesystem-gcc43-mt

Execution


./a.out ds_ch17.xml.gz
count:740
{ _id: "rs69624490", name: "rs69624490", seq5: "CTCCAGCCCGGGCCCACCCTACAGCCGACACCAAGTTGCGTCACCGTGATCTGGACACCCAGACGTACAT...", observed: "C/G", seq3: "GGTACTCCACCAGCAGGAGCAGGAGGCTCTCCCACCTCCACCTGCCACTGGCCGACAGCTCCCAGTGCGT..." }
{ _id: "rs69626771", name: "rs69626771", seq5: "GAAATTCCTACAAAAATCCATCTTTGTGGCATCATCAGGGGTGATCTGTCCTTGCAGGAATATCTCACGT...", observed: "G/T", seq3: "CTGTAGTTAAACTGGTAACGAAGGTTCCACCCTTCCTCCCAGGCCACAGCGCCCCCAGCAGTCTTGCACC..." }
{ _id: "rs69628798", name: "rs69628798", seq5: "ACCTCAGTAAAATGAAGATTATTACTATGTGCCCTATGCAGGACAGGGACTGTGTTCTGATACAGGCCCT...", observed: "A/G", seq3: "TTGAACAGATATCAGAAAAAGGGGGAGAGAGAAATCAGTTGGTTGGGAGGAGAATGAGGGGGGCAGGAGG..." }
{ _id: "rs69633675", name: "rs69633675", seq5: "GTTGTTGGTGAGGGGAGGGAGTGGGGCAAGAGGAACAGTGTGGTCTAGAGGATAGAGCAAGGGAATGGGA...", observed: "A/G", seq3: "TGGCTCAGTGGAAAGAACACGGGCTTTGGAGTCAGAGATCAGGGGTTCGAATCCCGGCTCTGCCACTTGG..." }
{ _id: "rs69638381", name: "rs69638381", seq5: "GGCCACGCTGCTTGTCATAGCGGCTTTCCAGCTCTGCCCTCCGCAGCAAGGGCAGCGCTCACCTAGGCAA...", observed: "G/T", seq3: "GGGCCGTACCGATGAGTTCTCCCGGGGAGAGACCAGGAGCTCTGAGTCAGGAAGGGAATCAAAAGGCGAC..." }
{ _id: "rs69640257", name: "rs69640257", seq5: "TTTCCTCTTCTCCTTGGCCCTGCATTATCCCCACTATTTGACTTGTCCAGGTCAGCCTTTGTAAATAAAA...", observed: "A/G", seq3: "ACTCTTTACCATCACTTATTCTCAAGGCTGGTCAAAACCTTCTCTGTCATTGTCATTTGTTTATTGAGCA..." }
{ _id: "rs69640260", name: "rs69640260", seq5: "CTCAAAAGAGAATGGTTCCTTTAGGTCCCTGAGGACACCCCAGGAAGGTTGGGTATCCCTTGCTTCAATT...", observed: "C/T", seq3: "AGCGTGTCTTAGTGGAAAAAGCACAAGTCTGGGAGTCAGAGTATCTGGGTTCTAATCCCCACCACCCAAT..." }
{ _id: "rs69641743", name: "rs69641743", seq5: "TCTATACATTGTTTAGTTCCTCTCCCCCACTAGACTGTAAACACCTTGAGGGCAAGGAGCATCTCTTCTG...", observed: "A/C", seq3: "CTAGCAAATAGAGCACGGGCCTGGGAATCTGAATAGGTTCTAATCCTGGCTCTGCCACTTGTCTGCTATG..." }
{ _id: "rs69642295", name: "rs69642295", seq5: "AGCATGCACGATATAAAAAATGCCTGAGCACGTTTTCGACCACCTGAGGGAAGCAGAGGAAAGAGTGAAA...", observed: "A/G", seq3: "GGGATAGCTCTTGTGTGGATCAGGAGTGTGGTTCAACAGCACAAATTCATTTACAGAGATCACTAGAAGC..." }
{ _id: "rs69643084", name: "rs69643084", seq5: "AATTCCCTGATACAGTGTTTTGCACAGGGTGCTTTACTGCTGAAAGACTGAATGCTGCCCTATAGCCTGC...", observed: "C/T", seq3: "GTTTTCTTAAGAACCCATGGGGCTCTGGGCTACTTCTTTTCTTTGTGACTCCATAGTAGTTTACAAAAGC..." }
(...)

Source code


#include <libxml/parser.h>
#include <string>
#include <cstring>
#include <iostream>
#include "client/dbclient.h"
#include "db/jsobj.h"

/**
* Hold the state of the parser
*/
class DBSNPHandler
{
public:
mongo::DBClientConnection connection;
// current rs### id
std::string rs_id;
// current 5' sequence
std::string seq5;
// current observed variation
std::string observed;
// current 3' sequence
std::string seq3;
// current string handler by the SAX handler
std::string* content;
//did we find the sequence ?
bool sequence_found;

DBSNPHandler():content(NULL),sequence_found(false)
{
}

~DBSNPHandler()
{
clear();
}

void clear()
{
if(content!=NULL) delete content;
content=NULL;
sequence_found=false;
rs_id.clear();
seq5.clear();
observed.clear();
seq3.clear();
content=NULL;
}
};

/** called when an TAG is opened */
static void startElement(void * ctx,
const xmlChar * localname,
const xmlChar * prefix,
const xmlChar * URI,
int nb_namespaces,
const xmlChar ** namespaces,
int nb_attributes,
int nb_defaulted,
const xmlChar ** attributes)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(strcmp( (char*) localname,"Rs")==0)
{
state->clear();

for(int i=0;i< nb_attributes;++i)
{
if(strcmp((char*) attributes[i*5],"rsId")!=0) continue;
int len=(char*) attributes[i*5+4]-(char*) attributes[i*5+3];
state->rs_id.assign("rs");
state->rs_id.append(
(char*) attributes[i*5+3],
len
);
break;
}
}
else if((strcmp( (char*) localname,"Seq5")==0 ||
strcmp( (char*) localname,"Observed")==0 ||
strcmp( (char*) localname,"Seq3")==0) &&
state->sequence_found==false
)
{
state->content=new std::string;
}
}


/** called when an TAG is closed */
static void endElement(void * ctx,
const xmlChar * localname,
const xmlChar * prefix,
const xmlChar * URI)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(strcmp( (char*) localname,"Rs")==0)
{
mongo::BSONObjBuilder b;
b.append("_id", state->rs_id);
b.append("name", state->rs_id);
b.append("seq5", state->seq5);
b.append("observed", state->observed);
b.append("seq3", state->seq3);
mongo::BSONObj p = b.obj();

state->connection.insert("ncbi.dbsnp", p);


//we're done with this SNP, clear the state
state->clear();
}
else if(state->content!=NULL && strcmp( (char*) localname,"Seq5")==0)
{
state->seq5.assign(*(state->content));
delete state->content;
state->content=NULL;
}
else if(state->content!=NULL && strcmp( (char*) localname,"Observed")==0)
{
state->observed.assign(*(state->content));
delete state->content;
state->content=NULL;
}
else if(state->content!=NULL && strcmp( (char*) localname,"Seq3")==0)
{
state->seq3.assign(*(state->content));
delete state->content;
state->content=NULL;
state->sequence_found=true;
}
}

static void handleCharacters(void * ctx, const xmlChar * ch, int len)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(state->content!=NULL)
{
state->content->append((char*)ch,len);
}
}


int main(int argc,char **argv)
{
int res;
LIBXML_TEST_VERSION
xmlSAXHandler handler;
DBSNPHandler state;

memset(&handler,0,sizeof(xmlSAXHandler));


handler.startElementNs= startElement;
handler.endElementNs=endElement;
handler.characters=handleCharacters;
handler.initialized = XML_SAX2_MAGIC;

state.connection.connect("localhost");

for(int i=1;i< argc;++i)
{
res=xmlSAXUserParseFile(&handler,&state,argv[i]);

if(res!=0)
{
std::cerr << "Error "<< res << argv[i] << std::endl;
}
}

//dump results
std::cout << "count:" << state.connection.count("ncbi.dbsnp") << std::endl;
mongo::BSONObj emptyObj;
std::auto_ptr<mongo::DBClientCursor> cursor = state.connection.query("ncbi.dbsnp", emptyObj);
while( cursor->more() )
{
std::cout << cursor->next().toString() << std::endl;
}

xmlCleanupParser();
xmlMemoryDump();
return(0);
}


Hey... here come the sandman.. :-)


That's it !
Pierre

Id rather bee sleeping: fast parsing of dbsnp/XML with libxml

Hello world ! I'm now in Japan for Biohackathon 2010 it's 00H30 here and, of course, I cannot sleep. I should be tired : I spent 11H00 in a plane, 1H30 in a train and three hours skating in Tokyo with my kick scooter :-)

I'm killing the time with libxml, the C library for parsing XML. In the current post I'll write a simple SAX parser extracting the flanking sequences from the SNPs of dbSNP/xml.
First we need to include a few files:
#include <libxml/parser.h>
#include <string>
#include <cstring>
#include <iostream>
The the C++ class DBSNPHandler that holds the current state of the parser:
class DBSNPHandler
{
public:

// current rs### id
std::string rs_id;
// current 5' sequence
std::string seq5;
// current observed variation
std::string observed;
// current 3' sequence
std::string seq3;
// current string handler by the SAX handler
std::string* content;
//did we find the sequence ?
bool sequence_found;

DBSNPHandler();
~DBSNPHandler();
void clear();
}
I defined a number of callback methods that will be called when events occur during parsing. 3 callbacks are needed here: when the parser opens an element, when the parse closes an element and when it finds some text. Each callback uses a DBSNPHandler as the value of the user data (ctx);
startElement is called when the parser opens an element. If it is a <Rs> tag then we go threw the attributes to find the rs## id. If it is a component of the sequence then we tell the DBSNPHandler that it should store the text in the 'content' variable.
/** called when an TAG is opened */
void startElement(void * ctx,
const xmlChar * localname,
const xmlChar * prefix,
const xmlChar * URI,
int nb_namespaces,
const xmlChar ** namespaces,
int nb_attributes,
int nb_defaulted,
const xmlChar ** attributes)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(strcmp( (char*) localname,"Rs")==0)
{
state->clear();
//loop over the attributes an find the rs###
for(int i=0;i< nb_attributes;++i)
{
if(strcmp((char*) attributes[i*5],"rsId")!=0) continue;
int len=(char*) attributes[i*5+4]-(char*) attributes[i*5+3];

state->rs_id.assign(
(char*) attributes[i*5+3],
len
);
break;
}
}
else if((strcmp( (char*) localname,"Seq5")==0 ||
strcmp( (char*) localname,"Observed")==0 ||
strcmp( (char*) localname,"Seq3")==0) &&
state->sequence_found==false
)
{
//tells the DBSNPHandler we need to get the next text content
state->content=new std::string;
}
}

handleCharacters is called by the SAX parser it finds some text. This text is appended to the current content if we are parsing the sequence component of the current rs####.
static void handleCharacters(void * ctx, const xmlChar * ch, int len)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(state->content!=NULL)
{
state->content->append((char*)ch,len);
}
}

endElement is called when the parser closes an element. If it is a component of the sequences, the fields (seq5,seq3...) of the current state are updated. It it is a <Rs> element, then its sequence is printed out and the state is re-initialized.
static void endElement(void * ctx,
const xmlChar * localname,
const xmlChar * prefix,
const xmlChar * URI)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(strcmp( (char*) localname,"Rs")==0)
{
std::cout << "rs" << state->rs_id <<
"\t" << state->seq5 <<
"[" << state->observed << "]" <<
state->seq3 << std::endl;

//we're done with this SNP, clear the state
state->clear();
}
else if(state->content!=NULL && strcmp( (char*) localname,"Seq5")==0)
{
state->seq5.assign(*(state->content));
delete state->content;
state->content=NULL;
}
else if(state->content!=NULL && strcmp( (char*) localname,"Observed")==0)
{
state->observed.assign(*(state->content));
delete state->content;
state->content=NULL;
}
else if(state->content!=NULL && strcmp( (char*) localname,"Seq3")==0)
{
state->seq3.assign(*(state->content));
delete state->content;
state->content=NULL;
state->sequence_found=true;
}
}
In the main part, the SAX handler is initialized with the previous callbacks and the xml files are parsed:
xmlSAXHandler handler;
DBSNPHandler state;
(...)
handler.startElementNs= startElement;
handler.endElementNs=endElement;
handler.characters=handleCharacters;

for(int i=1;i< argc;++i)
{
xmlSAXUserParseFile(&handler,&state,argv[i]);
}

Compiling

g++ `xml2-config --cflags --libs` dbsnp.c

Running

a.out ds_ch17.xml.gz
rs69624490 CTCCAGCCCGGGCCCACCCTACAGCCGACACCAAGTTGCGTCACCGTGATCTGGACACCCAGACGTACATTAGAGCTGCTTTCCTTGATGAGCTCAGAACC[C/G]GGTACTCCACCAGCAGGAGCAGGAGGCTCTCCCACCTCCACCTGCCACTGGCCGACAGCTCCCAGTGCGTTCTTCAGGCGCCACTTTCTCGCTGGAGAAAA
rs69626771 GAAATTCCTACAAAAATCCATCTTTGTGGCATCATCAGGGGTGATCTGTCCTTGCAGGAATATCTCACGTCCTCTGTTTGTACCTTGACACGCTTGGCTGA[G/T]CTGTAGTTAAACTGGTAACGAAGGTTCCACCCTTCCTCCCAGGCCACAGCGCCCCCAGCAGTCTTGCACCAGATTCGAATTTATGAACCCACAGCACTTGC
rs69628798 ACCTCAGTAAAATGAAGATTATTACTATGTGCCCTATGCAGGACAGGGACTGTGTTCTGATACAGGCCCTGATTAGGTTAGTACAGTGCCTGGCACATAGT[A/G]TTGAACAGATATCAGAAAAAGGGGGAGAGAGAAATCAGTTGGTTGGGAGGAGAATGAGGGGGGCAGGAGGGTCCTGCAGGGTGTTGCAGGTGAGAAATGAC
rs69633675 GTTGTTGGTGAGGGGAGGGAGTGGGGCAAGAGGAACAGTGTGGTCTAGAGGATAGAGCAAGGGAATGGGAGTCAGAAGGACCTGTGTTCTGGAGAAGCAGC[A/G]TGGCTCAGTGGAAAGAACACGGGCTTTGGAGTCAGAGATCAGGGGTTCGAATCCCGGCTCTGCCACTTGGCAGCTGTGTGACTGTGGGCAAGTCACTTCAC
rs69638381 GGCCACGCTGCTTGTCATAGCGGCTTTCCAGCTCTGCCCTCCGCAGCAAGGGCAGCGCTCACCTAGGCAAGCCCAGAGGGCTTAGGAGGGAGGGGCGGGGC[G/T]GGGCCGTACCGATGAGTTCTCCCGGGGAGAGACCAGGAGCTCTGAGTCAGGAAGGGAATCAAAAGGCGACAGGCTCTCATCATCAGTGTTGCCCAGGATCT
rs69640257 TTTCCTCTTCTCCTTGGCCCTGCATTATCCCCACTATTTGACTTGTCCAGGTCAGCCTTTGTAAATAAAAACAACAGGCTTACCAGCCTACCTTTCTCCCA[A/G]ACTCTTTACCATCACTTATTCTCAAGGCTGGTCAAAACCTTCTCTGTCATTGTCATTTGTTTATTGAGCATCCATCTTCTACAGAGCCCTGGAGTAAGCGT
rs69640260 CTCAAAAGAGAATGGTTCCTTTAGGTCCCTGAGGACACCCCAGGAAGGTTGGGTATCCCTTGCTTCAATTATCAGAGGTAAAAACTCTCTCTCTAAAGGAG[C/T]AGCGTGTCTTAGTGGAAAAAGCACAAGTCTGGGAGTCAGAGTATCTGGGTTCTAATCCCCACCACCCAATGCTTGCTGTGTGACCTTGGGCCTCAGTTACC
rs69641743 TCTATACATTGTTTAGTTCCTCTCCCCCACTAGACTGTAAACACCTTGAGGGCAAGGAGCATCTCTTCTGAATCTGTTATATTCTCACAGGAGCAGCATGG[A/C]CTAGCAAATAGAGCACGGGCCTGGGAATCTGAATAGGTTCTAATCCTGGCTCTGCCACTTGTCTGCTATGGGACCTTGGGCAAATGATTTAACTTCTCTGT
(...)

The code

#include <libxml/parser.h>
#include <string>
#include <cstring>
#include <iostream>
/**
* Holds the state of the parser
*/
class DBSNPHandler
{
public:
// current rs### id
std::string rs_id;
// current 5' sequence
std::string seq5;
// current observed variation
std::string observed;
// current 3' sequence
std::string seq3;
// current string handler by the SAX handler
std::string* content;
//did we find the sequence ?
bool sequence_found;

DBSNPHandler():content(NULL),sequence_found(false)
{
}

~DBSNPHandler()
{
clear();
}

void clear()
{
if(content!=NULL) delete content;
content=NULL;
sequence_found=false;
rs_id.clear();
seq5.clear();
observed.clear();
seq3.clear();
content=NULL;
}
};

/** called when an TAG is opened */
static void startElement(void * ctx,
const xmlChar * localname,
const xmlChar * prefix,
const xmlChar * URI,
int nb_namespaces,
const xmlChar ** namespaces,
int nb_attributes,
int nb_defaulted,
const xmlChar ** attributes)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(strcmp( (char*) localname,"Rs")==0)
{
state->clear();

for(int i=0;i< nb_attributes;++i)
{
if(strcmp((char*) attributes[i*5],"rsId")!=0) continue;
int len=(char*) attributes[i*5+4]-(char*) attributes[i*5+3];

state->rs_id.assign(
(char*) attributes[i*5+3],
len
);
break;
}
}
else if((strcmp( (char*) localname,"Seq5")==0 ||
strcmp( (char*) localname,"Observed")==0 ||
strcmp( (char*) localname,"Seq3")==0) &&
state->sequence_found==false
)
{
state->content=new std::string;
}
}


/** called when an TAG is closed */
static void endElement(void * ctx,
const xmlChar * localname,
const xmlChar * prefix,
const xmlChar * URI)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(strcmp( (char*) localname,"Rs")==0)
{
std::cout << "rs" << state->rs_id <<
"\t" << state->seq5 <<
"[" << state->observed << "]" <<
state->seq3 << std::endl;

//we're done with this SNP, clear the state
state->clear();
}
else if(state->content!=NULL && strcmp( (char*) localname,"Seq5")==0)
{
state->seq5.assign(*(state->content));
delete state->content;
state->content=NULL;
}
else if(state->content!=NULL && strcmp( (char*) localname,"Observed")==0)
{
state->observed.assign(*(state->content));
delete state->content;
state->content=NULL;
}
else if(state->content!=NULL && strcmp( (char*) localname,"Seq3")==0)
{
state->seq3.assign(*(state->content));
delete state->content;
state->content=NULL;
state->sequence_found=true;
}
}

static void handleCharacters(void * ctx, const xmlChar * ch, int len)
{
DBSNPHandler* state=(DBSNPHandler*)ctx;
if(state->content!=NULL)
{
state->content->append((char*)ch,len);
}
}


int main(int argc,char **argv)
{
int res;
LIBXML_TEST_VERSION
xmlSAXHandler handler;
DBSNPHandler state;

memset(&handler,0,sizeof(xmlSAXHandler));


handler.startElementNs= startElement;
handler.endElementNs=endElement;
handler.characters=handleCharacters;
handler.initialized = XML_SAX2_MAGIC;

for(int i=1;i< argc;++i)
{
res=xmlSAXUserParseFile(&handler,&state,argv[i]);

if(res!=0)
{
std::cerr << "Error "<< res << argv[i] << std::endl;
}
}
xmlCleanupParser();
xmlMemoryDump();

return(0);
}


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

22 June 2009

Event-driven XML parsing (SAX) with Java+JavaScript

I just wrote SAXScript, an event-driven SAX parser java program invoking some javascript callbacks. It can be used to quickly write a piece of code to parse a huge XML file.

Download


Download saxscript.jar from http://code.google.com/p/lindenb/downloads/list

Invoke


java -jar saxscript.jar (options) [file|url]s

Options


-h (help) this screen
-f read javascript script from file
-e 'script' read javascript script from argument
-D add a variable (as string) in the scripting context.
__FILENAME__ is the current uri.
-n SAX parser is NOT namespace aware (default true)
-v SAX parser is validating (default false)

Callbacks


function startDocument()
{println("Start doc");}
function endDocument()
{println("End doc");}
function startElement(uri,localName,name,atts)
{
print(""+__FILENAME__+" START uri: "+uri+" localName:"+localName);
for(var i=0;atts!=undefined && i< atts.getLength();++i)
{
print(" @"+atts.getQName(i)+"="+atts.getValue(i));
}
println("");
}
function characters(s)
{println("Characters :" +s);}
function endElement(uri,localName,name)
{println("END: uri: "+uri+" localName:"+localName);}

Source Code



Example


The following shell script invokes NCBI/ESearch to retrieve a key to get all the bibliographic references about the Rotaviruses (8793 references).
This key is then used to download each pubmed entry and we then count the number of time each journal (tag is "MedlineTA") was cited.
#!/bin/sh
JAVA=${JAVA_HOME}/bin/java
WEBENV=`${JAVA} -jar saxscript.jar \
-e '
var WebEnv=null;
function startElement(uri,localName,name,atts)
{
if(name=="WebEnv") WebEnv="";
}

function characters(s)
{
if(WebEnv!=null) WebEnv+=s;
}

function endElement(uri,localName,name)
{
if(WebEnv!=null)
{
print(WebEnv);
WebEnv=null;
}
}
' \
"http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&usehistory=y&retmode=xml&term=Rotavirus"`


${JAVA} -jar saxscript.jar -e '
var content=null;
var hash=new Array();
function startElement(uri,localName,name,atts)
{
if(name=="MedlineTA") content="";
}

function characters(s)
{
if(content!=null) content+=s;
}

function endElement(uri,localName,name)
{
if(content!=null)
{
var c=hash[content];
hash[content]=(c==null?1:c+1);
content=null;
}
}
function endDocument()
{
for(var content in hash)
{
println(content+"\t"+ hash[content]);
}
}
' "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&query_key=1&WebEnv=${WEBENV}&retmode=xml" |\
sort -t ' ' -k2n

Result


Acta Gastroenterol Latinoam 1
Acta Histochem Suppl 1
Acta Microbiol Acad Sci Hung 1
Acta Microbiol Hung 1
Acta Microbiol Immunol Hung 1
Acta Pathol Microbiol Scand C 1
Acta Vet Acad Sci Hung 1
Adv Neonatal Care 1
Adv Nurse Pract 1
Adv Ther 1
Adv Vet Med 1
Afr J Med Med Sci 1
Age Ageing 1
AIDS Res Hum Retroviruses 1
AJNR Am J Neuroradiol 1
AJR Am J Roentgenol 1
Akush Ginekol (Sofiia) 1
(...)
Appl Environ Microbiol 87
J Pediatr Gastroenterol Nutr 97
J Virol Methods 130
Lancet 130
Vaccine 158
Pediatr Infect Dis J 177
J Gen Virol 217
Arch Virol 254
J Med Virol 262
J Infect Dis 265
Virology 278
J Virol 460
J Clin Microbiol 514


That's it !
Pierre