13 April 2012

I get mail: result of the contest "Draw a PhD Candidate"

I recently submitted the following illustration for a contest organized by the French University of Nice : "Draw a PHD candidate".



'Draw a PhD' Contest by ~yokofakun on deviantART



and I received the following email today:

 
De : D. Rauch
Envoyé le : Jeudi 12 avril 2012 21h31
Objet : résultat concours "Dessine-moi un doctorant !"

Monsieur,

L’AJC 06  a l’honneur de vous informer que votre œuvre « Gwennan André,
doctorante » a été retenue par le jury du concours « Dessine-moi un doctorant !
» parmi les œuvres gagnantes.

Votre œuvre ayant été classée à la première place dans la catégorie "public",
vous avez gagné un Ipad 2.
(...)

(translation:)

Sir,

The association AJC 06 is proud to inform you that your work titled "Gwennan André,
PhD student " was chosen by the jury of the contest" Draw me a PhD!
"Among the winning entries.

Your work has been ranked in first place in the "public" section, you have won an IPad 2.
 

Unfortunately I cannot go to Nice next week !!!

I have to renounce and say "goodbye" to this IPAD2 :-((


Sniffffffff


Pierre


EDIT: Simple & stupid Google translation of the text:

"Zig-zag path to understanding myself that it was research health that interested me, because this is an area overlooked by leaving the baccalaureate (At least when the secondary school is far from major cities, large research centers and our families are not the domain!). So the first year of medical school. This is where I found my way, I loved to understand the disease but did not want to be in direct contact with the patient. So reorientation IUT biological engineering: essentially the stage "business" of the second year. Two months in a lab search for "Robert Gordon University" in Aberdeen, Scotland history to be quickly immersed in a research lab and into human health in English. Excellent and successful experience: motivation professional found. Then L3 cell biology, physiology and genetics in Rennes and Master of Nantes "Biology Health" search option with a therapeutic internship at the Institute of alternation in the Thorax team Gervaise Loirand.
In 2010, got a scholarship to complete my thesis Ministerial in this team on a project of Vincent Sauzeau. Our goal is to understand the mechanisms underlying vascular disease to better manage diseases often heavy consequences.
Projects after the thesis: do a postdoc abroad to live fully international aspect of the research and see Search elsewhere. Trying to reconcile this with the personal life, something not always evident with extension to these studies, especially for a woman. Then return to do research in France, ideally in a position of Research Fellow in the audience. But it's still very open in my head (private, public, engineer ... to see according to experience and opportinities in 5 years ...!)
"

06 April 2012

Apache SOLR and GeneOntology: Creating the JQUERY-UI client (with autocompletion)

In the previous post I showed how to use Apache SOLR to index the content of gene ontology. In this post I will show how the SOLR server can be used as a source of data for a "JQuery-UI autocomplete" HTML input.

The JQuery autocomplete Widget expects an JSON array of Objects with label and value properties :
[ { label: "Choice1", value: "value1" }, ... ]

We need to tell SOLR to return the GO terms formatted this way.

Fortunately SOLR provides a mechanism to captures the output of the normal XML response and applies a XSLT transformation to it.

The following XSLT stylesheet transforms the XML response to JSON, it is created in the conf/xslt directory (apache-solr-3.5.0/example/solr/conf/xslt/solr2json.xsl)


We can test the stylesheet by adding the tr parameter in the query URL:
$ curl "http://localhost:8983/solr/select?q=go_definition:cancer&wt=xslt&tr=solr2json.xsl"

[
    {
        "label": "Z-phenylacetaldoxime metabolic process",
        "value": "GO:0018983"
    },
    {
        "label": "epothilone metabolic process",
        "value": "GO:0050813"
    },
    {
        "label": "epothilone biosynthetic process",
        "value": "GO:0050814"
    },
    {
        "label": "aflatoxin biosynthetic process",
        "value": "GO:0045122"
    },
    {
        "label": "aflatoxin metabolic process",
        "value": "GO:0046222"
    },
    {
        "label": "aflatoxin catabolic process",
        "value": "GO:0046223"
    }
]

The jquery client

Let's create the HTML/JQuery Client, I have downloaded the component for jquery and jquery-ui in the apache-solr-3.5.0/example/work/Jetty_xxxx_solr.war__solr__xxxx/webapp/ (and that's a bad idea, everything will be deleted with the SOLR server will be stopped or re-deployed).

I've also created the following HTML page 'solr.html' in the same directory: Let's visit http://localhost:8983/solr/solr.html to test the page:


(image via openwetware)

 That worked ! :-)

That's it,

Pierre

Indexing the content of Gene Ontology with apache SOLR

Via Wikipedia:"Solr (http://lucene.apache.org/solr/) is an open source enterprise search platform from the Apache Lucene project. Its major features include powerful full-text search, hit highlighting, faceted search, dynamic clustering, database integration, and rich document (e.g., Word, PDF) handling. Providing distributed search and index replication, Solr is highly scalable." In the this post, I'll show how I've used SOLR to index the content of GeneOntology.

Download and install SOLR

Download from http://mirrors.ircam.fr/pub/apache/lucene/solr/3.5.0/apache-solr-3.5.0.tgz.
tar xvfz apache-solr-3.5.0.tgz
rm apache-solr-3.5.0.tgz

Configure schema.xml

We need to tell SOLR about the which fields of GO will be indexed, what are their type, how they should be tokenized and parsed. This information is defined in the schema.xml. The following components will be indexed: accession, name, synonym and definition. Edit apache-solr-3.5.0/example/solr/conf/schema.xml and add the following <fields>:

<field name="go_name" type="text_general" indexed="true" stored="true" multiValued="false"/>
<field name="go_synonym" type="text_general" indexed="true" stored="true" multiValued="true"/>
<field name="go_definition" type="text_general" indexed="true" stored="true" multiValued="false"/>

Start the SOLR server

In this example, the SOLR server is started using the simple Jetty server provided in the distribution:

$ cd apache-solr-3.5.0/example/example
$ java -jar start.jar

(...)

Indexing Gene Ontology

Go is downloaded as RDF/XML from http://archive.geneontology.org/latest-termdb/go_daily-termdb.rdf-xml.gz
 
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE go:go PUBLIC "-//Gene Ontology//Custom XML/RDF Version 2.0//EN" "http://www.geneontology.org/dtd/go.dtd">

<go:go xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
    <rdf:RDF>
        <go:term rdf:about="http://www.geneontology.org/go#GO:0000001">
            <go:accession>GO:0000001</go:accession>
            <go:name>mitochondrion inheritance</go:name>
            <go:synonym>mitochondrial inheritance</go:synonym>
            <go:definition>The distribution of mitochondria, including the mitochondrial genome, into daughter cells after mitosis or meiosis, mediated by interactions between mitochondria and the cytoskeleton.</go:definition>
            <go:is_a rdf:resource="http://www.geneontology.org/go#GO:0048308" />
            <go:is_a rdf:resource="http://www.geneontology.org/go#GO:0048311" />
        </go:term>
        <go:term rdf:about="http://www.geneontology.org/go#GO:0000002">
            <go:accession>GO:0000002</go:accession>
            <go:name>mitochondrial genome maintenance</go:name>
            <go:definition>The maintenance of the structure and integrity of the mitochondrial genome; includes replication and segregation of the mitochondrial chromosome.</go:definition>
            <go:is_a rdf:resource="http://www.geneontology.org/go#GO:0007005" />
            <go:dbxref rdf:parseType="Resource">
                <go:database_symbol>InterPro</go:database_symbol>
(...)
 
We now need to transform this XML file to another XML file that can be indexed by the SOLR server.  

"You can modify a Solr index by POSTing XML Documents containing instructions to add (or update) documents, delete documents, commit pending adds and deletes, and optimize your index."

The following XSLT stylesheet is used to transform the RDF/XML for GO:


$ xsltproc --novalid go2solr.xsl go_daily-termdb.rdf-xml.gz > add.xml
$ cat add.xml

Before indexing the current disk usage under apache-solr-3.5.0 is 136Mo. We can now use the java utiliy post.jar to index GeneOntology.

 $ cd  ~/package/apache-solr-3.5.0/example/exampledocs
 $ java -jar post.jar  add.xml

SimplePostTool: version 1.4
SimplePostTool: POSTing files to http://localhost:8983/solr/update..
SimplePostTool: POSTing file jeter.xml
SimplePostTool: COMMITting Solr index changes..

After indexing, the disk usage under apache-solr-3.5.0 is 153Mo.

Querying

Search for the GO terms having go:definition containing "cancer" a go:name containing "genome" but discard those having go:definition containing "metabolism".
 curl "http://localhost:8983/solr/select/?q=go_definition%3Acancer+go_name%3Agenome+-go_definition%3Ametabolism&version=2.2&start=0&rows=10&indent=on"
Same query, but return the result as a JSON structure:
 curl "http://localhost:8983/solr/select/?q=go_definition%3Acancer+go_name%3Agenome+-go_definition%3Ametabolism&version=2.2&start=0&rows=10&indent=on&wt=json"
That's it, Pierre

08 March 2012

My first walker for the GATK : my notebook

This is my first notebook for developping a new Walker for the Genome Analysis Toolkit. This post was mostly inspired by the following pdf: kvg_20_line_lifesavers_mad_v2.pptx.pdf.

Get the sources

git clone http://github.com/broadgsa/gatk.git GATK.dev
the javac compiler also requires the following library from google :http://code.google.com/p/cofoja/.

A first "Short-Reads" walker

The following class ReadWalker scans the reads and print them as fasta. The @Output annotation tells the GATK that we're going to channel our output through the java.io.PrintStream object. This field is automatically filled by the application runtime.

Compilation

javac -cp /path/to/GenomeAnalysisTK.jar:/path/to/cofoja-1.0-r139.jar:. \
 -sourcepath src \
 -d tmp src/mygatk/HelloRead.java
jar cvf HelloRead.jar -C tmp .

Running

Here I'm using a BAM from the 'examples' folder of samtools. (We need to pre-process this BAM with picard AddOrReplaceReadGroups). We then use our library as follow:
java -cp path/to/GenomeAnalysisTK.jar:HelloRead.jar \
org.broadinstitute.sting.gatk.CommandLineGATK -T HelloRead \
 -I test.bam \
 -R ${SAMTOOLS}/examples/ex1.fa 

Result:

The Makefile

That's it, Pierre

04 March 2012

Java Remote Method Invocation (RMI) for Bioinformatics

"Java Remote Method Invocation (Java RMI) enables the programmer to create distributed Java technology-based to Java technology-based applications, in which the methods of remote Java objects can be invoked from other Java virtual machines*, possibly on different hosts. "[Oracle] In the current post a java client will send a java class to the server that will analyze a DNA sequence fetched from the NCBI, using the RMI technology.

Files and directories

I In this example, my files are structured as defined below:
./sandbox/client/FirstBases.java
./sandbox/client/GCPercent.java
./sandbox/client/SequenceAnalyzerClient.java
./sandbox/server/SequenceAnalyzerServiceImpl.java
./sandbox/shared/SequenceAnalyzerService.java
./sandbox/shared/SequenceAnalyzer.java
./client.policy
./server.policy

The Service: SequenceAnalyzerService.java

The remote service provided by the server is defined as an interface named SequenceAnalyzerService: it fetches a DNA sequence for a given NCBI-gi, processes the sequence with an instance of SequenceAnalyzer (see below) and returns a serializable value (that is to say, we can transmit this value through the network).

Extract a value from a DNA sequence : SequenceAnalyzer

The interface SequenceAnalyzer defines how the remote service should parse a sequence. A SAX Parser will be used by the 'SequenceAnalyzerService' to process a TinySeq-XML document from the NCBI. The method characters is called each time a chunck of sequence is found. At the end, the remote server will return the value calculated from getResult:

Server side : an implementation of SequenceAnalyzerService

The class SequenceAnalyzerServiceImpl is an implementation of the service SequenceAnalyzerService. In the method analyse, a SAXParser is created and the given 'gi' sequence is downloaded from the NCBI. The instance of SequenceAnalyzer received from the client is invoked for each chunck of DNA. At the end, the "value" calculated by the instance of SequenceAnalyzer is returned to the client through the network. The 'main' method contains the code to bind this service to the RMI registry:

Client side

On the client side, we're going to connect to the SequenceAnalyzerService and send two distinct implementations of SequenceAnalyzer. What's interesting here: the server doesn't know anything about those implementations of SequenceAnalyzer. The client's java compiled classes have to be sent to the service.

GCPercent.java

A first implementation of 'SequenceAnalyzer' computes the GC% of a sequence:

FirstBases

The second implementation of 'SequenceAnalyzer' retrieves the first bases of a sequence.

The Client

And here is the java code for the client. The client connects to the RMI server and invokes 'analyse' with the two instances of SequenceAnalyzer for some NCBI-gi:

A note about security

As the server/client doesn't want to receive some malicious code, we have to use some policy files:
server.policy:

client.policy:

Compiling and Running

Compiling the client

javac -cp . sandbox/client/SequenceAnalyzerClient.java

Compiling the server

javac -cp . sandbox/server/SequenceAnalyzerServiceImpl.java

Starting the RMI registry

${JAVA_HOME}/bin/rmiregistry

Starting the SequenceAnalyzerServiceImpl

$ java \
 -Djava.security.policy=server.policy \
 -Djava.rmserver.codebase=file:///path/to/RMI/ \
 -cp . sandbox.server.SequenceAnalyzerServiceImpl

SequenceAnalyzerService bound.

Running the client

$ java  \
 -Djava.rmi.server.codebase=file:///path/to/RMI/ \
 -Djava.security.policy=client.policy  \
 -cp . sandbox.client.SequenceAnalyzerClient  localhost

gi=25 gc%=2.1530612244897958
gi=25 start=TAGTTATTC
gi=26 gc%=2.1443298969072164
gi=26 start=TAGTTATTAA
gi=27 gc%=2.3022222222222224
gi=27 start=AACCAGTATTA
gi=28 gc%=2.376543209876543
gi=28 start=TCGTA
gi=29 gc%=2.2014742014742015
gi=29 start=TCTTTG
That's it, Pierre

31 January 2012

Inside the Variation Toolkit: Tools for Gene Ontology

GeneOntologyDbManager is a C++ tool that is part of my experimental Variation Toolkit.
This program is a set of tools for GeneOntology, it is based on the sqlite3 library.

Download

Download the sources from Google-Code using subversion:....
svn checkout http://variationtoolkit.googlecode.com/svn/trunk/ variationtoolkit-read-only
... or update the sources of an existing installation...
cd variationtoolkit
svn update
... and edit the variationtoolkit/congig.mk file.

Dependencies

Compilation

Define "SQLITE_LIB" and "SQLITE_CFLAGS" in config.mk (see HowToInstall )
$ cd variationtoolkit/src/
$ make ../bin/godbmgr 

if ! [ -z "-lsqlite3" ] ;then g++ -o ../bin/godbmgr godatabasemgr.cpp xsqlite.cpp application.o xstream.o xxml.o -g `xml2-config --cflags `  /usr/include/sqlite3.h   -lz -lsqlite3 `xml2-config  --libs` ; else g++ -o ../bin/godbmgr godatabasemgr.cpp  -DNOSQLITE -O3 -Wall  ; fi

Usage

godbmgr (program-name) -f database.sqlite [options] (file1.vcf file2... | stdin )

Sub-Program: loadrdf

Loads the RDF/XML GO database (http://archive.geneontology.org/latest-termdb/go_daily-termdb.rdf-xml.gz) into the sqlite3 database.

Usage

godbmgr loadrdf -f database.sqlite (stdin|file)

Options

  • -f (filename) the sqlite3 database

Example

$ curl -s "http://archive.geneontology.org/latest-termdb/go_daily-termdb.rdf-xml.gz" |\
  gunzip -c |\
  godbmgr loadrdf -f database.sqlite
list the content of the database:
$ sqlite3 -separator '  ' -header  database.sqlite 'select * from TERM where acn="GO:0000007"'
acn xml
GO:0000007 <go:term xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" rdf:about="http://www.geneontology.org/go#GO:0000007">
            <go:accession xmlns:go="http://www.geneontology.org/dtds/go.dtd#">GO:0000007</go:accession>
            <go:name xmlns:go="http://www.geneontology.org/dtds/go.dtd#">low-affinity zinc ion transmembrane transporter activity</go:name>
            <go:definition xmlns:go="http://www.geneontology.org/dtds/go.dtd#">Catalysis of the transfer of a solute or solutes from one side of a membrane to the other according to the reaction: Zn2+ = Zn2+, probably powered by proton motive force. In low affinity transport the transporter is able to bind the solute only if it is present at very high concentrations.</go:definition>
            <go:is_a xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" rdf:resource="http://www.geneontology.org/go#GO:0005385"/>
        </go:term>


$ sqlite3 -separator '  ' -header  database.sqlite 'select * from TERM2REL where acn="GO:0000007"'
acn rel target
GO:0000007 is_a GO:0005385

Sub-Program: loadgoa

inserts the database for GOA into a sqlite3 database (e.g: ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/HUMAN/gene_association.goa_human.gz)

Usage

godbmgr loadgoa -f database.sqlite (stdin|file)

Options

  • -f (filename) the sqlite3 database

Examples

$  curl -s "ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/HUMAN/gene_association.goa_human.gz" |\
     gunzip -c |\
     godbmgr loadgoa -f database.sqlite
list the content of the database:
$ sqlite3 -line   database.sqlite 'select * from GOA where term="GO:0005385" limit 2' 
              DB = UniProtKB
    DB_Object_ID = B3KU87
DB_Object_Symbol = SLC30A6
            term = GO:0005385
  DB_Object_Name = cDNA FLJ45816 fis, clone NT2RP7019682, highly similar to Homo sapiens solute carrier family 30 (zinc transporter), member 6 (SLC30A6), mRNA
         Synonym = B3KU87_HUMAN|SLC30A6|hCG_23082|IPI01009565|B7WP49
  DB_Object_Type = protein

              DB = UniProtKB
    DB_Object_ID = B5MCR8
DB_Object_Symbol = SLC30A6
            term = GO:0005385
  DB_Object_Name = Solute carrier family 30 (Zinc transporter), member 6, isoform CRA_b
         Synonym = B5MCR8_HUMAN|SLC30A6|hCG_23082|IPI00894292
  DB_Object_Type = protein

Sub-Program: desc

print the descendants (children) of a given GO node.

Usage

godbmgr desc -f db.sqlite [options] term1 term2 ... termn

Options

Examples

Default output

$ godbmgr desc -f database.sqlite "GO:0005385"
GO:0000006
GO:0000007
GO:0005385
GO:0015341
GO:0015633
GO:0016463
GO:0022883

xml/rdf output

$ godbmgr desc -f database.sqlite  -t xml "GO:0005385" | head

<go:go xmlns:go='http://www.geneontology.org/dtds/go.dtd#' xmlns:rdf='http://www.w3.org/1999/02/22-rdf-syntax-ns#'>
 <rdf:RDF>
<go:term xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" rdf:about="http://www.geneontology.org/go#GO:0000006">
            <go:accession xmlns:go="http://www.geneontology.org/dtds/go.dtd#">GO:0000006</go:accession>
            <go:name xmlns:go="http://www.geneontology.org/dtds/go.dtd#">high affinity zinc uptake transmembrane transporter activity</go:name>
            <go:definition xmlns:go="http://www.geneontology.org/dtds/go.dtd#">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.</go:definition>
            <go:is_a xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" rdf:resource="http://www.geneontology.org/go#GO:0005385"/>
        </go:term>
<go:term xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" rdf:about="http://www.geneontology.org/go#GO:0000007">
            <go:accession xmlns:go="http://www.geneontology.org/dtds/go.dtd#">GO:0000007</go:accession>

GOA output

$godbmgr desc -f database.sqlite -t goa  "GO:0005385"

UniProtKB B3KU87 SLC30A6 GO:0005385 cDNA FLJ45816 fis, clone NT2RP7019682, highly similar to Homo sapiens solute carrier family 30 (zinc transporter), member 6 (SLC30A6), mRNA B3KU87_HUMAN|SLC30A6|hCG_23082|IPI01009565|B7WP49 protein
UniProtKB B5MCR8 SLC30A6 GO:0005385 Solute carrier family 30 (Zinc transporter), member 6, isoform CRA_b B5MCR8_HUMAN|SLC30A6|hCG_23082|IPI00894292 protein
(..)
UniProtKB Q99726 SLC30A3 GO:0015633 Zinc transporter 3 ZNT3_HUMAN|ZNT3|SLC30A3|IPI00293793|Q8TC03protein

TSV output

$ godbmgr desc -f database.sqlite  -t tsv "GO:0022857" |\
    cut -c 1-100 |\
    head
#go:accession go.name go.def
GO:0000006 high affinity zinc uptake transmembrane transporter activity Catalysis of the transfer of
GO:0000007 low-affinity zinc ion transmembrane transporter activity Catalysis of the transfer of a s
GO:0000064 L-ornithine transmembrane transporter activity Catalysis of the transfer of L-ornithine f
GO:0000095 S-adenosylmethionine transmembrane transporter activity Catalysis of the transfer of S-ad
GO:0000099 sulfur amino acid transmembrane transporter activity Catalysis of the transfer of sulfur 
GO:0000100 S-methylmethionine transmembrane transporter activity Catalysis of the transfer of S-meth
GO:0000102 L-methionine secondary active transmembrane transporter activity Catalysis of the transfe
GO:0000227 oxaloacetate secondary active transmembrane transporter activity Catalysis of the transfe
GO:0000269 toxin export channel activity Enables the energy independent passage of toxins, sized les
(...)

Sub-Program: asc

prints the ascendants (parents) of a given node.

Usage

godbmgr asc -f db.sqlite [options] term1 term2 ... termn

Options

Examples

Default output

$ godbmgr asc -f database.sqlite "GO:0022857"
GO:0003674
GO:0005215
GO:0022857
all

xml/rdf output

$ godbmgr asc -f database.sqlite  -t xml "GO:0022857" | head

<go:go xmlns:go='http://www.geneontology.org/dtds/go.dtd#' xmlns:rdf='http://www.w3.org/1999/02/22-rdf-syntax-ns#'>
 <rdf:RDF>
<go:term xmlns:go="http://www.geneontology.org/dtds/go.dtd#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" rdf:about="http://www.geneontology.org/go#GO:0003674">
            <go:accession xmlns:go="http://www.geneontology.org/dtds/go.dtd#">GO:0003674</go:accession>
            <go:name xmlns:go="http://www.geneontology.org/dtds/go.dtd#">molecular_function</go:name>
            <go:synonym xmlns:go="http://www.geneontology.org/dtds/go.dtd#">GO:0005554</go:synonym>
            <go:synonym xmlns:go="http://www.geneontology.org/dtds/go.dtd#">molecular function</go:synonym>
            <go:synonym xmlns:go="http://www.geneontology.org/dtds/go.dtd#">molecular function unknown</go:synonym>
            <go:definition xmlns:go="http://www.geneontology.org/dtds/go.dtd#">Elemental activities, such as catalysis or binding, describing the actions of a gene product at the molecular level. A given gene product may exhibit one or more molecular functions.</go:definition>
            <go:comment xmlns:go="http://www.geneontology.org/dtds/go.dtd#">Note that, in addition to forming the root of the molecular function ontology, this term is recommended for use for the annotation of gene products whose molecular function is unknown. Note that when this term is used for annotation, it indicates that no information was available about the molecular function of the gene product annotated as of the date the annotation was made; the evidence code ND, no data, is used to indicate this.</go:comment>

GOA output

$godbmgr desc -f database.sqlite -t goa  "GO:0005385"

UniProtKB B3KU87 SLC30A6 GO:0005385 cDNA FLJ45816 fis, clone NT2RP7019682, highly similar to Homo sapiens solute carrier family 30 (zinc transporter), member 6 (SLC30A6), mRNA B3KU87_HUMAN|SLC30A6|hCG_23082|IPI01009565|B7WP49 protein
UniProtKB B5MCR8 SLC30A6 GO:0005385 Solute carrier family 30 (Zinc transporter), member 6, isoform CRA_b B5MCR8_HUMAN|SLC30A6|hCG_23082|IPI00894292 protein
(..)
UniProtKB Q99726 SLC30A3 GO:0015633 Zinc transporter 3 ZNT3_HUMAN|ZNT3|SLC30A3|IPI00293793|Q8TC03protein

TSV output

$ godbmgr asc -f database.sqlite  -t tsv "GO:0022857"   
#go:accession go.name go.def
GO:0003674 molecular_function Elemental activities, such as catalysis or binding, describing the actions of a gene product at the molecular level. A given gene product may exhibit one or more molecular functions.
GO:0005215 transporter activity Enables the directed movement of substances (such as macromolecules, small molecules, ions) into, out of or within a cell, or between cells.
GO:0022857 transmembrane transporter activity Enables the transfer of a substance from one side of a membrane to the other.
all all .

Sub-program: goa

Annotate a TSV file with the GOA annotation.

Usage

godbmgr goa -f db.sqlite [options] (stdin|files)

Options

  • -f (filename) the sqlite3 database
  • -c (column index) REQUIRED. The observed column.

Example

$ echo -e "#MyGene\nHello\nNOTCH2" |\
   godbmgr goa -c 1 -f database.sqlite  |\
   head -n 4 |\
   verticalize

>>> 2
$1 #MyGene                Hello
$2 DB                     .
$3 DB_Object_ID           .
$4 DB_Object_Symbol       .
$5 term                   .
$6 DB_Object_Name,Synonym .
$7 DB_Object_Type         .
$8   ???                    .
<<< 2

>>> 3
$1 #MyGene                NOTCH2
$2 DB                     UniProtKB
$3 DB_Object_ID           Q04721
$4 DB_Object_Symbol       NOTCH2
$5 term                   GO:0001709
$6 DB_Object_Name,Synonym Neurogenic locus notch homolog protein 2
$7 DB_Object_Type         NOTC2_HUMAN|NOTCH2|IPI00297655|Q5T3X7|Q99734|Q9H240
$8   ???                    protein
<<< 3

>>> 4
$1 #MyGene                NOTCH2
$2 DB                     UniProtKB
$3 DB_Object_ID           Q04721
$4 DB_Object_Symbol       NOTCH2
$5 term                   GO:0004872
$6 DB_Object_Name,Synonym Neurogenic locus notch homolog protein 2
$7 DB_Object_Type         NOTC2_HUMAN|NOTCH2|IPI00297655|Q5T3X7|Q99734|Q9H240
$8   ???                    protein

Sub-Program: grep

filters the line having an identifier (gene...) that is a children of a given GO term.

Usage

godbmgr grep -f db.sqlite [options] (stdin|files)

Options

  • -f (filename) the sqlite3 database
  • -c (column index) REQUIRED. The column containing a GO:acn
  • -v inverse selection
  • -t (GO:acn) add a GO term in the filter (One is REQUIRED).
  • -r (rel) add a go relationship (http://www.obofoundry.org/ro/) (OPTIONAL, default: it adds "is_a").

Example

$ 
$ echo -e "#MyACN\nGO:0003674\nGO:0001618\n" |\
  godbmgr grep -f database.sqlite -c 1 -t GO:0004872 -t GO:0004879 

#MyACN
GO:0001618
$ echo -e "#MyACN\nGO:0003674\nGO:0001618\n" |\
  godbmgr grep -f database.sqlite -c 1 -t GO:0004872 -t GO:0004879 -v

#MyACN
GO:0003674


That's it,

Pierre

28 January 2012

Inside the variation toolkit: VCF2XML

vcf2xml is C++ tool that is part of my Variation Toolkit.
It transforms a "Variant Call Format document" to XML, so it can be later processed with xslt, xquery, etc...

Dependencies

Download

Download the sources from Google-Code using subversion:....
svn checkout http://variationtoolkit.googlecode.com/svn/trunk/ variationtoolkit-read-only
... or update the sources of an existing installation...
cd variationtoolkit
svn update
... and edit the variationtoolkit/congig.mk file.

Compiling:

$ cd variationtoolkit/src/
$ make ../bin/vcf2xml

g++ -o ../bin/vcf2xml vcf2xml.cpp application.o -O3 -Wall `xml2-config --cflags --libs` -lz

Usage:

vcf2xml (file.vcf | stdin)

Example:

$ vcf2xml input.vcf | xmllint --format -

<?xml version="1.0" encoding="UTF-8"?>
<vcf>
  <head>
    <meta key="fileformat">VCFv4.1</meta>
    <meta key="samtoolsVersion">0.1.17 (r973:277)</meta>
    <infos>
      <info>
        <id>DP</id>
        <number>1</number>
        <type>Integer</type>
        <description>Raw read depth</description>
      </info>
      <info>
        <id>DP4</id>
        <number>4</number>
        <type>Integer</type>
        <description># high-quality ref-forward bases</description>
      </info>
      <info>
        <id>MQ</id>
(...)
      </calls>
    </variation>
    <variation>
      <chrom>chr1</chrom>
      <pos>112697</pos>
      <ref>T</ref>
      <alt>G</alt>
      <qual>10.4</qual>
      <infos>
        <info key="DP">1</info>
        <info key="AF1">1</info>
        <info key="AC1">2</info>
        <info key="DP4">0,0,0,1</info>
        <info key="MQ">60</info>
        <info key="FQ">-30</info>
      </infos>
      <calls>
        <call sample="input.bam">
          <prop key="GT">1/1</prop>
          <prop key="PL">40,3,0</prop>
          <prop key="GQ">5</prop>
        </call>
      </calls>
    </variation>
  </body>
</vcf>
That's it,
Pierre

Insert your VCFs in a sqlite database.

vcf2sqlite is C++ tool that is part of my Variation Toolkit.
It inserts a "Variant Call Format document" (VCF) into a sqlite3 database.

Download

Download the sources from Google-Code using subversion:....
svn checkout http://variationtoolkit.googlecode.com/svn/trunk/ variationtoolkit-read-only
... or update the sources of an existing installation...
cd variationtoolkit
svn update
... and edit the variationtoolkit/congig.mk file.

Dependencies

http://www.sqlite.org/ : libraries and headers for sqlite3.

Compilation

Define "SQLITE_LIB" and "SQLITE_CFLAGS" in config.mk (see HowToInstall )
$ cd variationtoolkit/src/
$ make ../bin/vcf2sqlite 

if ! [ -z "$(SQLITE_LIB)" ] ;then g++ -o ../bin/vcf2sqlite vcf2sqlite.cpp xsqlite.cpp application.o -O3 -Wall -lz   ; else g++ -o ../bin/vcf2sqlite vcf2sqlite.cpp  -DNOSQLITE -O3 -Wall  ; fi

Usage

vcf2sqlite -f database.sqlite (file1.vcf file2... | stdin )

Options

  • -f (file) sqlite3 database (REQUIRED).

Schema


Example:

$ vcf2sqlite -f db.sqlite file.vcf
$ sqlite3 -line db.sqlite  "select * from VCFCALL LIMIT 4"

       id = 1
   nIndex = 0
vcfrow_id = 1
sample_id = 1
     prop = GT
    value = 1/1

       id = 2
   nIndex = 1
vcfrow_id = 1
sample_id = 1
     prop = PL
    value = 46,6,0

       id = 3
   nIndex = 2
vcfrow_id = 1
sample_id = 1
     prop = GQ
    value = 10

       id = 4
   nIndex = 0
vcfrow_id = 2
sample_id = 1
     prop = GT
    value = 1/1

$ sqlite3 -column -header  db.sqlite \
   "select SAMPLE.name,VCFCALL.value,count(*) from VCFCALL,SAMPLE where SAMPLE.id=VCFCALL.sample_id and prop='GT' group by SAMPLE.id,VCFCALL.value"

name         value       count(*)  
-----------  ----------  ----------
rmdup_1.bam  0/1         545       
rmdup_1.bam  1/1         429       
rmdup_2.bam  0/1         625       
rmdup_2.bam  1/1         349       
rmdup_3.bam  0/1         595       
rmdup_3.bam  1/1         379       
rmdup_4.bam  0/1         548       
rmdup_4.bam  1/1         426       
rmdup_5.bam  0/1         564       
rmdup_5.bam  1/1         410       
rmdup_6.bam  0/1         724       
rmdup_6.bam  1/1         250
That's it
Pierre

07 January 2012

A CGI-version of samtools tview.

I've created a lightweight CGI-based web-application for samtools tview. This C++ program named ngsproject.cgi uses the samtools api, it allows any user to visualize all the alignments in a given NGS project. The projects and their BAMS are defined on the server side using a simple XML document. e.g:

<?xml version="1.0"?>
<projects>
 <reference id="hg19">
  <path>/home/lindenb/samtools-0.1.18/examples/ex1.fa</path>
 </reference>
 <bam id="b1">
  <sample>Sample 1</sample>
  <path>/home/lindenb/samtools-0.1.18/examples/ex1.bam</path>
 </bam>
 <bam id="b2">
  <sample>Sample 2</sample>
  <path>/home/lindenb/samtools-0.1.18/examples/ex1.bam</path>
 </bam>
 <project id="1">
  <name>Test 1</name>
  <description>Test</description>
  <bam ref="b1"/>
  <bam ref="b2"/>
  <reference ref="hg19" />
 </project>
 <project id="2">
  <name>Test 2</name>
  <description>Test</description>
  <bam ref="b2"/>
  <reference ref="hg19" />
 </project>
</projects>

Once the CGI has been installed, the user can visualize the reads of each samples.

This tool is available in the variation toolkit at http://code.google.com/p/variationtoolkit/.

That's it.

Pierre

05 January 2012

The Variation Toolkit

During the last weeks, I've worked on an experimental C++ package named The Variation Toolkit (varkit). It was originally designed to provide some command lines equivalent to knime4bio but I've added more tools over time. Some of those tools are very simple-and-stupid ( fasta2tsv) , reinvent the wheel ("numericsplit"), are part of an answer to biostar, are some old tools (e.g. bam2wig) that have been moved to this package, but some others like "samplepersnp", "groupbygene" might be useful to people.
The package is available at : http://code.google.com/p/variationtoolkit/.

Here is the current documentation (05 Jan 2012):




That's it,

Pierre