27 October 2010

Coloring a black and white drawing with the GIMP: my notebook

I've been kindly asked to draw a proposal for the cover of http://www.cell.com/immunity (It's just a project and I'm not sure it will be accepted). So, here is my notebook for using GIMP (GNU Image Manipulation Program) to colorize a black and white drawing. Foreword:: just like everybody: I do hate gimp.


The Original Drawing


Scan the original drawing at 600dpi. Open it as a ".xcf" color image (the native file format of the Gimp).
(Image stored on OpenWetWare )

Wash Tint


Create a new layer with a white background under the original source (mode = multiply) and fill the shapes . The visible top layer is displayed over the layer to see the strokes.


Gradients


Add a new transparent layer (mode = normal ) and create a radial gradient (transparent-black) from the nucleus to the side of the sheet. Use the eraser to highlight the pyramid-shaped proteins.


Airbrush


Create a new transparent layer (mode = normal ) and use the airbrushes to draw some shadows.




Merge the layers


Et voila !


That's it,

Pierre

24 October 2010

Where are the alternative reading frames in the Human Genome ?

The following post was inspired by a question asked recently on Biostar :"Do exons ever have different reading frames in spliced variants?".

To find those alternative reading frames I've used the table KnownGene available at UCSC from : http://hgdownload.cse.ucsc.edu/goldenPath/hg18/database/knownGene.txt.gz. This file contains the positions of the exons for each transcript in the human genome:

mysql -h genome-mysql.cse.ucsc.edu -A -u genome -D hg18 -e 'select * from knownGene limit 10\G'

(...)
*************************** 7. row ***************************
name: uc009vis.1
chrom: chr1
strand: -
txStart: 4268
txEnd: 6628
cdsStart: 4268
cdsEnd: 4268
exonCount: 4
exonStarts: 4268,4832,5658,6469,
exonEnds: 4692,4901,5805,6628,
proteinID:
alignID: uc009vis.1
*************************** 8. row ***************************
name: uc009vit.1
chrom: chr1
strand: -
txStart: 4268
txEnd: 9622
cdsStart: 4268
cdsEnd: 4268
exonCount: 9
exonStarts: 4268,4832,5658,6469,6720,7095,7777,8130,8775,
exonEnds: 4692,4901,5810,6628,6918,7605,7924,8229,9622,
proteinID:
alignID: uc009vit.1


The following java program creates an array of bytes having a length greater than the length of the human chromosome chr1. This array is initialized with the constant 'NIL'. Then for each chromosome and each transcript, we loop over each exon and we record what was the reading frame (0, 1 or 2) at a given position. If this position was already flagged with another frame, a warning is printed to stdout.

Compilation & Execution


javac BioStar3034.java
java BioStar3034

Result


(...)
chr1:53286155-53286156 (+)
chr1:53286156-53286157 (+)
chr1:53286157-53286158 (+)
chr1:53286158-53286159 (+)
chr1:53286159-53286160 (+)
chr1:53286160-53286161 (+)
chr1:53286161-53286162 (+)
chr1:53286162-53286163 (+)
(...)
java BioStar3034 | sort | uniq | wc -l
300696



(Image from UCSC/OpenWetWare)


That's it
Pierre

11 October 2010

A custom JSP tag for mediawiki

I'm pretty happy with the following custom JSP tag for mediawiki I wrote last week (http://code.google.com/p/code915/source/browse/trunk/charpak/src/WEB-INF/src/fr/inserm/umr915/charpak/j2ee/tags/MediaWikiTag.java). For a given title, it calls the mediawiki api to test if an article exists in a mediawiki installation. In my case, this system will be useful to let my users write some custom annotations to some records from our database.

<u915:mediawiki title="${geneName}" preload="Template:Gene"/>
If the article does exist on the mediawiki installation a simple hyperlink is created. If the article does not exist, an hyperlink for creating/editing the new article is displayed with a predload option: Preloading wikitext presents the user with a partially created page rather than a blank page, possibly with inline instructions for content organization..For example, in the following web application, the blue and red icons point to some existing articles or some articles to be created:


Internals

This custom tag calls the mediawiki api (e.g api.php?action=query&format=xml&titles=Charles+Darwin )and get the XML response as a DOM document. An XPath expression is then evaluated to test if an attribute '@missing' was declared for the given article.
boolean pageFound=false;
String termUTF8 = URLEncoder.encode(term,"UTF-8");
URL url=new URL(baseurl+"/api.php?action=query&format=xml&titles="+termUTF8);
URLConnection con=url.openConnection();
in=con.getInputStream();
Document dom=this.domBuilder.parse(in);
Element e=(Element)this.xpathPage.evaluate(dom,XPathConstants.NODE);
if(e!=null && e.getAttributeNode("missing")==null)
{
pageFound=true;
}
in.close();


That's it,
Pierre

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

10 October 2010

The Geek is dead. Long live the Geek :-)

congratulations to Neil! ;-)



24 September 2010

Connecting to a MongoDB database from R using the C API for MongoDB

Today, Neil posted an article titled" Connecting to a MongoDB database from R using Java". In the current post, I'll show how to use the C API for MongoDB to fetch some MongoDB data from R. The code will be somehow similar to my previous post "A stateful C function for R: parsing Fasta sequences".

OK, First, let's add a few values in mongo:

for(i=1;i> 20;++i)
db.dbsnps.save({_id:"rs"+i,name:"rs"+i});


The C code contains 3 functions.

The first function mongoRconnect connects to the MongoDB server and put the pointer into a R variable.
SEXP mongoRconnect()
{
mongo_connection* conn; /* ptr */
mongo_connection_options opts[1];
mongo_conn_return status;

conn=(mongo_connection*)malloc(sizeof(mongo_connection));
strcpy( opts->host , "127.0.0.1" );
opts->port = 27017;
status = mongo_connect( conn, opts );

return R_MakeExternalPtr(conn, R_NilValue, R_NilValue);
}
The second method mongoRdiconnect closes the connection:
SEXP mongoRdiconnect(SEXP r_handle)
{
mongo_connection* conn;
conn = (mongo_connection*)R_ExternalPtrAddr(r_handle);;
if(conn==NULL) MONGO_ERROR("conn==NULL");
mongo_destroy( conn );
free(conn);
R_ClearExternalPtr(r_handle);
return ScalarInteger(0);
}
The last method mongoRquery scans the database test.dbsnps and inserts the name of the snps into an R array:
SEXP mongoRquery(SEXP r_handle)
{
SEXP values=NULL;
mongo_cursor *cursor;
bson empty[1];
bson_empty( empty );
int i;
mongo_connection* conn = R_ExternalPtrAddr(r_handle);
SEXP* array=NULL;
int array_size=0;
//the R value contains two objects


if(conn==NULL) MONGO_ERROR("handle==NULL");
cursor = mongo_find( conn,
"test.dbsnps",/* ns */
empty,/* fields */
empty,/* return */
0,/* return */
0,/* skip */
0 /* options */
);

while( mongo_cursor_next( cursor ) )
{
bson_iterator it[1];
if ( bson_find( it, &(cursor->current), "name" ))
{
array=(SEXP*)realloc(array,(array_size+1)*sizeof(SEXP));
if(array==NULL) error("out of memory");
array[array_size]=mkChar( bson_iterator_string( it ));
array_size++;
}
}
mongo_cursor_destroy( cursor );
PROTECT(values = allocVector(STRSXP, array_size));
for(i=0;i< array_size;++i)
{
SET_STRING_ELT(values, i, array[i]);
}
free(array);
UNPROTECT(1);
return values;
}
This C code is then be called from R:
mongo <- mongo.open()
mongo.snps(mongo)
mongo.close(mongo)

Result:
[1] "rs1" "rs2" "rs3" "rs4" "rs5" "rs6" "rs7" "rs8" "rs9" "rs10"
[11] "rs11" "rs12" "rs13" "rs14" "rs15" "rs16" "rs17" "rs18" "rs19"



Source code


Makefile

R_HOME=R-2.11.0
MONGO_HOME=mongo-c-driver
run:
gcc -fPIC -I -g -c -Wall -DMONGO_HAVE_STDINT -I ${R_HOME}/include -I ${MONGO_HOME}/src mongoR.c ${MONGO_HOME}/src/*.c
gcc -shared -Wl,-soname,rmongo.so.1 -o librmongo.so *.o
${R_HOME}/bin/R --no-save < mongo.R
clean:
rm *.o


mongoR.c

(again, I'm not sure about those PROTECT/UNPROTECT ...)
#include <ctype.h>
#include <errno.h>
#include <R.h>
#include <Rinternals.h>
#include <bson.h>
#include <mongo.h>

#define MONGO_ERROR(a) { error(a); fputs(a,stdout);exit(EXIT_FAILURE);}

/**
* connect to MONGO
*/
SEXP mongoRconnect()
{
mongo_connection* conn; /* ptr */
mongo_connection_options opts[1];
mongo_conn_return status;

conn=(mongo_connection*)malloc(sizeof(mongo_connection));
if(conn==NULL)
{
MONGO_ERROR("out of memory");
}

strcpy( opts->host , "127.0.0.1" );
opts->port = 27017;

status = mongo_connect( conn, opts );
if(status!= mongo_conn_success)
{
MONGO_ERROR("connection failed");
}

/** the handle is bound a R variable */
return R_MakeExternalPtr(conn, R_NilValue, R_NilValue);
}

/**
* close the mongo connection
*/
SEXP mongoRdiconnect(SEXP r_handle)
{
mongo_connection* conn;
conn = (mongo_connection*)R_ExternalPtrAddr(r_handle);;
if(conn==NULL) MONGO_ERROR("conn==NULL");
mongo_destroy( conn );
free(conn);
R_ClearExternalPtr(r_handle);
return ScalarInteger(0);
}

/**
* get all SNPS
*/
SEXP mongoRquery(SEXP r_handle)
{
SEXP values=NULL;
mongo_cursor *cursor;
bson empty[1];
bson_empty( empty );
int i;
mongo_connection* conn = R_ExternalPtrAddr(r_handle);
SEXP* array=NULL;
int array_size=0;
//the R value contains two objects


if(conn==NULL) MONGO_ERROR("handle==NULL");
cursor = mongo_find( conn,
"test.dbsnps",/* ns */
empty,/* fields */
empty,/* return */
0,/* return */
0,/* skip */
0 /* options */
);

while( mongo_cursor_next( cursor ) )
{
bson_iterator it[1];
if ( bson_find( it, &(cursor->current), "name" ))
{
array=(SEXP*)realloc(array,(array_size+1)*sizeof(SEXP));
if(array==NULL) error("out of memory");
array[array_size]=mkChar( bson_iterator_string( it ));
array_size++;
}
}
mongo_cursor_destroy( cursor );
PROTECT(values = allocVector(STRSXP, array_size));
for(i=0;i< array_size;++i)
{
SET_STRING_ELT(values, i, array[i]);
}
free(array);
UNPROTECT(1);
return values;
}


mongo.R

dyn.load(paste("librmongo", .Platform$dynlib.ext, sep=""))

mongo.open <- function()
{
.Call("mongoRconnect")
}

mongo.close <- function(handler)
{
.Call("mongoRdiconnect",handler)
}

mongo.snps <- function(handler)
{
.Call("mongoRquery",handler)
}

mongo <- mongo.open()
mongo.snps(mongo)
mongo.close(mongo)


That's it,
Pierre

22 September 2010

A Simple tool to get the sex ratio in pubmed.

Just for fun, I wrote a simple java tool to get the sex ratio of the authors in Pubmed. This program fetches a list of names/genders I found in the following perl module: http://cpansearch.perl.org/src/EDALY/Text-GenderFromName-0.33/GenderFromName.pm. The source code is available at

.

(In the following examples, the many names that couldn't be associated to a gender were ignored).

Bioinformatics


Here is the result for "Bioinformatics[journal]"
Women: 3178 (19%) Men: 13149 (80%)
Bioinformatics[Journal]


The 'Lancet' in 2009

Women: 579 (30%) Men: 1331 (69%)
Lancet[Journal] 2009[Date]


Nature in 2009

Women: 1616 (30%) Men: 3768 (69%)
Nature[Journal] 2009[Date]


Nursing in 2009

Women: 29 (70%) Men: 12 (29%)
Nursing[Journal] 2009[Date]



Articles about Charles Darwin

Women: 25 (17%) Men: 118 (82%)
"Darwin C"[PS]



etc... etc..

Source code

/**
* Author:
* Pierre Lindenbaum PhD
* plindenbaum@yahoo.fr
* Source of data:
* http://cpansearch.perl.org/src/EDALY/Text-GenderFromName-0.33/GenderFromName.pm
*/
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.net.URL;
import java.net.URLEncoder;
import java.text.Collator;
import java.util.Locale;
import java.util.Map;
import java.util.TreeMap;
import javax.xml.stream.XMLEventReader;
import javax.xml.stream.XMLInputFactory;
import javax.xml.stream.XMLOutputFactory;
import javax.xml.stream.XMLStreamException;
import javax.xml.stream.XMLStreamWriter;
import javax.xml.stream.events.XMLEvent;

/**
* PubmedGender
*/
public class PubmedGender
{
private Map<String,Float> males=null;
private Map<String,Float> females=null;
private int limit=1000;
private String query="";
private int canvasSize=200;
private boolean ignoreUndefined=false;
private PubmedGender()
{
Collator collator= Collator.getInstance(Locale.US);
collator.setStrength(Collator.PRIMARY);
this.males=new TreeMap<String, Float>(collator);
this.females=new TreeMap<String, Float>(collator);
}

private void loadNames()
throws IOException
{
BufferedReader in=new BufferedReader(new InputStreamReader(new URL("http://cpansearch.perl.org/src/EDALY/Text-GenderFromName-0.33/GenderFromName.pm").openStream()));
String line;
Map<String,Float> map=null;
int posAssign=-1;
while((line=in.readLine())!=null)
{
if(line.startsWith("$Males = {"))
{
map=this.males;
}
else if(line.startsWith("$Females = {"))
{
map=this.females;
}
else if(line.contains("}"))
{
map=null;
}
else if(map!=null && ((posAssign=line.indexOf("=>"))!=-1))
{
String name=line.substring(0,posAssign).replaceAll("'","").toLowerCase().trim();
Float freq=Float.parseFloat(line.substring(posAssign+2).replaceAll("[',]","").toLowerCase().trim());
map.put(name, freq);
}
else
{
map=null;
}
}
in.close();
}
private XMLEventReader newReader(URL url) throws IOException,XMLStreamException
{
XMLInputFactory f= XMLInputFactory.newInstance();
f.setProperty(XMLInputFactory.IS_COALESCING, Boolean.TRUE);
f.setProperty(XMLInputFactory.IS_NAMESPACE_AWARE,Boolean.FALSE);
f.setProperty(XMLInputFactory.IS_REPLACING_ENTITY_REFERENCES,Boolean.TRUE);
f.setProperty(XMLInputFactory.IS_VALIDATING,Boolean.FALSE);
f.setProperty(XMLInputFactory.SUPPORT_DTD,Boolean.FALSE);
XMLEventReader reader=f.createXMLEventReader(url.openStream());
return reader;
}

private void run() throws Exception
{
int countMales=0;
int countFemales=0;
int countUnknown=0;

URL url= new URL(
"http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&term="+
URLEncoder.encode(this.query, "UTF-8")+
"&retstart=0&retmax="+this.limit+"&usehistory=y&retmode=xml&email=plindenbaum_at_yahoo.fr&tool=gender");

XMLEventReader reader= newReader(url);
XMLEvent evt;
String QueryKey=null;
String WebEnv=null;
int countId=0;
while(!(evt=reader.nextEvent()).isEndDocument())
{
if(!evt.isStartElement()) continue;
String tag= evt.asStartElement().getName().getLocalPart();
if(tag.equals("QueryKey"))
{
QueryKey= reader.getElementText().trim();
}
else if(tag.equals("WebEnv"))
{
WebEnv= reader.getElementText().trim();
}
else if(tag.equals("Id"))
{
++countId;
}
}
reader.close();

if(countId!=0)
{
url= new URL("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&WebEnv="+
URLEncoder.encode(WebEnv,"UTF-8")+
"&query_key="+URLEncoder.encode(QueryKey,"UTF-8")+
"&retmode=xml&retmax="+this.limit+"&email=plindenbaum_at_yahoo.fr&tool=mail");

reader= newReader(url);


while(reader.hasNext())
{
evt=reader.nextEvent();
if(!evt.isStartElement()) continue;
if(!evt.asStartElement().getName().getLocalPart().equals("Author")) continue;
String firstName=null;
String initials=null;

while(reader.hasNext())
{
evt=reader.nextEvent();
if(evt.isStartElement())
{
String localName=evt.asStartElement().getName().getLocalPart();
if(localName.equals("ForeName") || localName.equals("FirstName"))
{
firstName=reader.getElementText().toLowerCase();
}
else if(localName.equals("Initials"))
{
initials=reader.getElementText().toLowerCase();
}
}
else if(evt.isEndElement())
{
if(evt.asEndElement().getName().getLocalPart().equals("Author")) break;
}
}
if( firstName==null ) continue;
if( firstName.length()==1 ||
firstName.equals(initials)) continue;

String tokens[]=firstName.split("[ ]+");
firstName="";
for(String s:tokens)
{
if(s.length()> firstName.length())
{
firstName=s;
}
}


if( firstName.length()==1 ||
firstName.equals(initials)) continue;

Float male= this.males.get(firstName);
Float female= this.females.get(firstName);

if(male==null && female==null)
{
//System.err.println("Undefined "+firstName+" / "+lastName);
countUnknown++;
}
else if(male!=null && female==null)
{
countMales++;
}
else if(male==null && female!=null)
{
countFemales++;
}
else if(male < female)
{
countFemales++;
}
else if(female < male)
{
countMales++;
}
else
{
//System.err.println("Undefined "+firstName+" / "+lastName);
countUnknown++;
}
}
reader.close();
}
if(ignoreUndefined) countUnknown=0;

float total= countMales+countFemales+countUnknown;

double radMale=(countMales/total)*Math.PI*2.0;
double radFemale=(countFemales/total)*Math.PI*2.0;
int radius= (canvasSize-2)/2;
String id= "ctx"+System.currentTimeMillis()+""+(int)(Math.random()*1000);
XMLOutputFactory xmlfactory= XMLOutputFactory.newInstance();
XMLStreamWriter w= xmlfactory.createXMLStreamWriter(System.out,"UTF-8");
w.writeStartElement("html");
w.writeStartElement("body");
w.writeStartElement("div");
w.writeAttribute("style","margin:10px;padding:10px;text-align:center;");
w.writeStartElement("div");
w.writeEmptyElement("canvas");
w.writeAttribute("width", String.valueOf(canvasSize+1));
w.writeAttribute("height", String.valueOf(canvasSize+1));
w.writeAttribute("id", id);
w.writeStartElement("script");
w.writeCharacters(
"function paint"+id+"(){var canvas=document.getElementById('"+id+"');"+
"if (!canvas.getContext) return;var c=canvas.getContext('2d');"+
"c.fillStyle='white';c.strokeStyle='black';"+
"c.fillRect(0,0,"+canvasSize+","+canvasSize+");"+
"c.fillStyle='gray';c.beginPath();c.arc("+(canvasSize/2)+","+(canvasSize/2)+","+radius+",0,Math.PI*2,true);c.fill();c.stroke();"+
"c.fillStyle='blue';c.beginPath();c.moveTo("+(canvasSize/2)+","+(canvasSize/2)+");c.arc("+(canvasSize/2)+","+(canvasSize/2)+","+radius+",0,"+radMale+",false);c.closePath();c.fill();c.stroke();"+
"c.fillStyle='pink';c.beginPath();c.moveTo("+(canvasSize/2)+","+(canvasSize/2)+");c.arc("+(canvasSize/2)+","+(canvasSize/2)+","+radius+","+radMale+","+(radMale+radFemale)+",false);c.closePath();c.fill();c.stroke();}"+
"window.addEventListener('load',function(){ paint"+id+"(); },true);"
);
w.writeEndElement();
w.writeEndElement();

w.writeStartElement("span");
w.writeAttribute("style","color:pink;");
w.writeCharacters("Women: "+countFemales+" ("+(int)((countFemales/total)*100.0)+"%)");
w.writeEndElement();
w.writeCharacters(" ");
w.writeStartElement("span");
w.writeAttribute("style","color:blue;");
w.writeCharacters("Men: "+countMales+" ("+(int)((countMales/total)*100.0)+"%)");
w.writeEndElement();
w.writeCharacters(" ");

if(!this.ignoreUndefined)
{
w.writeStartElement("span");
w.writeAttribute("style","color:gray;");
w.writeCharacters("Undefined : "+countUnknown+" ("+(int)((countUnknown/total)*100.0)+"%)");
w.writeEndElement();
}
w.writeEmptyElement("br");

w.writeStartElement("a");
w.writeAttribute("target","_blank");
w.writeAttribute("href","http://www.ncbi.nlm.nih.gov/sites/entrez?db=pubmed&amp;cmd=search&amp;term="+URLEncoder.encode(this.query,"UTF-8"));
w.writeCharacters(this.query);
w.writeEndElement();


w.writeEndElement();
w.writeEndElement();
w.writeEndElement();
w.flush();
w.close();
}

public static void main(String[] args)
{
try
{
PubmedGender app=new PubmedGender();

int optind=0;
while(optind< args.length)
{
if(args[optind].equals("-h") ||
args[optind].equals("-help") ||
args[optind].equals("--help"))
{
System.err.println("Options:");
System.err.println(" -h help; This screen.");
System.err.println(" -w <int> canvas size default:"+app.canvasSize);
System.err.println(" -L <int> limit number default:"+app.limit);
System.err.println(" -i ignore undefined default:"+app.ignoreUndefined);
System.err.println(" query terms...");
return;
}
else if(args[optind].equals("-L"))
{
app.limit=Integer.parseInt(args[++optind]);
}
else if(args[optind].equals("-w"))
{
app.canvasSize=Integer.parseInt(args[++optind]);
}
else if(args[optind].equals("-i"))
{
app.ignoreUndefined=true;
}
else if(args[optind].equals("--"))
{
optind++;
break;
}
else if(args[optind].startsWith("-"))
{
System.err.println("Unknown option "+args[optind]);
return;
}
else
{
break;
}
++optind;
}
if(optind==args.length)
{
System.err.println("Query missing");
return;
}
app.query="";
while(optind< args.length)
{
if(!app.query.isEmpty()) app.query+=" ";
app.query+=args[optind++];
}
app.query=app.query.trim();
if(app.query.trim().isEmpty())
{
System.err.println("Query is empty");
return;
}
app.loadNames();

app.run();

}
catch (Exception e)
{
e.printStackTrace();
}
}
}


That's it

Pierre

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

19 September 2010

Indexing some genomic positions with MongoDB: my benchmark

The aim of this post was to find a good way to index some genomic positions with mongodb. (Update: it was tested on my laptop without replication )

My initial dataset is a list of SNPs on the chromosomes chr22 and chrM from the UCSC.

mysql -N -h genome-mysql.cse.ucsc.edu -A -u genome -D hg18 \
-e 'select chrom,chromStart,name from snp130 where chrom in("chr22","chrM")'


From this dataset, a javascript input for mongodb was generated using the following AWK scripts:
awk -f data2mongo.awk dataset.xls > mongo.js
Each script runs a loop searching some SNPs in a random range on the chr22.
mongo bio mongo.js


The winner: Test 1

Indexing both fields:db.things.ensureIndex({chrom:1,position:1}):
BEGIN {
printf("db.snps.drop();\n");
}

{
printf("db.snps.save({chrom:\"%s\",position:%s,name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("db.snps.ensureIndex({chrom:1,position:1});\n");
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find({chrom:\"chr22\",position:{$gt:pos1,$lt:pos2}});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}

result: seconds:0.228

Test 2

Indexing each field : db.things.ensureIndex({chrom:1}); and db.things.ensureIndex({position:1});.
BEGIN {
printf("db.snps.drop();\n");
}

{
printf("db.snps.save({chrom:\"%s\",position:%s,name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("db.snps.ensureIndex({chrom:1});\n");
printf("db.snps.ensureIndex({position:1});\n");
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find({chrom:\"chr22\",position:{$gt:pos1,$lt:pos2}});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}

Result: seconds:0.25



Test 3


No index.


BEGIN {
printf("db.snps.drop();\n");
}

{
printf("db.snps.save({chrom:\"%s\",position:%s,name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find({chrom:\"chr22\",position:{$gt:pos1,$lt:pos2}});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}

Result:seconds:277.751



Test 4

String padding:the chromosome and the position are concatenated in a fixed-length string


BEGIN {
printf("db.snps.drop();\n");
printf("function pad2(s,L) { while(s.length<L) { s=\"0\"+s;} return s;}\n");
printf("function pad(chrom,position) { return pad2(chrom,2)+\":\"+pad2(\"\"+position,10);}\n");
}

{
gsub(/chr/,"",$1);
printf("db.snps.save({position:pad(\"%s\",%s),name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find({position:{$gt:pad(\"22\",pos1),$lt:pad(\"22\",pos2)}});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}
Result: seconds:169.028

Test 5

string padding + index


BEGIN {
printf("db.snps.drop();\n");
printf("function pad2(s,L) { while(s.length<L) { s=\"0\"+s;} return s;}\n");
printf("function pad(chrom,position) { return pad2(chrom,2)+\":\"+pad2(\"\"+position,10);}\n");
}

{
gsub(/chr/,"",$1);
printf("db.snps.save({position:pad(\"%s\",%s),name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("db.snps.ensureIndex({position:1});\n");
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find({position:{$gt:pad(\"22\",pos1),$lt:pad(\"22\",pos2)}});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}

Result: seconds:0.292



Test 6

_id as a padded string as _id

.
BEGIN {
printf("db.snps.drop();\n");
printf("function pad2(s,L) { while(s.length<L) { s=\"0\"+s;} return s;}\n");
printf("function pad(chrom,position) { return pad2(chrom,2)+\":\"+pad2(\"\"+position,10);}\n");
}

{
gsub(/chr/,"",$1);
printf("db.snps.save({_id:pad(\"%s\",%s),name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find({_id:{$gt:pad(\"22\",pos1),$lt:pad(\"22\",pos2)}});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}

result : seconds:1.252



Test 7

Using a padded string for _id and min()/max() for searching
BEGIN {
printf("db.snps.drop();\n");
printf("function pad2(s,L) { while(s.length<L) { s=\"0\"+s;} return s;}\n");
printf("function pad(chrom,position) { return pad2(chrom,2)+\":\"+pad2(\"\"+position,10);}\n");
}

{
gsub(/chr/,"",$1);
printf("db.snps.save({_id:pad(\"%s\",%s),name:\"%s\"});\n",$1,$2,$3);
}

END {
printf("var tStart=new Date().getTime();\n");
printf("for(i=0;i< 1000;++i)\n{\n");
printf("var pos1 =Math.floor(Math.random()*50000000);\n");
printf("var pos2 =pos1 + Math.floor(Math.random()*10000);\n");
printf("var c=db.snps.find().min({_id:pad(\"22\",pos1)}).max({_id:pad(\"22\",pos2)});\n");
printf("while(c.hasNext()) c.next();\n");
printf("}\n");
printf("print(\"seconds:\"+(new Date().getTime()-tStart)/1000)");
}

result Seconds:3.189



Test 8

I also used a composite _id: db.snps.save({_id:{chrom:"chr22",position:14430966},name:"rs2844899"}); but i was not able to query this table with $gt/$lt:


> db.snps.find({_id:{chrom:"chr22",position:{$gt:14430966}}}).count()
0
> db.snps.find({_id:{chrom:"chr22",position:14430966}}).count()
1


Any other idea for indexing those data ? Feel free to leave a message here or on biostar or stackoverflow.


That's it

Pierre