Showing posts with label parsing. Show all posts
Showing posts with label parsing. Show all posts

05 February 2013

Making use of Picard Metrics files using XML and XSLT. #ngs

Many tools in the Picard package produce some "Metrics File" (described at http://picard.sourceforge.net/picard-metric-definitions.shtml). The picard API contains a java parser "MetricsFile" parsing those metrics-file:

MetricsFile<MetricBase, Comparable<?>> metricsFile=new MetricsFile<MetricBase, Comparable<?>>();
metricsFile.read(new FileReader("metrics.txt"));
In order produce some custom reports from those files, I've created a tool that dump the content of the MetricsFile as a XML file. The source code is available at: http://code.google.com/p/jvarkit/source/browse/trunk/src/main/java/fr/inserm/umr1087/jvarkit/tools/picard/metrics2xml/PicardMetricsToXML.java.

Compilation

$ mkdir tmp
$ javac -d tmp -cp  /path/to/picard.jar:/path/to/sam.jar \
     -sourcepath  src/main/java \
     src/main/java/fr/inserm/umr1087/jvarkit/tools/picard/metrics2xml/PicardMetricsToXML.java
$ jar vcf picardmetrics2xml.jar -C tmp .

Usage

Say you have used the tool 'CollectInsertSizeMetrics.jar' from picard:
$ java -jar/path/to/CollectInsertSizeMetrics.jar \
 O=out.metrics \
 I=/path/to/samtools/examples/sorted.bam \
 AS=true \
 R=/path/to/samtools/ex1.fa \
 H=chart.pdf
The file out.metrics looks like this:
## net.sf.picard.metrics.StringHeader
# net.sf.picard.analysis.CollectInsertSizeMetrics HISTOGRAM_FILE=(...)
## net.sf.picard.metrics.StringHeader
# Started on: Tue Feb 05 12:51:30 CET 2013

## METRICS CLASS net.sf.picard.analysis.InsertSizeMetrics
MEDIAN_INSERT_SIZE MEDIAN_ABSOLUTE_DEVIATION MIN_INSERT_SIZE MAX_INSERT_SIZE MEAN_INSERT_SIZE STANDARD_DEVIATION READ_PAIRS
209 10 54 243 208.857506 13.614603 4716 FR 5 9 13 17 21 25 29 35 43 

## HISTOGRAM java.lang.Integer
insert_size All_Reads.fr_count
54 3
170 3
173 9
174 3
175 3
177 6
(...)
This file can be converted to XML using the following command:
$ java -cp /path/to/picard.jar:/path/to/sam.jar:picardmetrics2xml.jar file.metrics


<?xml version="1.0" encoding="UTF-8"?><picard-metrics xmlns="http://picard.sourc
eforge.net/" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"><metrics-file
 file="file.metrics"><headers><header class="net.sf.picard.metrics.StringHeader"
>net.sf.picard.analysis.CollectInsertSizeMetrics HISTOGRAM_FILE=jeter2 INPUT=/ho
me/lindenb/package/samtools-0.1.18/examples/sorted.bam OUTPUT=jeter REFERENCE_SE
QUENCE=/home/lindenb/package/samtools-0.1.18/examples/ex1.fa ASSUME_SORTED=true 
   DEVIATIONS=10.0 MINIMUM_PCT=0.05 METRIC_ACCUMULATION_LEVEL=[ALL_READS] STOP_A
FTER=0 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL
=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false</header><h
eader class="net.sf.picard.metrics.StringHeader">Started on: Tue Feb 05 12:51:30
 CET 2013</header></headers><metrics><thead class="net.sf.picard.analysis.Insert
SizeMetrics"><th class="double">MEDIAN_INSERT_SIZE</th><th class="double">MEDIAN
_ABSOLUTE_DEVIATION</th><th class="int">MIN_INSERT_SIZE</th><th class="int">MAX_
INSERT_SIZE</th><th class="double">MEAN_INSERT_SIZE</th><th class="double">STAND
ARD_DEVIATION</th><th class="long">READ_PAIRS</th><th class="net.sf.picard.sam.S
amPairUtil$PairOrientation">PAIR_ORIENTATION</th><th class="int">WIDTH_OF_10_PER
CENT</th><th class="int">WIDTH_OF_20_PERCENT</th><th class="int">WIDTH_OF_30_PER
CENT</th><th class="int">WIDTH_OF_40_PERCENT</th><th class="int">WIDTH_OF_50_PER
CENT</th><th class="int">WIDTH_OF_60_PERCENT</th><th class="int">WIDTH_OF_70_PER
CENT</th><th class="int">WIDTH_OF_80_PERCENT</th><th class="int">WIDTH_OF_90_PER
CENT</th><th class="int">WIDTH_OF_99_PERCENT</th><th class="java.lang.String">SA
MPLE</th><th class="java.lang.String">LIBRARY</th><th class="java.lang.String">R
EAD_GROUP</th></thead><tbody><tr><td>209.0</td><td>10.0</td><td>54</td><td>243</
td><td>208.857506</td><td>13.614603</td><td>4716</td><td>FR</td><td>5</td><td>9<
/td><td>13</td><td>17</td><td>21</td><td>25</td><td>29</td><td>35</td><td>43</td
><td>65</td><td xsi:nil="true"/><td xsi:nil="true"/><td xsi:nil="true"/></tr></t
body></metrics><histogram class="java.lang.Integer"><thead><th>insert_size</th><
th>All_Reads.fr_count</th></thead><tbody><tr><td>54</td><td>3.0</td></tr><tr><td
>170</td><td>3.0</td></tr><tr><td>173</td><td>9.0</td></tr><tr><td>174</td><td>3
.0</td></tr><tr><td>175</td><td>3.0</td></tr><tr><td>177</td><td>6.0</td></tr><t
r><td>178</td><td>6.0</td></tr><tr><td>179</td><td>9.0</td></tr><tr><td>180</td>
<td>6.0</td></tr><tr><td>181</td><td>6.0</td></tr><tr><td>182</td><td>21.0</td><
/tr><tr><td>183</td><td>9.0</td></tr><tr><td>184</td><td>15.0</td></tr><tr><td>1
85</td><td>33.0</td></tr><tr><td>186</td><td>15.0</td></tr><tr><td>187</td><td>3
(...)

Converting to JSON

Now, we can convert the XML to whatever we want using XSLT. I wrote a stylesheet picardmetrics2json.xsl converting the XML to JSON (though, I should escape the quotes in the strings ).
$ xsltproc picardmetrics2json.xsl metrics.xml


{
    "metrics.xml": {
        "headers": [
            {
                "class": "net.sf.picard.metrics.StringHeader",
                "value": "net.sf.picard.analysis.CollectInsertSizeMetrics HISTOGRAM_FILE=metrics.pdf INPUT=samtools-0.1.18/examples/sorted.bam OUTPUT=metrics.txt REFERENCE_SEQUENCE=/home/lindenb/package/samtools-0.1.18/examples/ex1.fa ASSUME_SORTED=true    DEVIATIONS=10.0 MINIMUM_PCT=0.05 METRIC_ACCUMULATION_LEVEL=[ALL_READS] STOP_AFTER=0 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false"
            },
            {
                "class": "net.sf.picard.metrics.StringHeader",
                "value": "Started on: Tue Feb 05 12:51:30 CET 2013"
            }
        ],
        "metrics": [
            {
                "MEDIAN_INSERT_SIZE": 209,
                "MEDIAN_ABSOLUTE_DEVIATION": 10,
                "MIN_INSERT_SIZE": 54,
                "MAX_INSERT_SIZE": 243,
                "MEAN_INSERT_SIZE": 208.857506,
                "STANDARD_DEVIATION": 13.614603,
                "READ_PAIRS": 4716,
                "PAIR_ORIENTATION": "FR",
                "WIDTH_OF_10_PERCENT": 5,
                "WIDTH_OF_20_PERCENT": 9,
                "WIDTH_OF_30_PERCENT": 13,
                "WIDTH_OF_40_PERCENT": 17,
                "WIDTH_OF_50_PERCENT": 21,
                "WIDTH_OF_60_PERCENT": 25,
                "WIDTH_OF_70_PERCENT": 29,
                "WIDTH_OF_80_PERCENT": 35,
                "WIDTH_OF_90_PERCENT": 43,
                "WIDTH_OF_99_PERCENT": 65,
                "SAMPLE": null,
                "LIBRARY": null,
                "READ_GROUP": null
            }
        ],
        "histogram": [
            {
                "insert_size": 54,
                "All_Reads.fr_count": 3
            },
            {
                "insert_size": 170,
                "All_Reads.fr_count": 3
            },(...)

Converting to HTML

Another stylesheet convert the XML to HTML. It also produces the javascript code to display the histograms using Google chart:
$ xsltproc picardmetrics2html.xsl metrics.xml > output.html


That's it,
Pierre

09 January 2013

A XML schema (xsd) for GeneOntology

The GeneOntology can be downloaded as a RDF/XML file from http://archive.geneontology.org/latest-termdb/go_daily-termdb.rdf-xml.gz.
Although it is a RDF file, the structure of the file remains the same. As a consequence, it is shipped with a DTD that describes the structure of the document ( http://www.geneontology.org/dtd/go.dtd ).
I've just written a XML schema (XSD) for this RDF file. This schema is available on github at:
https://github.com/lindenb/xsd-sandbox/tree/master/schemas/bio/go.

Validation with xmllint

The RDF file is successfully validated against my xsd schema:
$ curl "http://archive.geneontology.org/latest-termdb/go_daily-termdb.rdf-xml.gz" |\
 gunzip -c | grep -v "<!DOCTYPE " > go.xml

xmllint  --noout --schema go.xsd go.xml
go.xml validates
Note: I've ignored the elements defined in the DTD but absent in the RDF file.

Code Generation with XJC

XJC can be used to generate the java classes for this schema:
xjc go.xsd 
parsing a schema...
compiling a schema...
org/w3/_1999/_02/_22_rdf_syntax_ns_/ObjectFactory.java
org/w3/_1999/_02/_22_rdf_syntax_ns_/RDF.java
org/w3/_1999/_02/_22_rdf_syntax_ns_/package-info.java
org/geneontology/dtds/go/AbstractRelation.java
org/geneontology/dtds/go/Go.java
org/geneontology/dtds/go/IsA.java
org/geneontology/dtds/go/NegativelyRegulates.java
org/geneontology/dtds/go/ObjectFactory.java
org/geneontology/dtds/go/PartOf.java
org/geneontology/dtds/go/PositivelyRegulates.java
org/geneontology/dtds/go/Regulates.java
org/geneontology/dtds/go/package-info.java

Java Parsing

... and we can parse the terms of GO with java without writing a new parser and without any dependencies. For example, the following code parses the whole ontology and prints it to stdout as XML:
import java.io.InputStream;
import java.io.StringWriter;
import org.geneontology.dtds.go.*;
import org.w3._1999._02._22_rdf_syntax_ns_.*;
import javax.xml.namespace.QName;
import javax.xml.bind.JAXBContext;
import javax.xml.bind.JAXBElement;
import javax.xml.bind.Unmarshaller;
import javax.xml.bind.Marshaller;
import javax.xml.transform.stream.StreamSource;

public class TestGo
    {
    public static void main(String[] args) throws Exception
        {
 JAXBContext jaxbCtxt=JAXBContext.newInstance("org.geneontology.dtds.go:org.w3._1999._02._22_rdf_syntax_ns_");
 Marshaller marshaller = jaxbCtxt.createMarshaller();
 Unmarshaller unmarshaller=jaxbCtxt.createUnmarshaller();
        marshaller.setProperty("jaxb.formatted.output",true);
        Object go=unmarshaller.unmarshal(new java.io.File("go.xml"));
        marshaller.marshal(go, System.out);
        }
    }
compile and run:
$javac TestGo.java \
  org/w3/_1999/_02/_22_rdf_syntax_ns_/ObjectFactory.java \
  org/geneontology/dtds/go/ObjectFactory.java

$ java TestGo  | head -n 100
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<go xmlns="http://www.geneontology.org/dtds/go.dtd#" xmlns:ns2="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
    <ns2:RDF>
        <term ns2:about="http://www.geneontology.org/go#GO:0000001">
            <accession>GO:0000001</accession>
            <name>mitochondrion inheritance</name>
            <synonym>mitochondrial inheritance</synonym>
            <definition>The distribution of mitochondria, including the mitochondrial genome, into daughter cells after mitosis or meiosis, mediated by interactions between mitochondria and the cytoskeleton.</definition>
            <is_a ns2:resource="http://www.geneontology.org/go#GO:0048308"/>
            <is_a ns2:resource="http://www.geneontology.org/go#GO:0048311"/>
        </term>
        <term ns2:about="http://www.geneontology.org/go#GO:0000002">
            <accession>GO:0000002</accession>
            <name>mitochondrial genome maintenance</name>
            <definition>The maintenance of the structure and integrity of the mitochondrial genome; includes replication and segregation of the mitochondrial chromosome.</definition>
            <is_a ns2:resource="http://www.geneontology.org/go#GO:0007005"/>
            <dbxref ns2:parseType="Resource">
                <database_symbol>InterPro</database_symbol>
                <reference>IPR009446</reference>

That's it,

Pierre

PS: many thanks to @bdoughan for his help on SO.

21 November 2012

visualizing the dependencies in a Makefile

Update 2014: I wrote a C version at https://github.com/lindenb/makefile2graph.
I've just coded a tool to visualize the dependencies in a Makefile. The java source code is available on github at : https://github.com/lindenb/jsandbox/blob/master/src/sandbox/MakeGraphDependencies.java. This simple tool parses the ouput of
make -dq
( here option '-d' is 'Print lots of debugging information' and '-q' is 'Run no commands') and prints a graphiz-dot file.

Example

Below is a simple NGS workflow:
%.bam.bai : %.bam
 
file.vcf:  merged.bam.bai ref.fa
merged.bam : sorted1.bam sorted2.bam
sorted1.bam: lane1_1.fastq  lane1_2.fastq ref.fa
sorted2.bam: lane2_1.fastq  lane2_2.fastq ref.fa
Invoking the program:
make -d --dry-run | java -jar makegraphdependencies.jar
generates the following graphiz-dot file:
digraph G {
n9[label="sorted2.bam" ];
n3[label="merged.bam.bai" ];
n10[label="lane2_1.fastq" ];
n11[label="lane2_2.fastq" ];
n2[label="file.vcf" ];
n4[label="merged.bam" ];
n6[label="lane1_1.fastq" ];
n8[label="ref.fa" ];
n7[label="lane1_2.fastq" ];
n0[label="[ROOT]" ];
n5[label="sorted1.bam" ];
n1[label="Makefile" ];
n10->n9;
n11->n9;
n8->n9;
n4->n3;
n3->n2;
n8->n2;
n9->n4;
n5->n4;
n2->n0;
n1->n0;
n6->n5;
n8->n5;
n7->n5;
}
The result: (here using the google chart API for Graphviz)

That's it,
Pierre

13 July 2012

Parsing the Newick format in C using flex and bison.

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


Example:

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

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


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

   descendant_list ==> ( subtree { , subtree } )

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

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

                 label ==> unquoted_label
                       ==> quoted_label

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

         branch_length ==> signed_number
                       ==> unsigned_number

The Flex Lexer

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

The Bison Scanner

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

Makefile


Testing

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

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


That's it,

Pierre

06 September 2011

Parsing a BAM file with javascript, yes we can. (Node.js and V8)

Node.js is an event-driven I/O server-side JavaScript environment based on V8, Google's open source JavaScript engine. In the current post I will describe how I've used Node/V8 to parse a BAM file. Here I've used node v0.5.5 and my code is hosted in a new git repository bionode.

Designing a C++ native extension for Node.js wrapping the bgzf format

BAM files are stored using the bgzf format. We must first create a C++ extension wrapping the methods related to bgzf. This process is nicely described in "Writing Node.js Native Extensions". Here, a BGZF* pointer is wrapped into a class BGZFSupport that extends v8::ObjectWrapp:
class BGZFSupport: public ObjectWrap
 {
 private:
   BGZF* file;
 public:
    (...)
   int close()
    {
    int ret=0;
    if(file!=NULL) ret=::bgzf_close(file);
    file=NULL;
    return ret;
    }
   ~BGZFSupport()
    {
   if(file!=NULL) ::bgzf_close(file);
    }
 (...)
The javascript constructor for BGZFSupport opens the bgzfile and is implemented on the C++ side as:
  static Handle<Value> New(const Arguments& args)
    {
    HandleScope scope;
    if (args.Length() < 2)
      {
      RETURN_THROW("Expected two parameters for bgfz");
      }
    if(!args[0]->IsString())
     {
     RETURN_THROW("1st argument is not a string");
     }
    if(!args[1]->IsString())
     {
     RETURN_THROW("2nd argument is not a string");
     }
    
    v8::String::Utf8Value filename(args[0]);
    v8::String::Utf8Value mode(args[1]);
    BGZF* file= ::bgzf_open(ToCString(filename),ToCString(mode));
    if(file==NULL)
     {
     RETURN_THROW("Cannot open \"" << ToCString(filename) <<  "\"");
     }
    BGZFSupport* instance = new BGZFSupport(file);
    instance->Wrap(args.This());
    return args.This();
    }
... and so on for the other functions...

Implementing the javascript-based BAM-Reader

Next, we can embbed this BGZFSupport in a javascript file that will read a BAM file:
var bgzf=require("bgzf");
and we create a javascript class/function BamReader that will open the file as bgzf and will read the BAM header:
var bgzf=require("bgzf");
var Buffer = require('buffer').Buffer;


function BamReader(path)
 {
 this.fd= new bgzf.bgzf(path,"r");
 var b=new Buffer(4);
 var n = this.fd.read(b,0,4);
 if(n!=4) throw new Error("Cannot read 4 bytes");
 if(b[0]!=66)  throw new Error("Error MAGIC[0]");
 if(b[1]!=65)  throw new Error("Error MAGIC[1] got"+b[1]);
 if(b[2]!="M".charCodeAt(0))  throw new Error("Error MAGIC[2]");
 if(b[3]!="\1".charCodeAt(0))  throw new Error("Error MAGIC[3]");
 
 /* l_text */
 n = this.fd.read(b,0,4);
 if(n!=4) throw new Error("Cannot read 4 bytes");
 var l_text=b.readInt32LE(0);
 b=new Buffer(l_text);
 n = this.fd.read(b,0,l_text);
 if(n!=l_text) throw new Error("Cannot read "+l_text+" bytes (l_text)");
 this.text=b.toString('utf-8', 0, l_text);
 
 /* n_seq */
 b=new Buffer(4);
 n = this.fd.read(b,0,4);
 if(n!=4) throw new Error("Cannot read 4 bytes");
 var n_ref=b.readInt32LE(0);
 this.references=[];
 this.name2seq={};
 for(var i=0;i< n_ref;++i)
  {
  var refseq={};
  /* l_name */
  b=new Buffer(4);
  n = this.fd.read(b,0,4);
  if(n!=4) throw new Error("Cannot read 4 bytes");
  var l_name=b.readInt32LE(0);
  /* name */
  b=new Buffer(l_name);
  n = this.fd.read(b,0,l_name);
  if(n!=l_name) throw new Error("Cannot read "+l_name+" bytes (name)");
  refseq.name=b.toString('utf-8', 0,l_name-1);//\0 terminated
  /* l_ref */
  b=new Buffer(4);
  n = this.fd.read(b,0,4);
  if(n!=4) throw new Error("Cannot read 4 bytes");
  refseq.l_ref=b.readInt32LE(0);
  this.references.push(refseq);
  this.name2seq[refseq.name]=refseq;
  }
 //console.log(this.name2seq);
 }
Another function next() reads the next alignment or returns null ( see the code ).

Testing

$ export NODE_PATH=/path/to/bionode/build

the script reads a simple BAM file and prints the positions of the reads:
(...)

var r= new BamReader("/path/to/samtools-0.1.17/examples/toy.bam");
var align;
while((align=r.next())!=null)
 {
 console.log(
  r.references[align.refID].name+"\t"+
  align.read_name+"\t"+
  align.pos
  );
 }
r.close();

Result

$ node bgzf.js
ref r001 6
ref r002 8
ref r003 8
ref r004 15
ref r003 28
ref r001 36
ref2 x1 0
ref2 x2 1
ref2 x3 5
ref2 x4 9
ref2 x5 11
ref2 x6 13


Remaining questions:

At the moment, I don't know how to correctly package the C++ and javascript files for node.js, how to correctly include the files, how to group the different files under a common 'namespace', etc...

That's It,
Pierre

30 March 2011

Parsing a genomic position with javacc

Parsing a genomic position (chrom:start-end) is an easy task but I've always been too lazy to create a library for this. Today I wrote a Java-CC-based parser for analyzing the various syntaxes of a genomic position. Here is the grammar I used:

COMMA: ","
LETTER: (["a"-"z"]|["A"-"Z"]|"_") ;
DIGIT: ["0"-"9"];
INT:<DIGIT> ( (<DIGIT>|<COMMA>)* <DIGIT>)? ;
BP: "b" ("p")? ;
KB: ("k") ("B")? ;
MB: ("m") ("B")? ;
GB: ("g") ("B")? ;
IDENTIFIER: <LETTER> (<DIGIT>|<LETTER>)* ;
COLON: ":" ;
DASH: "-" ;
PLUS: "+" ;
DELIM: ("|"|";") ;


java.util.List<Segment> many(): segment() ((<DELIM>)? segment() )* )? <EOF>);
Segment one(): segment() <EOF>;
Segment segment(): chromName() <COLON> position() (<DASH> position()| <PLUS> position())? );
BigInteger position():integer() (factor())?;
BigInteger factor(): ( <BP> | <KB>| <MB> | <GB> );
BigInteger integer():<INT> ;
String chromName():( integer() | identifier());
String identifier(): <IDENTIFIER> ;

Source code



Compiling

javacc SegmentParser.jj
javac SegmentParser.java

Running

echo " chrM:1-100,000"| java SegmentParser
chrM:1-100000
echo " c1:1000"| java SegmentParser
c1:1000-1001
echo "2:1Gb+1 " | java SegmentParser
chr2:999999999-1000000002
echo "chr2:10+100" | java SegmentParser
ParseException: -90 < 0)
echo "chrX:3147483647" | java SegmentParser
ParseException: 3147483647 > 2147483647 (int-max)
echo "2:1Gb+a azd " | java SegmentParser
ParseException: Encountered "a" at line 1, column 7


That's it,

Pierre

22 February 2011

A flex scanner extracting the metadata from a PDF file.

4 years ago, I played with the adobe XMP library to extract the XMP metadata from a set of PDF files.

Today, I was suprised to simply display the XMP data contained in a PDF from Nature by using the following command line:

curl -s "http://www.nature.com/nrcardio/journal/v8/n2/pdf/nrcardio.2010.184.pdf" |\
strings |\
grep -A 100 "<x:xmpmeta"


I've generalized this process by implementing a GNU-flex scanner that prints the XML content between two <xmp:xmpmeta/> tags. The source code is available on github at: https://github.com/lindenb/ccsandbox/blob/master/src/xmpextractor.l.

Compilation

flex -f -B --read xmpextractor.l
gcc -o xmpextractor lex.yy.c

Testing with "PLOS"

Harper MA, Chen Z, Toy T, Machado IMP, Nelson SF, et al. (2011) Phenotype Sequencing: Identifying the Genes That Cause a Phenotype Directly from Pooled Sequencing of Independent Mutants. PLoS ONE 6(2): e16517. doi:10.1371/journal.pone.0016517:

curl -s "http://www.plosone.org/article/fetchObjectAttachment.action?uri=info%3Adoi%2F10.1371%2Fjournal.pone.0016517&representation=PDF" |\
./xmpextractor

<?xml version="1.0" encoding="UTF-8"?>
<XMP>
<x:xmpmeta xmlns:x="adobe:ns:meta/" x:xmptk="XMP toolkit 2.9.1-13, framework 1.6">
<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:iX="http://ns.adobe.com/iX/1.0/">
<rdf:Description xmlns:pdf="http://ns.adobe.com/pdf/1.3/" rdf:about="uuid:a0b7f786-d005-411e-b6e8-61fac5a69e23" pdf:Producer="Acrobat Distiller 7.0 (Windows)"/>
<rdf:Description xmlns:xap="http://ns.adobe.com/xap/1.0/" rdf:about="uuid:a0b7f786-d005-411e-b6e8-61fac5a69e23" xap:CreateDate="2011-02-14T08:21:58+08:00" xap:CreatorTool="3B2 Total Publishing System 7.51n/W" xap:ModifyDate="2011-02-17T14:32:08+08:00" xap:MetadataDate="2011-02-17T14:32:08+08:00"/>
<rdf:Description xmlns:xapMM="http://ns.adobe.com/xap/1.0/mm/" rdf:about="uuid:a0b7f786-d005-411e-b6e8-61fac5a69e23" xapMM:DocumentID="uuid:f98295b5-0980-42f2-884e-6cecd2d75c90" xapMM:InstanceID="uuid:d226393a-bfaf-4a2b-8c28-08215008694e"/>
<rdf:Description xmlns:dc="http://purl.org/dc/elements/1.1/" rdf:about="uuid:a0b7f786-d005-411e-b6e8-61fac5a69e23" dc:format="application/pdf">
<dc:title>
<rdf:Alt>
<rdf:li xml:lang="x-default">pone.0016517 1..16</rdf:li>
</rdf:Alt>
</dc:title>
</rdf:Description>
</rdf:RDF>
</x:xmpmeta>
</XMP>

Hum, nothing really interesting here.

Testing with "Nature Reviews Cardiology"

A test with Percutaneous coronary intervention in the elderly. Nature Reviews Cardiology 8, 79 (2010). doi:10.1038/nrcardio.2010.184
curl -s "http://www.nature.com/nrcardio/journal/v8/n2/pdf/nrcardio.2010.184.pdf" |\
./xmpextractor


<?xml version="1.0" encoding="UTF-8"?>
<XMP><x:xmpmeta xmlns:x="adobe:ns:meta/" x:xmptk="Adobe XMP Core 4.2.1-c041 52.342996, 2008/05/07-20:48:00 ">
<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
<rdf:Description rdf:about="doi:10.1038/nrcardio.2010.184"
xmlns:dc="http://purl.org/dc/elements/1.1/">
<dc:format>application/pdf</dc:format>
<dc:identifier>doi:10.1038/nrcardio.2010.184</dc:identifier>
<dc:creator>
<rdf:Seq>
<rdf:li>Tracy Y. Wang</rdf:li>
<rdf:li>Antonio Gutierrez</rdf:li>
<rdf:li>Eric D. Peterson</rdf:li>
</rdf:Seq>
</dc:creator>
<dc:description>
<rdf:Alt>
<rdf:li xml:lang="x-default">Nature Reviews Cardiology 8, 79 (2010). doi:10.1038/nrcardio.2010.184</rdf:li>
</rdf:Alt>
</dc:description>
<dc:publisher>
<rdf:Bag>
<rdf:li>Nature Publishing Group</rdf:li>
</rdf:Bag>
</dc:publisher>
<dc:rights>
<rdf:Alt>
<rdf:li xml:lang="x-default">&#xA; © 2010 Nature Publishing Group, a division of Macmillan Publishers Limited. All Rights Reserved.</rdf:li>
</rdf:Alt>
</dc:rights>
<dc:title>
<rdf:Alt>
<rdf:li xml:lang="x-default">Percutaneous coronary intervention in the elderly</rdf:li>
</rdf:Alt>
</dc:title>
</rdf:Description>
<rdf:Description rdf:about="doi:10.1038/nrcardio.2010.184"
xmlns:pdf="http://ns.adobe.com/pdf/1.3/">
<pdf:Producer>Adobe PDF Library 8.0</pdf:Producer>
</rdf:Description>
<rdf:Description rdf:about="doi:10.1038/nrcardio.2010.184"
xmlns:prism="http://prismstandard.org/namespaces/basic/2.0/">
<prism:copyright>© 2010 Nature Publishing Group</prism:copyright>
<prism:doi>10.1038/nrcardio.2010.184</prism:doi>
<prism:eIssn>1759-5010</prism:eIssn>
<prism:endingPage>90</prism:endingPage>
<prism:issn>1759-5002</prism:issn>
<prism:number>2</prism:number>
<prism:publicationName>Nature Publishing Group</prism:publicationName>
<prism:rightsAgent>permissions@nature.com</prism:rightsAgent>
<prism:startingPage>79</prism:startingPage>
<prism:volume>8</prism:volume>
<prism:publicationDate>
<rdf:Bag>
<rdf:li>2010-12-07</rdf:li>
</rdf:Bag>
</prism:publicationDate>
<prism:url>
<rdf:Bag>
<rdf:li>http://dx.doi.org/10.1038/nrcardio.2010.184</rdf:li>
</rdf:Bag>
</prism:url>
</rdf:Description>
<rdf:Description rdf:about="doi:10.1038/nrcardio.2010.184"
xmlns:xmp="http://ns.adobe.com/xap/1.0/">
<xmp:CreateDate>2011-01-10T10:09:23+05:30</xmp:CreateDate>
<xmp:CreatorTool/>
<xmp:Label>Nature Reviews Cardiology 8, 79 (2010). doi:10.1038/nrcardio.2010.184</xmp:Label>
<xmp:MetadataDate>2011-01-14T18:25:45+05:30</xmp:MetadataDate>
<xmp:ModifyDate>2011-01-14T18:25:45+05:30</xmp:ModifyDate>
<xmp:Identifier>
<rdf:Bag>
<rdf:li>doi:10.1038/nrcardio.2010.184</rdf:li>
</rdf:Bag>
</xmp:Identifier>
</rdf:Description>
<rdf:Description rdf:about="doi:10.1038/nrcardio.2010.184"
xmlns:xmpRights="http://ns.adobe.com/xap/1.0/rights/">
<xmpRights:Marked>True</xmpRights:Marked>
</rdf:Description>
<rdf:Description rdf:about="doi:10.1038/nrcardio.2010.184"
xmlns:xmpMM="http://ns.adobe.com/xap/1.0/mm/">
<xmpMM:DocumentID>uuid:51421664-c9c6-4657-9fbe-318ac969ca26</xmpMM:DocumentID>
<xmpMM:InstanceID>uuid:c65fa1fb-b6f3-4a61-8615-07b04871749e</xmpMM:InstanceID>
</rdf:Description>
</rdf:RDF>
</x:xmpmeta>
</XMP>

That's more interesting isn't it ?

That's it,

Pierre

18 February 2011

A Data Scraper for Amazonia (expression)


Today, we had a lecture about the "human induced pluripotent stem cells", presented by John De Vos. He introduced Amazonia, a free web atlas that allows an easy query of public human transcriptome data. Although there is no web service (REST/SOAP) to access this data, I was interested in getting some profiles of expression from this database as it is something I've failed to achieve with NCBI/GEO.

I wrote the following java scraper:

  • Line 84: we search for a gene name
  • 88: if there is a http redirection, the gene has been found
  • 96: the HTML page is downloaded
  • 100-112: fix the HTML to create a valid XML document
  • 133: transform the HTML page to a DOM document
  • 135-151: use XPATH to find the images and the labels
  • 189-211; put the data into a java/SWING Dialog

Compilation

javac AmazoniaRobot.java

Execution

java AmazoniaRobot EIF4G1
Et voilà:


That's it !

Pierre

17 November 2010

BLAST/XML+Annotations

I recently asked on Biostar if it would be possible to align two sequences while displaying their respective annotations.

As both answers I received (SPICE and jalview ) require a graphical interface, I quickly wrote a command-line java program doing the job. This program reads a NCBI/BLAST XML output and, if the 'query' or the 'hit' definition lines start with "gi|....", it fetches the Genbank records and the annotations for the sequence and map them onto the alignments.

The program is available on github at https://github.com/lindenb/jsandbox/blob/master/src/sandbox/BlastAnnotation.java.

Do we need an external library parsing Blast?

No, the java binding compiler, ${JAVA_HOME}/bin/xjc, can generate a java parser from BLAST DTD:
xjc -d src -p sandbox.ncbi.blast -dtd http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd

parsing a schema...
compiling a schema...
sandbox/ncbi/blast/BlastOutput.java
sandbox/ncbi/blast/BlastOutputIterations.java
sandbox/ncbi/blast/BlastOutputMbstat.java
sandbox/ncbi/blast/BlastOutputParam.java
sandbox/ncbi/blast/Hit.java
sandbox/ncbi/blast/HitHsps.java
sandbox/ncbi/blast/Hsp.java
sandbox/ncbi/blast/Iteration.java
sandbox/ncbi/blast/IterationHits.java
sandbox/ncbi/blast/IterationStat.java
sandbox/ncbi/blast/ObjectFactory.java
sandbox/ncbi/blast/Parameters.java
sandbox/ncbi/blast/Statistics.java


And do we need an external library parsing Genbank?

No, again xjc did the job:
xjc -d src -p sandbox.ncbi.gbc -dtd http://www.ncbi.nlm.nih.gov/dtd/INSD_INSDSeq.dtd

parsing a schema...
compiling a schema...
sandbox/ncbi/gbc/INSDAltSeqData.java
sandbox/ncbi/gbc/INSDAltSeqDataItems.java
sandbox/ncbi/gbc/INSDAltSeqItem.java
sandbox/ncbi/gbc/INSDAltSeqItemInterval.java
sandbox/ncbi/gbc/INSDAltSeqItemIsgap.java
sandbox/ncbi/gbc/INSDAuthor.java
sandbox/ncbi/gbc/INSDComment.java
sandbox/ncbi/gbc/INSDCommentItem.java
sandbox/ncbi/gbc/INSDCommentParagraph.java
sandbox/ncbi/gbc/INSDCommentParagraphItems.java
sandbox/ncbi/gbc/INSDCommentParagraphs.java
sandbox/ncbi/gbc/INSDFeature.java
sandbox/ncbi/gbc/INSDFeatureIntervals.java
sandbox/ncbi/gbc/INSDFeaturePartial3.java
sandbox/ncbi/gbc/INSDFeaturePartial5.java
sandbox/ncbi/gbc/INSDFeatureQuals.java
sandbox/ncbi/gbc/INSDFeatureSet.java
sandbox/ncbi/gbc/INSDFeatureSetFeatures.java
sandbox/ncbi/gbc/INSDFeatureXrefs.java
sandbox/ncbi/gbc/INSDInterval.java
sandbox/ncbi/gbc/INSDIntervalInterbp.java
sandbox/ncbi/gbc/INSDIntervalIscomp.java
sandbox/ncbi/gbc/INSDKeyword.java
sandbox/ncbi/gbc/INSDQualifier.java
sandbox/ncbi/gbc/INSDReference.java
sandbox/ncbi/gbc/INSDReferenceAuthors.java
sandbox/ncbi/gbc/INSDReferenceXref.java
sandbox/ncbi/gbc/INSDSecondaryAccn.java
sandbox/ncbi/gbc/INSDSeq.java
sandbox/ncbi/gbc/INSDSeqAltSeq.java
sandbox/ncbi/gbc/INSDSeqCommentSet.java
sandbox/ncbi/gbc/INSDSeqFeatureSet.java
sandbox/ncbi/gbc/INSDSeqFeatureTable.java
sandbox/ncbi/gbc/INSDSeqKeywords.java
sandbox/ncbi/gbc/INSDSeqOtherSeqids.java
sandbox/ncbi/gbc/INSDSeqReferences.java
sandbox/ncbi/gbc/INSDSeqSecondaryAccessions.java
sandbox/ncbi/gbc/INSDSeqStrucComments.java
sandbox/ncbi/gbc/INSDSeqid.java
sandbox/ncbi/gbc/INSDSet.java
sandbox/ncbi/gbc/INSDStrucComment.java
sandbox/ncbi/gbc/INSDStrucCommentItem.java
sandbox/ncbi/gbc/INSDStrucCommentItems.java
sandbox/ncbi/gbc/INSDXref.java
sandbox/ncbi/gbc/ObjectFactory.java

Example


As an example I've aligned the "human eif4G1" (gi|303227906) with "Mus musculus eif4G1" (gi|56699433).
The very first lines of the BLAST report are:
<?xml version="1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "NCBI_BlastOutput.dt
<BlastOutput>
<BlastOutput_program>blastn</BlastOutput_program>
<BlastOutput_version>BLASTN 2.2.24+</BlastOutput_version>
<BlastOutput_reference>Stephen F. Altschul, Thomas L. Madden, Alejandro A. Sch
<BlastOutput_db>n/a</BlastOutput_db>
<BlastOutput_query-ID>gi|303227906|ref|NM_198241.2|</BlastOutput_query-ID>
<BlastOutput_query-def>Homo sapiens eukaryotic translation initiation factor 4
<BlastOutput_query-len>5538</BlastOutput_query-len>
<BlastOutput_param>
<Parameters>
<Parameters_expect>10</Parameters_expect>
<Parameters_sc-match>2</Parameters_sc-match>
<Parameters_sc-mismatch>-3</Parameters_sc-mismatch>
<Parameters_gap-open>5</Parameters_gap-open>
<Parameters_gap-extend>2</Parameters_gap-extend>
<Parameters_filter>L;m;</Parameters_filter>
</Parameters>
</BlastOutput_param>
<BlastOutput_iterations>
<Iteration>
<Iteration_iter-num>1</Iteration_iter-num>
<Iteration_query-ID>gi|303227906|ref|NM_198241.2|</Iteration_query-ID>
<Iteration_query-def>Homo sapiens eukaryotic translation initiation factor 4 g
<Iteration_query-len>5538</Iteration_query-len>
<Iteration_hits>
<Hit>
<Hit_num>1</Hit_num>
<Hit_id>gi|56699433|ref|NM_001005331.1|</Hit_id>
<Hit_def>Mus musculus eukaryotic translation initiation factor 4, gamma 1 (Eif
<Hit_accession>NM_001005331</Hit_accession>
<Hit_len>5460</Hit_len>
<Hit_hsps>
<Hsp>
<Hsp_num>1</Hsp_num>
<Hsp_bit-score>6818.02</Hsp_bit-score>
<Hsp_score>7560</Hsp_score>
<Hsp_evalue>0</Hsp_evalue>
<Hsp_query-from>53</Hsp_query-from>
<Hsp_query-to>5538</Hsp_query-to>
<Hsp_hit-from>1</Hsp_hit-from>
<Hsp_hit-to>5418</Hsp_hit-to>
<Hsp_query-frame>1</Hsp_query-frame>
<Hsp_hit-frame>1</Hsp_hit-frame>
<Hsp_identity>4820</Hsp_identity>
<Hsp_positive>4820</Hsp_positive>
<Hsp_gaps>138</Hsp_gaps>
<Hsp_align-len>5521</Hsp_align-len>
<Hsp_qseq>GGCGCCGGCTGCGCCTGCGGAGAAGCGGTGGCCGCCGAGCGGGATCTGTGCGGGGAGCCGGAAA...
<Hsp_hseq>GGCGCTGGCTGCGCCTGCGGAGAAGCGGTGGCCGCCGAGCGGGATCTGTGCGGGGAGCCGGAAA...
<Hsp_midline>||||| |||||||||||||||||||||||||||||||||||||||||||||||||||||||...
</Hsp>
</Hit_hsps>
</Hit>
</Iteration_hits>
<Iteration_stat>
<Statistics>
<Statistics_db-num>0</Statistics_db-num>
<Statistics_db-len>0</Statistics_db-len>
<Statistics_hsp-len>0</Statistics_hsp-len>
<Statistics_eff-space>0</Statistics_eff-space>
<Statistics_kappa>-1</Statistics_kappa>
<Statistics_lambda>-1</Statistics_lambda>
<Statistics_entropy>-1</Statistics_entropy>
</Statistics>
</Iteration_stat>
</Iteration>
</BlastOutput_iterations>
</BlastOutput>

And here is the output of my program:

java -jar dist/blastannot.jar ~/jeter.blast.xml

QUERY: Homo sapiens eukaryotic translation initiation factor 4 gamma, 1 (EIF4G1), transcript variant 2, mRNA
ID:gi|303227906|ref|NM_198241.2| Len:5538
>Mus musculus eukaryotic translation initiation factor 4, gamma 1 (Eif4g1), transcript variant 2, mRNA
NM_001005331
id:gi|56699433|ref|NM_001005331.1| len:5460

e-value:0 gap:138 bitScore:6818.02

#####:############################################ exon 1..180 gene:EIF4G1
QUERY 000000053 GGCGCCGGCTGCGCCTGCGGAGAAGCGGTGGCCGCCGAGCGGGATCTGTG 000000102
||||| ||||||||||||||||||||||||||||||||||||||||||||
HIT 000000001 GGCGCTGGCTGCGCCTGCGGAGAAGCGGTGGCCGCCGAGCGGGATCTGTG 000000050
#####:############################################ exon 1..128 gene:Eif4g1



################################################## exon 1..180 gene:EIF4G1
QUERY 000000103 CGGGGAGCCGGAAATGGTTGTGGACTACGTCTGTGCGGCTGCGTGGGGCT 000000152
||||||||||||||||||||||||||||||||||||||||||||||||||
HIT 000000051 CGGGGAGCCGGAAATGGTTGTGGACTACGTCTGTGCGGCTGCGTGGGGCT 000000100
################################################## exon 1..128 gene:Eif4g1



############::::::::::###### exon 1..180 gene:EIF4G1
#:::::::::::::###::::: exon 181..237 gene:EIF4G1
QUERY 000000153 CGGCCGCGCGGACTGAAGGAGACTGAAGGCCCTCGGATGCCCAGAACCTG 000000202
|||||||||||| ||||||| |||
HIT 000000101 CGGCCGCGCGGA----------CTGAAGG-------------AGA----- 000000122
############----------#######-------------### gene 1..5460 gene:Eif4g1
############----------#######-------------### exon 1..128 gene:Eif4g1



::::::::::::::::::::::##:##:::::::# exon 181..237 gene:EIF4G1
############### exon 238..331 gene:EIF4G1
QUERY 000000203 TAGGCCGCACCGTGGACTTGTTCTTAATCGAGGGGGTGCTGGGGGGACCC 000000252
|| || ||||||||||||||||
HIT 000000123 ----------------------CTGAA-------GGTGCTGGGGGGACCC 000000143
----------------------##:##-------# exon 1..128 gene:Eif4g1
############### exon 129..222 gene:Eif4g1



#:###############################:###:############ exon 238..331 gene:EIF4G1
##############:###:############ CDS 272..5071 gene:EIF4G1
QUERY 000000253 TGATGTGGCACCAAATGAAATGAACAAAGCTCCACAGTCCACAGGCCCCC 000000302
| ||||||||||||||||||||||||||||||| ||| ||||||||||||
HIT 000000144 TAATGTGGCACCAAATGAAATGAACAAAGCTCCCCAGCCCACAGGCCCCC 000000193
#:###############################:###:############ exon 129..222 gene:Eif4g1
##############:###:############ CDS 163..4944 gene:Eif4g1


(...)


############:#:#:#####:######:#:########:##:###### exon 4890..5521 gene:EIF4G1
############:#:#:#####:######:#:########:##:###### STS 4948..5505 gene:EIF4G1
############:#:#:#####:######:#:########:##:###### STS 5174..5403 gene:EIF4G1
QUERY 000005319 TTGGTGTGTCTTGGGGTGGGGAGGGGCACCAACGCCTGCCCCTGGGGTCC 000005368
|||||||||||| | | ||||| |||||| | |||||||| || ||||||
HIT 000005201 TTGGTGTGTCTTTGCGGGGGGAAGGGCACTACCGCCTGCCTCTAGGGTCC 000005250
############:#:#:#####:######:#:########:##:###### exon 4760..5396 gene:Eif4g1



::##############:##########:###################### exon 4890..5521 gene:EIF4G1
::##############:##########:###################### STS 4948..5505 gene:EIF4G1
::##############:##########:####### STS 5174..5403 gene:EIF4G1
QUERY 000005369 TTTTTTTTATTTTCTGAAAATCACTCTCGGGACTGCCGTCCTCGCTGCTG 000005418
|||||||||||||| |||||||||| ||||||||||||||||||||||
HIT 000005251 --TTTTTTATTTTCTG-AAATCACTCTTGGGACTGCCGTCCTCGCTGCTG 000005297
--##############-##########:###################### exon 4760..5396 gene:Eif4g1



######################:#############:############# exon 4890..5521 gene:EIF4G1
######################:#############:############# STS 4948..5505 gene:EIF4G1
QUERY 000005419 GGGGCATATGCCCCAGCCCCTGTACCACCCCTGCTGTTGCCTGGGCAGGG 000005468
|||||||||||||||||||||| ||||||||||||| |||||||||||||
HIT 000005298 GGGGCATATGCCCCAGCCCCTGCACCACCCCTGCTGCTGCCTGGGCAGGG 000005347
######################:#############:############# exon 4760..5396 gene:Eif4g1



#:##-############################################: exon 4890..5521 gene:EIF4G1
#:##-################################# STS 4948..5505 gene:EIF4G1
###### polyA_signal 5496..5501 gene:EIF4G1
# polyA_site 5516 gene:EIF4G1
QUERY 000005469 GGAA-GGGGGGGCACGGTGCCTGTAATTATTAAACATGAATTCAATTAAG 000005517
| || ||||||||||||||||||||||||||||||||||||||||||||
HIT 000005348 GAAAGGGGGGGGCACGGTGCCTGTAATTATTAAACATGAATTCAATTAAA 000005397
#:##:############################################ exon 4760..5396 gene:Eif4g1



:::# exon 4890..5521 gene:EIF4G1
# polyA_site 5521 gene:EIF4G1
QUERY 000005518 CTCAAAAAAAAAAAAAAAAAA 000005538
||||||||||||||||||
HIT 000005398 AAAAAAAAAAAAAAAAAAAAA 000005418



That's it,
Pierre

19 April 2010

A stateful C function for R: parsing Fasta sequences

In the following post, I'll create a C extension for R. This extension will iterate over all the FASTA sequences in a file and will return a pair(name,sequence) for each sequence, that is to say that I won't store all the sequences in memory.

The C code

The structure FastaHandler holds the state of the fasta parser. It contains a pointer to the FILE, the current sequence and fasta header. We also need to save the previous line.
/** stateful structure for the fastafile */
typedef struct fastaHandler_t
{
/** input file */
FILE* in;
/** save the previous line, it will contains the next fasta header */
char* previous_line;
/** sequence name */
char *seq_name;
/** sequence dna */
char *seq_dna;
/** sequence_length */
int seq_length;
}FastaHandler,*FastaHandlerPtr;

The function fasta_open opens the fasta file and initialize a new FastaHandlerPtr. This pointer is then persisted into a 'R' variable using R_MakeExternalPtr.
/**
* open a fasta file, init the FastaHandler
*/
SEXP fasta_open(SEXP filename)
{
FastaHandlerPtr handle=NULL;
if(!isString(filename)) error("filename is not a string");
if(length(filename)!=1) error("expected only one filename");

handle=(FastaHandlerPtr)calloc(1,sizeof(FastaHandler));
if(handle==NULL)
{
error("Cannot alloc FastaHandler");
}
const char* c_filename=CHAR(STRING_ELT(filename,0));
errno=0;
handle->in= fopen( c_filename,"r");
if(handle->in==NULL)
{
error("Cannot open \"%s\": %s",c_filename,strerror(errno));
}
/** the handle is bound a R variable */
return R_MakeExternalPtr(handle, R_NilValue, R_NilValue);
}
We also need a function to close the structure. The pointer to the FastaHandlerPtr is retrieved using R_ExternalPtrAddr(R_variable):
/**
* close the FastaHandlerPtr
*/
SEXP fasta_close(SEXP r_handle)
{
FastaHandlerPtr handle = R_ExternalPtrAddr(r_handle);;
if(handle==NULL) error("handle==NULL");
if(handle->in!=NULL)
{
fclose(handle->in);
handle->in=NULL;
}
if(handle->previous_line!=NULL)
{
free(handle->previous_line);
handle->previous_line=NULL;
}
free(handle->seq_name);
handle->seq_name=NULL;
free(handle->seq_dna);
handle->seq_dna=NULL;
R_ClearExternalPtr(r_handle);
return ScalarInteger(0);
}
The function fasta_next return the next pair(name,sequence):
/**
* read the next fasta sequence
*/
SEXP fasta_next(SEXP r_handle)
{
int count=0;
char* ptr=NULL;
FastaHandlerPtr handle = R_ExternalPtrAddr(r_handle);;
if(handle==NULL) error("handle==NULL");
if(handle->in==NULL) error("handle->in==NULL");

while((ptr=readLine(handle))!=NULL)
{
if(ptr[0]=='>')
{
//this is the header of a fasta sequence
if(handle->seq_name!=NULL || handle->seq_dna!=NULL)
{
handle->previous_line=ptr;
SEXP return_value=NULL;
//if there a sequence in memory ?
if( handle->seq_name!=NULL &&
handle->seq_dna!=NULL &&
handle->seq_length>0)
{
//create the sequence
return_value=make_sequence(
&handle->seq_name[1],
handle->seq_dna
);
}
free(handle->seq_name);
free(handle->seq_dna);
handle->seq_name=NULL;
handle->seq_dna=NULL;
handle->seq_length=0;
//return the sequence if any
if(return_value!=NULL) return return_value;
}
//cleanup
handle->seq_name=ptr;
handle->seq_length=0;
handle->seq_dna=NULL;
}
else
{
int j=0;
int len=0;
//remove blank characters
while(ptr[j]!=0)
{
if(!isspace(ptr[j]))
{
ptr[len++]=ptr[j];
}
++j;
}
//enlarge the dna sequence
handle->seq_dna=realloc(
handle->seq_dna,
sizeof(char)*(handle->seq_length+1+len)
);
if(handle->seq_dna==NULL) error("cannot realloc seq_dna");
//append the new line
memcpy(&handle->seq_dna[handle->seq_length],ptr,sizeof(char)*len);
handle->seq_length+=len;
handle->seq_dna[handle->seq_length]=0;
free(ptr);
ptr=NULL;
}
}
//last sequence
SEXP return_value=NULL;
if( handle->seq_name!=NULL &&
handle->seq_dna!=NULL &&
handle->seq_length>0)
{
return_value=make_sequence(
&handle->seq_name[1],
handle->seq_dna
);
}
free(handle->seq_name);
free(handle->seq_dna);
handle->seq_name=NULL;
handle->seq_dna=NULL;
handle->seq_length=0;
return (return_value==NULL?R_NilValue:return_value);
}
This function calls readLine reading the next non-empty line:
/**
* return the next non empty line from a FastaHandlerPtr
* should be free()
*/
static char* readLine(FastaHandlerPtr handle)
{
int c;
char* ptr=NULL;
int length=0;
/* if there was a saved line, return it */
if(handle->previous_line!=NULL)
{
ptr=handle->previous_line;
handle->previous_line=NULL;
return ptr;
}
/** while the line is empty */
while(length==0)
{
int capacity=0;
//no more to read ?
if(feof(handle->in)) return NULL;
ptr=NULL;
//read each char
while((c=fgetc(handle->in))!=EOF)
{
//if this is an EOL, break
if(c=='\n') break;
//enlarge the buffer if it is too small
if(ptr==NULL || (length+2)>=capacity)
{
capacity+=500;
ptr=(char*)realloc(ptr,capacity*sizeof(char));
if(ptr==NULL) error("cannot realloc to %d bytes",capacity);
}
//append the char to the buffer
ptr[(length)++]=c;
}
//nothing was read (empty line) ? continue
if(length==0)
{
free(ptr);
ptr=NULL;
}
}
//add a end-of-string char
ptr[length]=0;
//return the line
return ptr;
}
It also invokes the function make_sequence(name,dna). This function creates a new R pair(name,sequence) :
/** create a pair list containing the name and the dna */
static SEXP make_sequence(const char* seq_name,const char* length)
{
SEXP values,names;
//the R value contains two objects
PROTECT(values = allocVector(VECSXP, 2));
//first item is the sequence name
SET_VECTOR_ELT(values, 0, mkString(seq_name));
//second item is the sequence dna
SET_VECTOR_ELT(values, 1, mkString(length));

//the labels
PROTECT(names = allocVector(STRSXP, 2));
//first label is the name
SET_STRING_ELT(names, 0, mkChar("name"));
//2nd label is the name
SET_STRING_ELT(names, 1, mkChar("sequence"));

//bind the labels to the values
setAttrib(values, R_NamesSymbol, names);
UNPROTECT(2);
return values;
}

The 'R' code

The 'C' dynamic library is loaded and each C function is bound to R using ".Call":
dyn.load(paste("libfasta", .Platform$dynlib.ext, sep=""))

fasta.open <- function(filename)
{
.Call("fasta_open", filename)
}

fasta.close <- function(handler)
{
.Call("fasta_close",handler)
}

fasta.next <- function(handler)
{
.Call("fasta_next",handler)
}

Compiling

Here is my makefile:
run:libfasta.so rs_chLG1.fas
${R_HOME}/bin/R --no-save < fasta.R

rs_chLG1.fas:
wget -O $@.gz ftp://ftp.ncbi.nih.gov/snp/organisms/bee_7460/rs_fasta/rs_chLG1.fas.gz
gunzip $@.gz

libfasta.so:fasta.c
gcc -fPIC -I -g -c -Wall -I ${R_HOME}/include fasta.c
gcc -shared -Wl,-soname,fasta.so.1 -o $@ fasta.o

Testing

I've downloaded some fasta sequences from dbSNP/Apis Mellifera. This file is open with fasta.open, we loop over each sequence and we print those sequences having a length lower than 300bp, at the end we close the stream with fasta.close:
f<-fasta.open("rs_chLG1.fas")

while(!is.null(snp <-fasta.next(f)))
{
if(nchar(as.character(snp$sequence))< 300)
print(snp)
}

fasta.close(f)

Result:
$name
[1] "gnl|dbSNP|rs44106732 rs=44106732|pos=48|len=298|taxid=7460|mol=\"genomic\"|class=1|alleles=\"C/T\"|build=127"

$sequence
[1] "GAAAATCAAAGACAATTTTTGGAATGGAACTAAATAACTTTATTCTTYCTTTCTTTCTCGTCGAGAATATCGTTATCGTTGCATGGGTTATGGAATAGCGTTCGTTAAAAATGTTATATTTCGAGGAAATATCGAAGATAGGCTTTGCGAAAGTCTGTTTCTCTAGAATTAAGATTTATATGTGTTGCAAGGGGAAGTTCAAAGAGAAATCGTGGCCAGTTCGAACATTATTATGTCTATCAATGATCGAGAAGTGTCATTGATGCAAGAGAAAAGTTTTCTCTTGCATTTAAGTATC"

$name
[1] "gnl|dbSNP|rs44108824 rs=44108824|pos=251|len=268|taxid=7460|mol=\"genomic\"|class=1|alleles=\"C/T\"|build=127"

$sequence
[1] "GTTAAATGGGAATTTTGGGGATTATTGGAGGAGGATTTACGTTTCGAGGATTGTTGATGATCTTAGGATTGTTTTCAGTTTGGAATTTTTTCTTCTTCTTCAACGTTGTACAACATTCTCCTCGAATTTTGTGTAACGAGGAGGATTAAACCTTTGGAAAATCACGTAAATTAGAGGACGATATATTGGGTTGGCAACTAAGTAATTGCGGATTTTTTTTAGAAAATCAAAGACAATTTTTGGAATGGAAYTAAATAACTTTATTCTT"

$name
[1] "gnl|dbSNP|rs44150807 rs=44150807|pos=251|len=274|taxid=7460|mol=\"genomic\"|class=1|alleles=\"C/T\"|build=127"

$sequence
[1] "TTCGATCTTCGTTTCACGCTCACGTTTCACGTTTCTCGTCCGGCGTAACGAACGGTATTCCGCTGATCACGAATAATTTCTTTCTGGAGAGTCCATTAGGGGACCGTCCCCTCCCCCCTCTCGCGCACACAGACACCCATCGCTTTCGACGCCTCGTCCATTCGAGGGAGAAACGAACGACTATTAGAAAAAAATCTTCTTTATATCTCTATAAATTCAATTTGCAGAGAAGCAAAGAGCTTTAAAATATYAACCATTATAACCGAACTTGTTG"

(....)


Full Source code


Makefile


run:libfasta.so rs_chLG1.fas
${R_HOME}/bin/R --no-save < fasta.R

rs_chLG1.fas:
wget -O $@.gz ftp://ftp.ncbi.nih.gov/snp/organisms/bee_7460/rs_fasta/rs_chLG1.fas.gz
gunzip $@.gz

libfasta.so:fasta.c
gcc -fPIC -I -g -c -Wall -I ${R_HOME}/include fasta.c
gcc -shared -Wl,-soname,fasta.so.1 -o $@ fasta.o

fasta.c


#include <ctype.h>
#include <errno.h>
#include <R.h>
#include <Rinternals.h>

/** stateful structure for the fastafile */
typedef struct fastaHandler_t
{
/** input file */
FILE* in;
/** save previous line, it will contains the next fasta header */
char* previous_line;
/** sequence name */
char *seq_name;
/** sequence dna */
char *seq_dna;
/** sequence_length */
int seq_length;

}FastaHandler,*FastaHandlerPtr;


/**
* return the next non empty line from a FastaHandlerPtr
* should be free()
*/
static char* readLine(FastaHandlerPtr handle)
{
int c;
char* ptr=NULL;
int length=0;
/* if there was a saved line, return it */
if(handle->previous_line!=NULL)
{
ptr=handle->previous_line;
handle->previous_line=NULL;
return ptr;
}
/** while the line is empty */
while(length==0)
{
int capacity=0;
//no more to read ?
if(feof(handle->in)) return NULL;
ptr=NULL;
//read each char
while((c=fgetc(handle->in))!=EOF)
{
//if this is an EOL, break
if(c=='\n') break;
//enlarge the buffer if it is too small
if(ptr==NULL || (length+2)>=capacity)
{
capacity+=500;
ptr=(char*)realloc(ptr,capacity*sizeof(char));
if(ptr==NULL) error("cannot realloc to %d bytes",capacity);
}
//append the char to the buffer
ptr[(length)++]=c;
}
//nothing was read (empty line) ? continue
if(length==0)
{
free(ptr);
ptr=NULL;
}
}
//add a end-of-string char
ptr[length]=0;
//return the line
return ptr;
}

/**
* open a fasta file, init the FastaHandler
*/
SEXP fasta_open(SEXP filename)
{
FastaHandlerPtr handle=NULL;
if(!isString(filename)) error("filename is not a string");
if(length(filename)!=1) error("expected only one filename");

handle=(FastaHandlerPtr)calloc(1,sizeof(FastaHandler));
if(handle==NULL)
{
error("Cannot alloc FastaHandler");
}
const char* c_filename=CHAR(STRING_ELT(filename,0));
errno=0;
handle->in= fopen( c_filename,"r");
if(handle->in==NULL)
{
error("Cannot open \"%s\": %s",c_filename,strerror(errno));
}
/** the handle is bound a R variable */
return R_MakeExternalPtr(handle, R_NilValue, R_NilValue);
}

/** create a pair list containing the name and the dna */
static SEXP make_sequence(const char* seq_name,const char* length)
{
SEXP values,names;
//the R value contains two objects
PROTECT(values = allocVector(VECSXP, 2));
//first item is the sequence name
SET_VECTOR_ELT(values, 0, mkString(seq_name));
//second item is the sequence dna
SET_VECTOR_ELT(values, 1, mkString(length));

//the labels
PROTECT(names = allocVector(STRSXP, 2));
//first label is the name
SET_STRING_ELT(names, 0, mkChar("name"));
//2nd label is the name
SET_STRING_ELT(names, 1, mkChar("sequence"));

//bind the labels to the values
setAttrib(values, R_NamesSymbol, names);
UNPROTECT(2);
return values;
}

/**
* read the next fasta sequence
*/
SEXP fasta_next(SEXP r_handle)
{
int count=0;
char* ptr=NULL;
FastaHandlerPtr handle = R_ExternalPtrAddr(r_handle);;
if(handle==NULL) error("handle==NULL");
if(handle->in==NULL) error("handle->in==NULL");

while((ptr=readLine(handle))!=NULL)
{
if(ptr[0]=='>')
{
//this is the header of a fasta sequence
if(handle->seq_name!=NULL || handle->seq_dna!=NULL)
{
handle->previous_line=ptr;
SEXP return_value=NULL;
//if there a sequence in memory ?
if( handle->seq_name!=NULL &&
handle->seq_dna!=NULL &&
handle->seq_length>0)
{
//create the sequence
return_value=make_sequence(
&handle->seq_name[1],
handle->seq_dna
);
}
free(handle->seq_name);
free(handle->seq_dna);
handle->seq_name=NULL;
handle->seq_dna=NULL;
handle->seq_length=0;
//return the sequence if any
if(return_value!=NULL) return return_value;
}
//cleanup
handle->seq_name=ptr;
handle->seq_length=0;
handle->seq_dna=NULL;
}
else
{
int j=0;
int len=0;
//remove blank characters
while(ptr[j]!=0)
{
if(!isspace(ptr[j]))
{
ptr[len++]=ptr[j];
}
++j;
}
//enlarge the dna sequence
handle->seq_dna=realloc(
handle->seq_dna,
sizeof(char)*(handle->seq_length+1+len)
);
if(handle->seq_dna==NULL) error("cannot realloc seq_dna");
//append the new line
memcpy(&handle->seq_dna[handle->seq_length],ptr,sizeof(char)*len);
handle->seq_length+=len;
handle->seq_dna[handle->seq_length]=0;
free(ptr);
ptr=NULL;
}
}
//last sequence
SEXP return_value=NULL;
if( handle->seq_name!=NULL &&
handle->seq_dna!=NULL &&
handle->seq_length>0)
{
return_value=make_sequence(
&handle->seq_name[1],
handle->seq_dna
);
}
free(handle->seq_name);
free(handle->seq_dna);
handle->seq_name=NULL;
handle->seq_dna=NULL;
handle->seq_length=0;
return (return_value==NULL?R_NilValue:return_value);
}

/**
* close the FastaHandlerPtr
*/
SEXP fasta_close(SEXP r_handle)
{
FastaHandlerPtr handle = R_ExternalPtrAddr(r_handle);;
if(handle==NULL) error("handle==NULL");
if(handle->in!=NULL)
{
fclose(handle->in);
handle->in=NULL;
}
if(handle->previous_line!=NULL)
{
free(handle->previous_line);
handle->previous_line=NULL;
}
free(handle->seq_name);
handle->seq_name=NULL;
free(handle->seq_dna);
handle->seq_dna=NULL;
R_ClearExternalPtr(r_handle);
return ScalarInteger(0);
}

fasta.R


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

fasta.open <- function(filename)
{
.Call("fasta_open", filename)
}
fasta.close <- function(handler)
{
.Call("fasta_close",handler)
}

fasta.next <- function(handler)
{
.Call("fasta_next",handler)
}

f<-fasta.open("rs_chLG1.fas")

while(!is.null(snp <-fasta.next(f)))
{
if(nchar(as.character(snp$sequence))< 300)
print(snp)
}

fasta.close(f)

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