18 September 2009

Translating DNA with XALAN: A custom extension for XSLT

In this post, I'll show how I've create a custom extension for XALAN , a java-based XSLT engine. My favorite XSLT processor has always been xsltproc but I was missing the capacity to create a custom function to process the XML document: Here I show how a java class can be plugged to XALAN to translate a DNA sequence to a peptide.

The Java class

The following class test.Translate translate a DNA to an amino acid sequence, the argument of the constructor is the transl_table given by the NCBI (see http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi). This is not Rocket Science.
package test;
public class Translate
{
static String geneticCode;

public Translate(String geneticCode)
{
this.geneticCode=geneticCode;
}

private static int base2index(char c)
{
switch(Character.toLowerCase(c))
{
case 't': return 0;
case 'c': return 1;
case 'a': return 2;
case 'g': return 3;
default: return -1;
}
}

public String translate(String sequence) {
StringBuilder b= new StringBuilder(1+sequence.length()/3);
for(int i=0;i+2< sequence.length();i+=3)
{
int base1= base2index(sequence.charAt(i));
int base2= base2index(sequence.charAt(i+1));
int base3= base2index(sequence.charAt(i+2));
if(base1==-1 || base2==-1 || base3==-1)
{
b.append('?');
}
else
{
b.append(geneticCode.charAt(base1*16+base2*4+base3));
}
}
return b.toString();
}
}

Compiling and packaging the source


javac test/Translate.java
jar cvf translate.jar test

The XML source


The XML source is a set of INSDSeq sequences downloaded from the NCBI/Genbank
<?xml version="1.0"?>
<!DOCTYPE INSDSet PUBLIC "-//NCBI//INSD INSDSeq/EN" "http://ww
w.ncbi.nlm.nih.gov/dtd/INSD_INSDSeq.dtd">
<INSDSet>
<INSDSeq>
<INSDSeq_locus>NM_004953</INSDSeq_locus>
<INSDSeq_length>4888</INSDSeq_length>
<INSDSeq_strandedness>single</INSDSeq_strandedness>
<INSDSeq_moltype>mRNA</INSDSeq_moltype>
<INSDSeq_topology>linear</INSDSeq_topology>
<INSDSeq_division>PRI</INSDSeq_division>
<INSDSeq_update-date>03-SEP-2009</INSDSeq_update-date>
<INSDSeq_create-date>14-MAY-1999</INSDSeq_create-date>
<INSDSeq_definition>Homo sapiens eukaryotic translation initiation factor 4 gamma, 1 (EIF4G1), transcript variant 5, mRNA</INSDSeq_definition>
<INSDSeq_primary-accession>NM_004953</INSDSeq_primary-accession>
<INSDSeq_accession-version>NM_004953.3</INSDSeq_accession-version>
<INSDSeq_other-seqids>
<INSDSeqid>ref|NM_004953.3|</INSDSeqid>
<INSDSeqid>gi|148277098</INSDSeqid>
</INSDSeq_other-seqids>
<INSDSeq_source>Homo sapiens (human)</INSDSeq_source>
<INSDSeq_organism>Homo sapiens</INSDSeq_organism>
<INSDSeq_taxonomy>Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Eu
teleostomi; Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini; Catarrh
ini; Hominidae; Homo</INSDSeq_taxonomy>
(...)
<INSDFeature_key>CDS</INSDFeature_key>
<INSDFeature_location>207..4418</INSDFeature_location>
<INSDFeature_intervals>
<INSDInterval>
<INSDInterval_from>207</INSDInterval_from>
<INSDInterval_to>4418</INSDInterval_to>
<INSDInterval_accession>NM_004953.3</INSDInterval_accession>
</INSDInterval>
</INSDFeature_intervals>
<INSDFeature_quals>
<INSDQualifier>
<INSDQualifier_name>gene</INSDQualifier_name>
<INSDQualifier_value>EIF4G1</INSDQualifier_value>
</INSDQualifier>
<INSDQualifier>
<INSDQualifier_name>gene_synonym</INSDQualifier_name>
<INSDQualifier_value>DKFZp686A1451; EIF4F; EIF4G; p220</INSDQualifier_value>
</INSDQualifier>
(...)
cttttaatgatgagggtaactatttcagttgtgagccttctagggccccaggctgggaggctcagaggactgaatctgggacctgtgttccccccggcaggcagggacaagatggcatggcaagcatgggggcggggtgggtggggagggatgctgcatttctcagctgggcagtaatcaatttaatggtcctttaaaatgtctgtgtattaaaaatttaagaataccacactttaatattaaatattcataaggtctagtatcttgataataatgtagatgttttaataacaatttttgtccttcttaaaataaaatgaaagaaacttgcttcccttagcctttgttctagaaaataaacttgtgcactttga</INSDSeq_sequence>
</INSDSeq>
</INSDSet>

The XSLT stylesseet

In the following XSLT stylesheet:
in the header, the prefix bio is associated with our Translate class
a parameter named GENETICCODE defines the default transl_table
A new Translate object named code is created with bio:new($GENETICCODE)
This code is then called on the sub-string of DNA containing the CDS: bio:translate($code,substring($dna,$start,1+($end - $start)))


<xsl:stylesheet version="1.0"
xmlns:xsl='http://www.w3.org/1999/XSL/Transform'
xmlns:bio="xalan://test.Translate"
extension-element-prefixes="bio">

<xsl:output method="xml" indent="yes"/>

<xsl:param name="GENETICCODE">FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG</xsl:param>

<xsl:template match="/INSDSet">
<html><body>
<xsl:apply-templates select="INSDSeq"/>
</body></html>
</xsl:template>


<xsl:template match="INSDSeq">
<xsl:variable name="dna" select="INSDSeq_sequence"/>
<div>
<h2><xsl:value-of select="INSDSeq_locus"/><xsl:text> : </xsl:text><xsl:value-of select="INSDSeq_definition"/></h2>
<xsl:for-each select="INSDSeq_feature-table/INSDFeature[INSDFeature_key='CDS']/INSDFeature_intervals">
<xsl:variable name="code" select="bio:new($GENETICCODE)"/>
<xsl:variable name="start" select="number(INSDInterval/INSDInterval_from)"/>
<xsl:variable name="end" select="number(INSDInterval/INSDInterval_to)"/>
<div style="font-family:monospace;word-wrap:break-word;width:400px;background-color:rgb(230,230,230);"><xsl:value-of select="bio:translate($code,substring($dna,$start,1+($end - $start)))"/></div>
</xsl:for-each>
</div>
</xsl:template>


</xsl:stylesheet>

Running this XSLT stylesheet with XALAN


java -cp ${XALAN_PATH}/org.apache.xalan_2.7.1.v200905122109.jar:${XALAN_PATH}/org.apache.xml.serializer_2.7.1.v200902170519.jar:translate.jar org.apache.xalan.xslt.Process -IN sequences.gbc -XSL seq2html.xsl

Result


NM_004953 : Homo sapiens eukaryotic translation initiation factor 4 gamma, 1
(EIF4G1), transcript variant 5, mRNA


MSGARTASTPTPPQTGGGLEPQANGETPQVAVIVRPDDRSQGAIIADRPGLPGPEHSPSESQPSSPSPTPSPSPVLEPGSEPNLAVLSIPGDTMTTIQMSVEESTPISRETGEPYRLSPEPTPLAEPILEVEVTLSKPVPESEFSSSPLQAPTPLASHTVEIHEPNGMVPSEDLEPEVESSPELAPPPACPSESPVPIAPTAQPEELLNGAPSPPAVDLSPVSEPEEQAKEVTASMAPPTIPSATPATAPSATSPAQEEEMEEEEEEEEGEAGEAGEAESEKGGEELLPPESTPIPANLSQNLEAAAATQVAVSVPKRRRKIKELNKKEAVGDLLDAFKEANPAVPEVENQPPAGSNPGPESEGSGVPPRPEEADETWDSKEDKIHNAENIQPGEQKYEYKSDQWKPLNLEEKKRYDREFLLGFQFIFASMQKPEGLPHISDVVLDKANKTPLRPLDPTRLQGINCGPDFTPSFANLGRTTLSTRGPPRGGPGGELPRGPAGLGPRRSQQGPRKEPRKIIATVLMTEDIKLNKAEKAWKPSSKRTAADKDRGEEDADGSKTQDLFRRVRSILNKLTPQMFQQLMKQVTQLAIDTEERLKGVIDLIFEKAISEPNFSVAYANMCRCLMALKVPTTEKPTVTVNFRKLLLNRCQKEFEKDKDDDEVFEKKQKEMDEAATAEERGRLKEELEEARDIARRRSLGNIKFIGELFKLKMLTEAIMHDCVVKLLKNHDEESLECLCRLLTTIGKDLDFEKAKPRMDQYFNQMEKIIKEKKTSSRIRFMLQDVLDLRGSNWVPRRGDQGPKTIDQIHKEAEMEEHREHIKVQQLMAKGSDKRRGGPPGPPISRGLPLVDDGGWNTVPISKGSRPIDTSRLTKITKPGSIDSNNQLFAPGGRLSWGKGSSGGSGAKPSDAASEAARPATSTLNRFSALQQAVPTESTDNRRVVQRSSLSRERGEKAGDRGDRLERSERGGDRGDRLDRARTPATKRSFSKEVEERSRERPSQPEGLRKAASLTEDRDRGRDAVKREAALPPVSPLKAALSEEELEKKSKAIIEEYLHLNDMKEAVQCVQELASPSLLFIFVRHGVESTLERSAIAREHMGQLLHQLLCAGHLSTAQYYQGLYEILELAEDMEIDIPHVWLYLAELVTPILQEGGVPMGELFREITKPLRPLGKAASLLLEILGLLCKSMGPKKVGTLWREAGLSWKEFLPEGQDIGAFVAEQKVEYTLGEESEAPGQRALPSEELNRQLEKLLKEGSSNQRVFDWIEANLSEQQIVSNTLVRALMTAVCYSAIIFETPLRVDVAVLKARAKLLQKYLCDEQKELQALYALQALVVTLEQPPNLLRMFFDALYDEDVVKEDAFYSWESSKDPAEQQGKGVALKSVTAFFKWLREAEEESDHN*



NM_182917 : Homo sapiens eukaryotic translation initiation factor 4 gamma, 1
(EIF4G1), transcript variant 1, mRNA


MNKAPQSTGPPPAPSPGLPQPAFPPGQTAPVVFSTPQATQMNTPSQPRQHFYPSRAQPPSSAASRVQSAAPARPGPAAHVYPAGSQVMMIPSQISYPASQGAYYIPGQGRSTYVVPTQQYPVQPGAPGFYPGASPTEFGTYAGAYYPAQGVQQFPTGVAPAPVLMNQPPQIAPKRERKTIRIRDPNQGGKDITEEIMSGARTASTPTPPQTGGGLEPQANGETPQVAVIVRPDDRSQGAIIADRPGLPGPEHSPSESQPSSPSPTPSPSPVLEPGSEPNLAVLSIPGDTMTTIQMSVEESTPISRETGEPYRLSPEPTPLAEPILEVEVTLSKPVPESEFSSSPLQAPTPLASHTVEIHEPNGMVPSEDLEPEVESSPELAPPPACPSESPVPIAPTAQPEELLNGAPSPPAVDLSPVSEPEEQAKEVTASMAPPTIPSATPATAPSATSPAQEEEMEEEEEEEEGEAGEAGEAESEKGGEELLPPESTPIPANLSQNLEAAAATQVAVSVPKRRRKIKELNKKEAVGDLLDAFKEANPAVPEVENQPPAGSNPGPESEGSGVPPRPEEADETWDSKEDKIHNAENIQPGEQKYEYKSDQWKPLNLEEKKRYDREFLLGFQFIFASMQKPEGLPHISDVVLDKANKTPLRPLDPTRLQGINCGPDFTPSFANLGRTTLSTRGPPRGGPGGELPRGPAGLGPRRSQQGPRKEPRKIIATVLMTEDIKLNKAEKAWKPSSKRTAADKDRGEEDADGSKTQDLFRRVRSILNKLTPQMFQQLMKQVTQLAIDTEERLKGVIDLIFEKAISEPNFSVAYANMCRCLMALKVPTTEKPTVTVNFRKLLLNRCQKEFEKDKDDDEVFEKKQKEMDEAATAEERGRLKEELEEARDIARRRSLGNIKFIGELFKLKMLTEAIMHDCVVKLLKNHDEESLECLCRLLTTIGKDLDFEKAKPRMDQYFNQMEKIIKEKKTSSRIRFMLQDVLDLRGSNWVPRRGDQGPKTIDQIHKEAEMEEHREHIKVQQLMAKGSDKRRGGPPGPPISRGLPLVDDGGWNTVPISKGSRPIDTSRLTKITKPGSIDSNNQLFAPGGRLSWGKGSSGGSGAKPSDAASEAARPATSTLNRFSALQQAVPTESTDNRRVVQRSSLSRERGEKAGDRGDRLERSERGGDRGDRLDRARTPATKRSFSKEVEERSRERPSQPEGLRKAASLTEDRDRGRDAVKREAALPPVSPLKAALSEEELEKKSKAIIEEYLHLNDMKEAVQCVQELASPSLLFIFVRHGVESTLERSAIAREHMGQLLHQLLCAGHLSTAQYYQGLYEILELAEDMEIDIPHVWLYLAELVTPILQEGGVPMGELFREITKPLRPLGKAASLLLEILGLLCKSMGPKKVGTLWREAGLSWKEFLPEGQDIGAFVAEQKVEYTLGEESEAPGQRALPSEELNRQLEKLLKEGSSNQRVFDWIEANLSEQQIVSNTLVRALMTAVCYSAIIFETPLRVDVAVLKARAKLLQKYLCDEQKELQALYALQALVVTLEQPPNLLRMFFDALYDEDVVKEDAFYSWESSKDPAEQQGKGVALKSVTAFFKWLREAEEESDHN*



Running this XSLT stylesheet with XALAN and another genetic code


We set the param GENETICCODE to the Scenedesmus obliquus mitochondrial Code:
java -cp ${XALAN_PATH}/org.apache.xalan_2.7.1.v200905122109.jar:${XALAN_PATH}/org.apache.xml.serializer_2.7.1.v200902170519.jar:translate.jar org.apache.xalan.xslt.Process -IN sequences.gbc -XSL seq2html.xsl -PARAM GENETICCODE "'FFLLSSSSYY*QCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG'"

Result



NM_004953 : Homo sapiens eukaryotic translation initiation factor 4 gamma, 1
(EIF4G1), transcript variant 5, mRNA


ILGARMASTPTLPQTGGELELHVTGETPQRVVRVRPADRSQGAIRVDRPGLLGPEPSLSDSQLSSLLPTPSPSPVLDPGLELTLAVLLRLGDMITMIHILVDDSTPISQDMGEPSRLLPDPMLLADPILDVDVTWSNPRPDLE'LSKLLQVPTLLALHTVDRPELTGIVPLDALDPEVESSPEWVLPPVCPSDSLVPRVPMAHLEDLLNGAPSPPVVDFSPVKEPEEQAKEVTASIAPPTIPLVMPVTVLSVMSPVQEEDIDDDDDEDDGDAGDAGDVEKENGGEDLLPPEKTLRPANLLQTLEAAAAMHVAVLVPKRSRNRKELTKKEVRGDWLAAFKEANPAVPEVDTQLLAGSTPGPELEGKEVPPQLEDAAET*DSKDDNRHTVENIQPGDQKSDSKSAQ*KLLNLEENNQYDQEFLWE'QFI'AKIQKPEGLPPIKDVVLDKATNTPLRPLAPMSLHGITQGPDFMPS'ANWGRTTWSTQGPPREGPEGELPQGPVGLGPRRLQQGPRNDPRKIRATVFITDAINLNNAENA*NPSSKRTAVAKARGDDAVAGSNTQDLFRRVRSILTNLTPQIFQQLIKHVTQLAIDTEDRLNGVRDLR'EKARSEPNFLVASANICRCLIALNVPMTDKPTVMVNFRKLLLTRQQKE'ENDNAAAER'EKKHNEIADVVTAEDRGRLKDELDEVRDIARRRLFGTIK'RGELFNLKIFTEAIIPDQVVNLWKNPADESWECWQQLLTTRGNDLD'DNAKPRIAQSFNQIDNIRNDKKTSSRIR'ILQDVLALRGST*VPRRGAQEPKTRDQIPKEVEIDDPREHINVQQLIAKGKDKQRGELPGLPISQGWPWVAAEG*NTRPISNESRPRDTSRLTKITKLGSIALNNQL'ALGGRLS*GKGSSGGSGAKPSDAASDVVRPVMKMLTRFSAWHHAVPTDSTATSQVVQRKSLSRDRGENVGDRGDRLERKDRGGDQGDRWAQARTLVTKRSFSKDVEERKSDRPSQLEGLRKAVSLTEARDQGRAAVKRDVALPPVSPLKAVLLEEEFEKNSKVIREDSLPLTDINEAVQCVQELASPSLLFI'VRPEVELTLERKARVQEPIGQLLHQLLQVGPLLMVQYYHGLSDILDLVEDIDRDIPHV*LYLADLVTPRLQDEGVPIGELFRERTKLLSPLGNVVSLLLEILGLLCNSIELNKVGTL*RDAGWS*KD'LLDGQDREAFVVDQKVESTLGEESDALGQRALPSEELNRQLEKLLKEGSKNQRVFD*IEANLKEQQIVSNTFRRALITVVCSLARR'EMPLRVDRAVLNARAKLLQNYLQDEQKELQALYALQAWVVTFDQLPNLLRIF'DALSDEDVVKEAAFYK*EKSKDPVEQQGKEVAWNLVTAFFK*LQDAEEELDHNC



NM_182917 : Homo sapiens eukaryotic translation initiation factor 4 gamma, 1
(EIF4G1), transcript variant 1, mRNA


INNVPQSTGPPPAPSPGLPQPA'PPGQTAPVVFKTPHATHINTLLQPRQHFYLSRAQPPSKAASRVQKAALARLGPVAPVYLVGSHVIIILSQISYPASQGAYYILGQGQSTYRVPTQQYLVQPGAPGFSPEASLTD'GTYVGAYSPAHGVQQ'PMGVAPAPRLINQPPQRVPKREQKTIRRRAPNHGGKAITEEIILGARMASTPTLPQTGGELELHVTGETPQRVVRVRPADRSQGAIRVDRPGLLGPEPSLSDSQLSSLLPTPSPSPVLDPGLELTLAVLLRLGDMITMIHILVDDSTPISQDMGEPSRLLPDPMLLADPILDVDVTWSNPRPDLE'LSKLLQVPTLLALHTVDRPELTGIVPLDALDPEVESSPEWVLPPVCPSDSLVPRVPMAHLEDLLNGAPSPPVVDFSPVKEPEEQAKEVTASIAPPTIPLVMPVTVLSVMSPVQEEDIDDDDDEDDGDAGDAGDVEKENGGEDLLPPEKTLRPANLLQTLEAAAAMHVAVLVPKRSRNRKELTKKEVRGDWLAAFKEANPAVPEVDTQLLAGSTPGPELEGKEVPPQLEDAAET*DSKDDNRHTVENIQPGDQKSDSKSAQ*KLLNLEENNQYDQEFLWE'QFI'AKIQKPEGLPPIKDVVLDKATNTPLRPLAPMSLHGITQGPDFMPS'ANWGRTTWSTQGPPREGPEGELPQGPVGLGPRRLQQGPRNDPRKIRATVFITDAINLNNAENA*NPSSKRTAVAKARGDDAVAGSNTQDLFRRVRSILTNLTPQIFQQLIKHVTQLAIDTEDRLNGVRDLR'EKARSEPNFLVASANICRCLIALNVPMTDKPTVMVNFRKLLLTRQQKE'ENDNAAAER'EKKHNEIADVVTAEDRGRLKDELDEVRDIARRRLFGTIK'RGELFNLKIFTEAIIPDQVVNLWKNPADESWECWQQLLTTRGNDLD'DNAKPRIAQSFNQIDNIRNDKKTSSRIR'ILQDVLALRGST*VPRRGAQEPKTRDQIPKEVEIDDPREHINVQQLIAKGKDKQRGELPGLPISQGWPWVAAEG*NTRPISNESRPRDTSRLTKITKLGSIALNNQL'ALGGRLS*GKGSSGGSGAKPSDAASDVVRPVMKMLTRFSAWHHAVPTDSTATSQVVQRKSLSRDRGENVGDRGDRLERKDRGGDQGDRWAQARTLVTKRSFSKDVEERKSDRPSQLEGLRKAVSLTEARDQGRAAVKRDVALPPVSPLKAVLLEEEFEKNSKVIREDSLPLTDINEAVQCVQELASPSLLFI'VRPEVELTLERKARVQEPIGQLLHQLLQVGPLLMVQYYHGLSDILDLVEDIDRDIPHV*LYLADLVTPRLQDEGVPIGELFRERTKLLSPLGNVVSLLLEILGLLCNSIELNKVGTL*RDAGWS*KD'LLDGQDREAFVVDQKVESTLGEESDALGQRALPSEELNRQLEKLLKEGSKNQRVFD*IEANLKEQQIVSNTFRRALITVVCSLARR'EMPLRVDRAVLNARAKLLQNYLQDEQKELQALYALQAWVVTFDQLPNLLRIF'DALSDEDVVKEAAFYK*EKSKDPVEQQGKEVAWNLVTAFFK*LQDAEEELDHNC



That's it
Pierre

04 September 2009

My Notes about Bowtie


Image via wikipedia

Bowtie is an ultrafast, memory-efficient alignment program for aligning short DNA sequence reads to large genomes. For the human genome, Burrows-Wheeler indexing allows Bowtie to align more than 25 million reads per CPU hour with a memory footprint of approximately 1.3 gigabytes. Bowtie is open source http://bowtie.cbcb.umd.edu (Langmead B, & al. "Ultrafast and memory-efficient alignment of short DNA sequences to the human genome". Genome Biol. 2009;10(3):R25.)

Building an indexed Database


The fasta sequences retrieved from the query
"ROTAVIRUS[organism] VP1 Human"
were saved to the file ./genomes/rotavirus.fna.
The file was then indexed with ./bowtie-build
./bowtie-build genomes/rotavirus.fna indexes/rotavirus

Settings:
Output files: "rotavirus.*.ebwt"
Line rate: 6 (line is 64 bytes)
Lines per side: 1 (side is 64 bytes)
Offset rate: 5 (one in 32)
FTable chars: 10
Strings: unpacked
Max bucket size: default
Max bucket size, sqrt multiplier: default
Max bucket size, len divisor: 4
Difference-cover sample period: 1024
Reference base cutoff: none
Endianness: little
Actual local endianness: little
Sanity checking: disabled
Assertions: disabled
Random seed: 0
Sizeofs: void*:4, int:4, long:4, size_t:4
Input files DNA, FASTA:
genomes/rotavirus.fna
Reading reference sizes
Time reading reference sizes: 00:00:00
Calculating joined length
= 316434 (0 characters of padding)
Writing header
Reserving space for joined string
Joining reference sequences
Time to join reference sequences: 00:00:00
bmax according to bmaxDivN setting: 79108
Using parameters --bmax 59331 --dcv 1024
Doing ahead-of-time memory usage test
Passed! Constructing with these parameters: --bmax 59331 --dcv 1024
Constructing suffix-array element generator
Building DifferenceCoverSample
Building sPrime
Building sPrimeOrder
V-Sorting samples
V-Sorting samples time: 00:00:00
Allocating rank array
Ranking v-sort output
Ranking v-sort output time: 00:00:00
Invoking Larsson-Sadakane on ranks
Invoking Larsson-Sadakane on ranks time: 00:00:00
Sanity-checking and returning
Building samples
Reserving space for 12 sample suffixes
Generating random suffixes
QSorting 12 sample offsets, eliminating duplicates
QSorting sample offsets, eliminating duplicates time: 00:00:00
Multikey QSorting 12 samples
(Using difference cover)
Multikey QSorting samples time: 00:00:00
Calculating bucket sizes
Binary sorting into buckets
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Binary sorting into buckets time: 00:00:00
Splitting and merging
Splitting and merging time: 00:00:00
Split 1, merged 6; iterating...
Binary sorting into buckets
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Binary sorting into buckets time: 00:00:00
Splitting and merging
Splitting and merging time: 00:00:00
Split 1, merged 0; iterating...
Binary sorting into buckets
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Binary sorting into buckets time: 00:00:00
Splitting and merging
Splitting and merging time: 00:00:00
Avg bucket size: 39553.4 (target: 59330)
Converting suffix-array elements to index image
Allocating ftab, absorbFtab
Entering Ebwt loop
Getting block 1 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 38593
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 38594
Getting block 2 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 53650
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 53651
Getting block 3 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 45711
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 45712
Getting block 4 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 18201
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 18202
Getting block 5 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 56112
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 56113
Getting block 6 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 24508
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 24509
Getting block 7 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 52958
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 52959
Getting block 8 of 8
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 26694
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 26695
Exited Ebwt loop
fchr[A]: 0
fchr[C]: 117471
fchr[G]: 165360
fchr[T]: 220740
fchr[$]: 316434
Exiting Ebwt::buildToDisk()
Returning from initFromVector
Wrote 4302635 bytes to primary EBWT file: rotavirus.1.ebwt
Wrote 39560 bytes to secondary EBWT file: rotavirus.2.ebwt
Re-opening _in1 and _in2 as input streams
Returning from Ebwt constructor
Headers:
len: 316434
bwtLen: 316435
sz: 79109
bwtSz: 79109
lineRate: 6
linesPerSide: 1
offRate: 5
offMask: 0xffffffe0
isaRate: -1
isaMask: 0xffffffff
ftabChars: 10
eftabLen: 20
eftabSz: 80
ftabLen: 1048577
ftabSz: 4194308
offsLen: 9889
offsSz: 39556
isaLen: 0
isaSz: 0
lineSz: 64
sideSz: 64
sideBwtSz: 56
sideBwtLen: 224
numSidePairs: 707
numSides: 1414
numLines: 1414
ebwtTotLen: 90496
ebwtTotSz: 90496
Total time for call to driver() for forward index: 00:00:01
Reading reference sizes
Time reading reference sizes: 00:00:00
Calculating joined length
= 316434 (0 characters of padding)
Writing header
Reserving space for joined string
Joining reference sequences
Time to join reference sequences: 00:00:00
bmax according to bmaxDivN setting: 79108
Using parameters --bmax 59331 --dcv 1024
Doing ahead-of-time memory usage test
Passed! Constructing with these parameters: --bmax 59331 --dcv 1024
Constructing suffix-array element generator
Building DifferenceCoverSample
Building sPrime
Building sPrimeOrder
V-Sorting samples
V-Sorting samples time: 00:00:00
Allocating rank array
Ranking v-sort output
Ranking v-sort output time: 00:00:00
Invoking Larsson-Sadakane on ranks
Invoking Larsson-Sadakane on ranks time: 00:00:00
Sanity-checking and returning
Building samples
Reserving space for 12 sample suffixes
Generating random suffixes
QSorting 12 sample offsets, eliminating duplicates
QSorting sample offsets, eliminating duplicates time: 00:00:00
Multikey QSorting 12 samples
(Using difference cover)
Multikey QSorting samples time: 00:00:00
Calculating bucket sizes
Binary sorting into buckets
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Binary sorting into buckets time: 00:00:00
Splitting and merging
Splitting and merging time: 00:00:00
Split 1, merged 7; iterating...
Binary sorting into buckets
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Binary sorting into buckets time: 00:00:00
Splitting and merging
Splitting and merging time: 00:00:00
Avg bucket size: 45204 (target: 59330)
Converting suffix-array elements to index image
Allocating ftab, absorbFtab
Entering Ebwt loop
Getting block 1 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 55484
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 55485
Getting block 2 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 58214
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 58215
Getting block 3 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 58224
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 58225
Getting block 4 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 33666
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 33667
Getting block 5 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 48159
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 48160
Getting block 6 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 45260
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 45261
Getting block 7 of 7
Reserving size (59331) for bucket
Calculating Z arrays
Calculating Z arrays time: 00:00:00
Entering block accumulator loop:
10%
20%
30%
40%
50%
60%
70%
80%
90%
100%
Block accumulator loop time: 00:00:00
Sorting block of length 17421
(Using difference cover)
Sorting block time: 00:00:00
Returning block of 17422
Exited Ebwt loop
fchr[A]: 0
fchr[C]: 117471
fchr[G]: 165360
fchr[T]: 220740
fchr[$]: 316434
Exiting Ebwt::buildToDisk()
Returning from initFromVector
Wrote 4302635 bytes to primary EBWT file: rotavirus.rev.1.ebwt
Wrote 39560 bytes to secondary EBWT file: rotavirus.rev.2.ebwt
Re-opening _in1 and _in2 as input streams
Returning from Ebwt constructor
Headers:
len: 316434
bwtLen: 316435
sz: 79109
bwtSz: 79109
lineRate: 6
linesPerSide: 1
offRate: 5
offMask: 0xffffffe0
isaRate: -1
isaMask: 0xffffffff
ftabChars: 10
eftabLen: 20
eftabSz: 80
ftabLen: 1048577
ftabSz: 4194308
offsLen: 9889
offsSz: 39556
isaLen: 0
isaSz: 0
lineSz: 64
sideSz: 64
sideBwtSz: 56
sideBwtLen: 224
numSidePairs: 707
numSides: 1414
numLines: 1414
ebwtTotLen: 90496
ebwtTotSz: 90496
Total time for backward call to driver() for mirror index: 00:00:00


Searching


The short sequences found with the query "NCBI: ROTAVIRUS[organism] 10:150[SLEN] bovine" were downloaded to ./reads/short-rota.fasta.
Those sequences were then mapped on the indexed sequences:
./bowtie -f rotavirus reads/short-rota.fasta | java -jar ~/lindenb/build/verticalize.jar -n
Reported 14 alignments to 1 output stream(s)
>>1
$1 ? : gi|62781718|gb|AR643436.1| Sequence 6 from patent US 6867021
$2 ? : -
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 7504
$5 ? : GTGGCTTCCATTAGAAGCATG
$6 ? : IIIIIIIIIIIIIIIIIIIII
$7 ? : 1
<<1
>>2
$1 ? : gi|62781717|gb|AR643435.1| Sequence 5 from patent US 6867021
$2 ? : +
$3 ? : gi|66774335|gb|AH014892.1|SEG_DQ005104S Human rotavirus A strain DRC88 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 7230
$5 ? : ACCACCAAATATGACACCAGC
$6 ? : IIIIIIIIIIIIIIIIIIIII
$7 ? : 1
<<2
>>3
$1 ? : gi|56638523|gb|AR590758.1| Sequence 32 from patent US 6805867
$2 ? : -
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 5633
$5 ? : GGATGGCCAACAGGATC
$6 ? : IIIIIIIIIIIIIIIII
$7 ? : 1
$8 ? : 2:T>A,5:T>A
<<3
>>4
$1 ? : gi|56638522|gb|AR590757.1| Sequence 31 from patent US 6805867
$2 ? : +
$3 ? : gi|5706619|gb|AF106283.1|AF106283 Human rotavirus G2 strain TA25 VP7 protein (VP7) gene, complete cds
$4 ? : 22
$5 ? : GTATGGTATTGAATATACCAC
$6 ? : IIIIIIIIIIIIIIIIIIIII
$7 ? : 16
<<4
>>5
$1 ? : gi|56638509|gb|AR590744.1| Sequence 18 from patent US 6805867
$2 ? : -
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 7658
$5 ? : TAGTGAGAGGATGTGACC
$6 ? : IIIIIIIIIIIIIIIIII
$7 ? : 1
<<5
>>6
$1 ? : gi|56638508|gb|AR590743.1| Sequence 17 from patent US 6805867
$2 ? : +
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 6320
$5 ? : GGCTTTTAAACGAAGTC
$6 ? : IIIIIIIIIIIIIIIII
$7 ? : 1
<<6
>>7
$1 ? : gi|56638505|gb|AR590740.1| Sequence 14 from patent US 6805867
$2 ? : -
$3 ? : gi|66774335|gb|AH014892.1|SEG_DQ005104S Human rotavirus A strain DRC88 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 15298
$5 ? : TGTGGAGATATGACC
$6 ? : IIIIIIIIIIIIIII
$7 ? : 1
<<7
>>8
$1 ? : gi|56638504|gb|AR590739.1| Sequence 13 from patent US 6805867
$2 ? : +
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 12626
$5 ? : GGCTATTAAAGGT
$6 ? : IIIIIIIIIIIII
$7 ? : 1
$8 ? : 12:C>T
<<8
>>9
$1 ? : gi|10003339|gb|AR076593.1|AR076593 Sequence 32 from patent US 5959093
$2 ? : -
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 5633
$5 ? : GGATGGCCAACAGGATC
$6 ? : IIIIIIIIIIIIIIIII
$7 ? : 1
$8 ? : 2:T>A,5:T>A
<<9
>>10
$1 ? : gi|10003338|gb|AR076592.1|AR076592 Sequence 31 from patent US 5959093
$2 ? : +
$3 ? : gi|5706619|gb|AF106283.1|AF106283 Human rotavirus G2 strain TA25 VP7 protein (VP7) gene, complete cds
$4 ? : 22
$5 ? : GTATGGTATTGAATATACCAC
$6 ? : IIIIIIIIIIIIIIIIIIIII
$7 ? : 16
<<10
>>11
$1 ? : gi|10003325|gb|AR076579.1|AR076579 Sequence 18 from patent US 5959093
$2 ? : -
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 7658
$5 ? : TAGTGAGAGGATGTGACC
$6 ? : IIIIIIIIIIIIIIIIII
$7 ? : 1
<<11
>>12
$1 ? : gi|10003324|gb|AR076578.1|AR076578 Sequence 17 from patent US 5959093
$2 ? : +
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 6320
$5 ? : GGCTTTTAAACGAAGTC
$6 ? : IIIIIIIIIIIIIIIII
$7 ? : 1
<<12
>>13
$1 ? : gi|10003321|gb|AR076575.1|AR076575 Sequence 14 from patent US 5959093
$2 ? : -
$3 ? : gi|66774335|gb|AH014892.1|SEG_DQ005104S Human rotavirus A strain DRC88 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 15298
$5 ? : TGTGGAGATATGACC
$6 ? : IIIIIIIIIIIIIII
$7 ? : 1
<<13
>>14
$1 ? : gi|10003320|gb|AR076574.1|AR076574 Sequence 13 from patent US 5959093
$2 ? : +
$3 ? : gi|78064499|gb|AH014893.1|SEG_DQ005115S Human rotavirus A strain DRC86 NSP5, NSP4, NSP3, NSP2, NSP1, VP7, VP6, VP4, VP3, VP2, and VP1 genes, complete cds
$4 ? : 12626
$5 ? : GGCTATTAAAGGT
$6 ? : IIIIIIIIIIIII
$7 ? : 1
$8 ? : 12:C>T
<<14
(the verticalize is a small utility I wrote, it transforms a horizontal output to a vertical one).


In the 9th result says : the reverse sequence of "gi|10003339" is GGATGGCCAACAGGATC . This sequence matches "gi|78064499" at position '5633' there are two mismatches:2:T>A,5:T>A. Control: running a blas2seq on those two sequences gives the following result:


>gb|AR076593.1|AR076593 Sequence 32 from patent US 5959093
Length=17

Score = 22.9 bits (24), Expect = 0.014
Identities = 15/17 (88%), Gaps = 0/17 (0%)
Strand=Plus/Minus


Query 5634 GGATGGCCAACTGGTTC 5650
||||||||||| || ||
Sbjct 17 GGATGGCCAACAGGATC 1


That's it
Pierre

03 September 2009

Building a naive Interactome Database with Hibernate.

This is my notebook for building a naive database of protein-protein interactions with Hibernate (a java object/relational persistence and query service).

Files and Directories


./project
./project/lib
./project/src
./project/src/hibernate.cfg.xml
./project/src/org
./project/src/org/lindenb
./project/src/org/lindenb/hbn01
./project/src/org/lindenb/hbn01/Journal.java
./project/src/org/lindenb/hbn01/mapping.hbm.xml
./project/src/org/lindenb/hbn01/Article.java
./project/src/org/lindenb/hbn01/Main.java
./project/src/org/lindenb/hbn01/PMID.java
./project/src/org/lindenb/hbn01/Interactor.java
./project/src/org/lindenb/hbn01/Complex.java
./project/src/org/lindenb/hbn01/Protein.java
./project/src/org/lindenb/hbn01/PMIDType.java
./project/src/log4j.properties
./project/build
./project/bin
./derby.log
./build
./build/db
./Makefile

The components / Java Classes


Interactor


An abstract class defining a protein or a complex: Just a name and an ID.
package org.lindenb.hbn01;

public class Interactor
implements java.io.Serializable
{
private Long id;
private String name;
protected Interactor()
{
}

protected Interactor(String name)
{
setName(name);
}

private void setId(Long id)
{
this.id=id;
}

public Long getId()
{
return this.id;
}
public String getName()
{
return this.name;
}
public void setName(String name)
{
this.name=name;
}

@Override
public boolean equals(Object o)
{
if(o==this) return true;
if(o==null || o.getClass()!=getClass()) return false;
return getId().equals(Interactor.class.cast(o).getId());
}

@Override
public String toString()
{
return getClass().getName()+":"+getName()+"("+getId()+")";
}
}

Protein

Protein is a concrete subclass of Interactor. This could be an Unigene entry.
package org.lindenb.hbn01;

public class Protein
extends Interactor
{
public Protein()
{
}

public Protein(String name)
{
super(name);
}

@Override
public String toString()
{
return "Protein:"+getName();
}
}

Complex

Complex is a concrete subclass of Interactor. It is a Set of Interactors. It also contains a Set of Articles holding the references for those interactions.
package org.lindenb.hbn01;
import java.util.*;

public class Complex
extends Interactor
{
private Set<Interactor> partners= new HashSet<Interactor>();
private Set<Article> articles= new HashSet<Article>();
public Complex()
{
}

public Complex(String name)
{
super(name);
}

public Set<Interactor> getPartners()
{
return this.partners;
}
public void setPartners(Set<Interactor> partners)
{
this.partners = partners;
}
public Set<Article> getArticles()
{
return this.articles;
}

public void setArticles(Set<Article> articles)
{
this.articles = articles;
}
@Override
public String toString()
{
String s="Complex:"+getName()+". ID:"+getId()+" interacts with";
for(Interactor i: getPartners())
{
s+=" "+i.getName();
}
return s;
}
}

Article

an Article is a reference to a paper in Pubmed. I wanted to use the custom dataType in hibernate, so I used the class PMID rather than an Integer. Each Article is linked to a Journal.
package org.lindenb.hbn01;

public class Article
implements java.io.Serializable
{
private PMID pmid;
private String title;
private Integer year;
private String doi;
private Journal journal;

public Article()
{
}

public Article(PMID pmid,Journal journal,Integer year,String title)
{
setPmid(pmid);
setJournal(journal);
setYear(year);
setTitle(title);
}

public Journal getJournal()
{
return journal;
}

public String getDoi()
{
return this.doi;
}

public void setDoi(String doi)
{
this.doi=doi;
}

public void setJournal(Journal journal)
{
this.journal=journal;
}

private void setPmid(PMID pmid)
{
this.pmid=pmid;
}
public PMID getPmid()
{
return this.pmid;
}
public void setTitle(String title)
{
this.title=title;
}
public String getTitle()
{
return this.title;
}
public void setYear(Integer year)
{
this.year=year;
}
public Integer getYear()
{
return this.year;
}

public String toString()
{
return "("+getYear()+")\""+getTitle()+"\"."+getJournal().getTitle();
}
}

PMID

A custom type holding a Pubmed identifier
package org.lindenb.hbn01;
import org.hibernate.*;
import java.io.Serializable;

public class PMID
implements java.io.Serializable
{
private long pmid;

public PMID(String pmid)
{
this(new Long(pmid));
}

public PMID(long pmid)
{
this.pmid=pmid;
}

public long value()
{
return this.pmid;
}

public int hashCode()
{
return 31+(int)this.pmid;
}

public boolean equals(Object o)
{
if(o==this) return true;
if(o==null || !(o instanceof PMID)) return false;
return PMID.class.cast(o).pmid==this.pmid;
}

public String toString()
{
return String.valueOf(this.pmid);
}
}

Journal

A Journal is a NLM-Id and a title
package org.lindenb.hbn01;

public class Journal
implements java.io.Serializable
{
private long nlmId;
private String title;


public Journal()
{
}

public Journal(long nlmId,String title)
{
setNlmId(nlmId);
setTitle(title);
}

private void setNlmId(long nlmId)
{
this.nlmId=nlmId;
}
public Long getNlmId()
{
return this.nlmId;
}
public void setTitle(String title)
{
this.title=title;
}
public String getTitle()
{
return this.title;
}

public String toString()
{
return getTitle()+"["+getNlmId()+"]";
}
}

Using a Custom Type

PMIDType implements EnhancedUserType. Hibernate will use this class to manage the class PMID (how to read/write it from/to the database).
package org.lindenb.hbn01;
import org.hibernate.*;
import org.hibernate.usertype.EnhancedUserType;
import java.sql.*;
import java.io.Serializable;

public class PMIDType
implements EnhancedUserType
{


public int[] sqlTypes() {
return new int[]{Types.INTEGER};
}


public Object assemble(Serializable cached,
Object owner)
throws HibernateException
{
return cached;
}

public Serializable disassemble(Object value)
throws HibernateException
{
return Serializable.class.cast(value);
}

public boolean isMutable() { return false;}
public Object deepCopy(Object value)
{
return value;
}
public boolean equals(Object a, Object b)
{
return a==null?b==null:a.equals(b);
}

public int hashCode(Object x) throws HibernateException
{
return x==null?0:x.hashCode();
}

public Object nullSafeGet(ResultSet rs,
String[] names,
Object owner)
throws HibernateException,
SQLException
{
Object o = rs.getObject( names[0] );
if(rs.wasNull()) return null;
if(o instanceof Number)
{
return new PMID(Number.class.cast(o).longValue());
}
else if(o instanceof String)
{
return new PMID(String.class.cast(o));
}
throw new IllegalArgumentException("Bad class "+o.getClass());
}

public void nullSafeSet(PreparedStatement st,
Object value,
int index)
throws HibernateException, SQLException
{
if(value==null)
{
st.setNull( index, Types.INTEGER );
}
else
{
st.setLong(index,PMID.class.cast(value).value());
}
}


public Object replace(Object original,
Object target,
Object owner)
throws HibernateException
{
return original;
}

public Class<?> returnedClass()
{
return PMID.class;
}


public Object fromXMLString(String xmlValue)
{
return xmlValue==null? null: new PMID(new Long(xmlValue));
}
public String objectToSQLString(Object value)
{
return value==null? null: String.valueOf(PMID.class.cast(value).
value());
}
public String toXMLString(Object value)
{
return value==null? null: String.valueOf(PMID.class.cast(value).
value());
}
}

The mapping file

The file mapping.hbm.xml tells hibernate how the classes are linked to each others.
<!DOCTYPE hibernate-mapping PUBLIC
"-//Hibernate/Hibernate Mapping DTD 3.0//EN"
"http://hibernate.sourceforge.net/hibernate-mapping-3.0.dtd">
<hibernate-mapping package="org.lindenb.hbn01" default-cascade="none" default-access="property" default-lazy="true" auto-import="true">

<class name="org.lindenb.hbn01.Article" table="Article" mutable="true" polymorphism="implicit" dynamic-update="false" dynamic-insert="false" select-before-update="false" optimistic-lock="version">
<meta attribute="class-description" inherit="true">A pubmed Article</meta>
<id name="pmid" column="pmid" type="org.lindenb.hbn01.PMIDType">
<meta attribute="field-description" inherit="true">Pubmed Identifier</meta>
<generator class="assigned"/>
</id>
<property name="title" not-null="true" unique="false" optimistic-lock="true" lazy="false" generated="never"/>
<property name="doi" unique="true" type="string" optimistic-lock="true" lazy="false" generated="never"/>
<property name="year" column="yearDate" type="integer" not-null="true" unique="false" optimistic-lock="true" lazy="false" generated="never"/>
<many-to-one name="journal" column="nlmId" not-null="true" unique="false" update="true" insert="true" optimistic-lock="true" not-found="exception" embed-xml="true"/>
</class>

<class name="Journal" mutable="true" polymorphism="implicit" dynamic-update="false" dynamic-insert="false" select-before-update="false" optimistic-lock="version">
<id name="nlmId" column="nlmId" type="long">
<generator class="assigned"/>
</id>
<property name="title" not-null="true" unique="false" optimistic-lock="true" lazy="false" generated="never"/>
</class>

<class name="Interactor" mutable="true" polymorphism="implicit" dynamic-update="false" dynamic-insert="false" select-before-update="false" optimistic-lock="version">
<id name="id" type="long">
<generator class="native"/>
</id>
<property name="name" not-null="true" unique="false" optimistic-lock="true" lazy="false" generated="never"/>

<joined-subclass name="Protein" dynamic-update="false" dynamic-insert="false" select-before-update="false">
<key column="interactorId" on-delete="noaction"/>
</joined-subclass>

<joined-subclass name="Complex" dynamic-update="false" dynamic-insert="false" select-before-update="false">
<key column="interactorId" on-delete="noaction"/>

<set name="partners" table="interactions" sort="unsorted" inverse="false" mutable="true" optimistic-lock="true" embed-xml="true">
<key column="complex_id" on-delete="noaction"/>
<many-to-many column="interactor_id" class="Interactor" embed-xml="true" not-found="exception" unique="false"/>
</set>

</joined-subclass>
</class>

</hibernate-mapping>

Configuring Hibernate


The file hibernate.cfg.xml describes the database we are using for persisting all the entities (driver, uri, login, password...). Here I've used JavaDB/Derby (Note: the 'unique' directive was ignored by Derby (?) ).
<!DOCTYPE hibernate-configuration PUBLIC
"-//Hibernate/Hibernate Configuration DTD 3.0//EN"
"http://hibernate.sourceforge.net/hibernate-configuration-3.0.dtd">
<hibernate-configuration>

<session-factory>

<!-- Database connection settings -->
<property name="connection.driver_class">org.apache.derby.jdbc.EmbeddedDriver</property>
<property name="connection.url">jdbc:derby:build/db/derby/hibernate;create=true</property>
<property name="connection.username">sa</property>
<property name="connection.password"/>

<!-- JDBC connection pool (use the built-in) -->
<property name="connection.pool_size">1</property>

<!-- SQL dialect -->
<property name="dialect">org.hibernate.dialect.DerbyDialect</property>

<!-- Enable Hibernate's automatic session context management -->
<property name="current_session_context_class">thread</property>

<!-- Disable the second-level cache -->
<property name="cache.provider_class">org.hibernate.cache.NoCacheProvider</property>

<!-- Echo all executed SQL to stdout -->
<property name="show_sql">true</property>

<!-- Drop and re-create the database schema on startup -->
<property name="hbm2ddl.auto">create</property>

<mapping resource="org/lindenb/hbn01/mapping.hbm.xml"/>

</session-factory>

</hibernate-configuration>

Running


Building a Session factory

sessionFactory = new Configuration().configure().buildSessionFactory();

Creating an Interacome

Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
Journal journal = new Journal(1L,"PNAS");
session.save(journal);
Article article= new Article(new PMID(12234),journal,1988,"Article title 1");
article.setDoi("1");
session.save(article);
article= new Article(new PMID(456789),journal,1989,"Article title 2");
article.setDoi("2");
session.save(article);


Protein prot1= new Protein("prot1");
session.save(prot1);
Protein prot2= new Protein("prot2");
session.save(prot2);
Protein prot3= new Protein("prot3");
session.save(prot3);
Complex c1= new Complex("cplx1");
c1.getPartners().add(prot1);
c1.getPartners().add(prot2);
session.save(c1);
Complex c2= new Complex("cplx2");
c2.getPartners().add(prot1);
c2.getPartners().add(c1);
c2.getArticles().add(article);
session.save(c2);

session.getTransaction().commit();

Querying

Listing the Journals
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
List list = session.createQuery("from Journal").list();

for(Object o:list)
{
System.out.println(o);
}
session.getTransaction().commit();

Listing the Articles
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
List list = session.createQuery("from Article").list();

for(Object o:list)
{
System.out.println(o);
}
session.getTransaction().commit();
Listing the Interactors
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
List list = session.createQuery("from Interactor").list();

for(Object o:list)
{
System.out.println("\n\n###\t"+o+"\n\n");
}
session.getTransaction().commit();

Full code

package org.lindenb.hbn01;

import org.hibernate.*;
import org.hibernate.cfg.*;
import java.util.*;

public class Main
{
private static final SessionFactory sessionFactory;

static
{
try
{
sessionFactory = new Configuration().configure().buildSessionFactory();
}
catch(Throwable err)
{
err.printStackTrace();
throw new ExceptionInInitializerError(err);
}
}

public static SessionFactory getSessionFactory()
{
return Main.sessionFactory;
}

private void listJournals()
{
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
List list = session.createQuery("from Journal").list();

for(Object o:list)
{
System.out.println(o);
}
session.getTransaction().commit();
}

private void listArticles()
{
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
List list = session.createQuery("from Article").list();

for(Object o:list)
{
System.out.println(o);
}
session.getTransaction().commit();
}

private void listInteractors()
{
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
List list = session.createQuery("from Interactor").list();

for(Object o:list)
{
System.out.println("\n\n###\t"+o+"\n\n");
}
session.getTransaction().commit();
}

public void run()
{
Session session= getSessionFactory().getCurrentSession();
session.beginTransaction();
Journal journal = new Journal(1L,"PNAS");
session.save(journal);
Article article= new Article(new PMID(12234),journal,1988,"Article title 1");
article.setDoi("1");
session.save(article);
article= new Article(new PMID(456789),journal,1989,"Article title 2");
article.setDoi("2");
session.save(article);


Protein prot1= new Protein("prot1");
session.save(prot1);
Protein prot2= new Protein("prot2");
session.save(prot2);
Protein prot3= new Protein("prot3");
session.save(prot3);
Complex c1= new Complex("cplx1");
c1.getPartners().add(prot1);
c1.getPartners().add(prot2);
session.save(c1);
Complex c2= new Complex("cplx2");
c2.getPartners().add(prot1);
c2.getPartners().add(c1);
c2.getArticles().add(article);
session.save(c2);

session.getTransaction().commit();

listJournals();
listArticles();
listInteractors();
}

public static void main(String args[])
{
try
{
Main app= new Main();
app.run();
}
catch(Throwable err)
{
err.printStackTrace();
}
finally
{
if(Main.sessionFactory!=null) Main.sessionFactory.close();
}
System.out.println("Done.");
}
}

Compiling

LIB=${HIBERNATE_HOME}/lib
LIBS=${LIB}/antlr-2.7.6.jar:${LIB}/cglib-2.1.3.jar:${LIB}/asm.jar:${LIB}/asm-attrs.jar:${LIB}/commons-collections-2.1.1.jar:${LIB}/commons-logging-1.0.4.jar:${HIBERNATE_HOME}/hibernate3.jar:${LIB}/jta.jar:${LIB}/dom4j-1.6.1.jar:${LIB}/log4j-1.2.11.jar:${DERBY_HOME}/derby.jar
test:
cp -r project/src/* project/build
javac -cp ${LIBS} -d project/build -sourcepath project/build project/build/org/lindenb/hbn01/*.java
jar cvf project/bin/project.jar -C project/build .
java -cp ${LIBS}:project/bin/project.jar org.lindenb.hbn01.Main

Output

21:49:15,396 INFO Environment:514 - Hibernate 3.2.6
21:49:15,402 INFO Environment:547 - hibernate.properties not found
21:49:15,405 INFO Environment:681 - Bytecode provider name : cglib
21:49:15,409 INFO Environment:598 - using JDK 1.4 java.sql.Timestamp handling
21:49:15,457 INFO Configuration:1432 - configuring from resource: /hibernate.cfg.xml
21:49:15,458 INFO Configuration:1409 - Configuration resource: /hibernate.cfg.xml
21:49:15,546 INFO Configuration:559 - Reading mappings from resource : org/lindenb/hbn01/mapping.hbm.xml
21:49:15,682 INFO HbmBinder:300 - Mapping class: org.lindenb.hbn01.Article -> Article
21:49:15,751 INFO HbmBinder:300 - Mapping class: org.lindenb.hbn01.Journal -> Journal
21:49:15,752 INFO HbmBinder:300 - Mapping class: org.lindenb.hbn01.Interactor -> Interactor
21:49:15,778 INFO HbmBinder:873 - Mapping joined-subclass: org.lindenb.hbn01.Protein -> Protein
21:49:15,780 INFO HbmBinder:873 - Mapping joined-subclass: org.lindenb.hbn01.Complex -> Complex
21:49:15,781 INFO HbmBinder:1419 - Mapping collection: org.lindenb.hbn01.Complex.partners -> interactions
21:49:15,783 INFO Configuration:1547 - Configured SessionFactory: null
21:49:15,802 INFO DriverManagerConnectionProvider:41 - Using Hibernate built-in connection pool (not for production use!)
21:49:15,803 INFO DriverManagerConnectionProvider:42 - Hibernate connection pool size: 1
21:49:15,803 INFO DriverManagerConnectionProvider:45 - autocommit mode: false
21:49:16,036 INFO DriverManagerConnectionProvider:80 - using driver: org.apache.derby.jdbc.EmbeddedDriver at URL: jdbc:derby:build/db/derby/hibernate;create=true
21:49:16,036 INFO DriverManagerConnectionProvider:86 - connection properties: {user=sa, password=****}
21:49:18,144 INFO SettingsFactory:89 - RDBMS: Apache Derby, version: 10.2.2.1 - (538595)
21:49:18,145 INFO SettingsFactory:90 - JDBC driver: Apache Derby Embedded JDBC Driver, version: 10.2.2.1 - (538595)
21:49:18,158 INFO Dialect:152 - Using dialect: org.hibernate.dialect.DerbyDialect
21:49:18,165 INFO TransactionFactoryFactory:31 - Using default transaction strategy (direct JDBC transactions)
21:49:18,167 INFO TransactionManagerLookupFactory:33 - No TransactionManagerLookup configured (in JTA environment, use of read-write or transactional second-level cache is not recommended)
21:49:18,167 INFO SettingsFactory:143 - Automatic flush during beforeCompletion(): disabled
21:49:18,167 INFO SettingsFactory:147 - Automatic session close at end of transaction: disabled
21:49:18,168 INFO SettingsFactory:162 - Scrollable result sets: enabled
21:49:18,168 INFO SettingsFactory:170 - JDBC3 getGeneratedKeys(): disabled
21:49:18,169 INFO SettingsFactory:178 - Connection release mode: auto
21:49:18,169 INFO SettingsFactory:205 - Default batch fetch size: 1
21:49:18,170 INFO SettingsFactory:209 - Generate SQL with comments: disabled
21:49:18,170 INFO SettingsFactory:213 - Order SQL updates by primary key: disabled
21:49:18,170 INFO SettingsFactory:217 - Order SQL inserts for batching: disabled
21:49:18,170 INFO SettingsFactory:386 - Query translator: org.hibernate.hql.ast.ASTQueryTranslatorFactory
21:49:18,172 INFO ASTQueryTranslatorFactory:24 - Using ASTQueryTranslatorFactory
21:49:18,173 INFO SettingsFactory:225 - Query language substitutions: {}
21:49:18,173 INFO SettingsFactory:230 - JPA-QL strict compliance: disabled
21:49:18,173 INFO SettingsFactory:235 - Second-level cache: enabled
21:49:18,173 INFO SettingsFactory:239 - Query cache: disabled
21:49:18,174 INFO SettingsFactory:373 - Cache provider: org.hibernate.cache.NoCacheProvider
21:49:18,174 INFO SettingsFactory:254 - Optimize cache for minimal puts: disabled
21:49:18,174 INFO SettingsFactory:263 - Structured second-level cache entries: disabled
21:49:18,178 INFO SettingsFactory:283 - Echoing all SQL to stdout
21:49:18,178 INFO SettingsFactory:290 - Statistics: disabled
21:49:18,178 INFO SettingsFactory:294 - Deleted entity synthetic identifier rollback: disabled
21:49:18,178 INFO SettingsFactory:309 - Default entity-mode: pojo
21:49:18,179 INFO SettingsFactory:313 - Named query checking : enabled
21:49:18,206 INFO SessionFactoryImpl:161 - building session factory
21:49:18,472 INFO SessionFactoryObjectFactory:82 - Not binding factory to JNDI, no JNDI name configured
21:49:18,477 INFO SchemaExport:154 - Running hbm2ddl schema export
21:49:18,477 DEBUG SchemaExport:170 - import file not found: /import.sql
21:49:18,478 INFO SchemaExport:179 - exporting generated schema to database
21:49:18,482 DEBUG SchemaExport:303 - alter table Article drop constraint FK379164D684F84236
21:49:18,754 DEBUG SchemaExport:303 - alter table Complex drop constraint FK9BDFFC90D86C24B8
21:49:18,798 DEBUG SchemaExport:303 - alter table Protein drop constraint FK50CD6F63D86C24B8
21:49:18,834 DEBUG SchemaExport:303 - alter table interactions drop constraint FK4F6EF4A127BFBBC5
21:49:18,903 DEBUG SchemaExport:303 - alter table interactions drop constraint FK4F6EF4A1EBB9AEEF
21:49:19,029 DEBUG SchemaExport:303 - drop table Article
21:49:19,173 DEBUG SchemaExport:303 - drop table Complex
21:49:19,276 DEBUG SchemaExport:303 - drop table Interactor
21:49:19,410 DEBUG SchemaExport:303 - drop table Journal
21:49:19,537 DEBUG SchemaExport:303 - drop table Protein
21:49:19,654 DEBUG SchemaExport:303 - drop table interactions
21:49:19,814 DEBUG SchemaExport:303 - drop table hibernate_unique_key
21:49:19,896 DEBUG SchemaExport:303 - create table Article (pmid integer not null, title varchar(255) not null, doi varchar(255), yearDate integer not null, nlmId bigint not null, primary key (pmid))
21:49:20,063 DEBUG SchemaExport:303 - create table Complex (interactorId bigint not null, primary key (interactorId))
21:49:20,224 DEBUG SchemaExport:303 - create table Interactor (id bigint not null, name varchar(255) not null, primary key (id))
21:49:20,359 DEBUG SchemaExport:303 - create table Journal (nlmId bigint not null, title varchar(255) not null, primary key (nlmId))
21:49:20,580 DEBUG SchemaExport:303 - create table Protein (interactorId bigint not null, primary key (interactorId))
21:49:20,725 DEBUG SchemaExport:303 - create table interactions (complex_id bigint not null, interactor_id bigint not null, primary key (complex_id, interactor_id))
21:49:20,858 DEBUG SchemaExport:303 - alter table Article add constraint FK379164D684F84236 foreign key (nlmId) references Journal
21:49:20,997 DEBUG SchemaExport:303 - alter table Complex add constraint FK9BDFFC90D86C24B8 foreign key (interactorId) references Interactor
21:49:21,055 DEBUG SchemaExport:303 - alter table Protein add constraint FK50CD6F63D86C24B8 foreign key (interactorId) references Interactor
21:49:21,086 DEBUG SchemaExport:303 - alter table interactions add constraint FK4F6EF4A127BFBBC5 foreign key (interactor_id) references Interactor
21:49:21,202 DEBUG SchemaExport:303 - alter table interactions add constraint FK4F6EF4A1EBB9AEEF foreign key (complex_id) references Complex
21:49:21,331 DEBUG SchemaExport:303 - create table hibernate_unique_key ( next_hi integer )
21:49:21,376 DEBUG SchemaExport:303 - insert into hibernate_unique_key values ( 0 )
21:49:21,500 INFO SchemaExport:196 - schema export complete
21:49:21,501 WARN JDBCExceptionReporter:54 - SQL Warning: 10000, SQLState: 01J01
21:49:21,501 WARN JDBCExceptionReporter:55 - Database 'build/db/derby/hibernate' not created, connection made to existing database instead.
21:49:21,724 WARN JDBCExceptionReporter:54 - SQL Warning: 10000, SQLState: 01J01
21:49:21,724 WARN JDBCExceptionReporter:55 - Database 'build/db/derby/hibernate' not created, connection made to existing database instead.
Hibernate: insert into Journal (title, nlmId) values (?, ?)
Hibernate: insert into Article (title, doi, yearDate, nlmId, pmid) values (?, ?, ?, ?, ?)
Hibernate: insert into Article (title, doi, yearDate, nlmId, pmid) values (?, ?, ?, ?, ?)
Hibernate: insert into Interactor (name, id) values (?, ?)
Hibernate: insert into Protein (interactorId) values (?)
Hibernate: insert into Interactor (name, id) values (?, ?)
Hibernate: insert into Protein (interactorId) values (?)
Hibernate: insert into Interactor (name, id) values (?, ?)
Hibernate: insert into Protein (interactorId) values (?)
Hibernate: insert into Interactor (name, id) values (?, ?)
Hibernate: insert into Complex (interactorId) values (?)
Hibernate: insert into Interactor (name, id) values (?, ?)
Hibernate: insert into Complex (interactorId) values (?)
Hibernate: insert into interactions (complex_id, interactor_id) values (?, ?)
Hibernate: insert into interactions (complex_id, interactor_id) values (?, ?)
Hibernate: insert into interactions (complex_id, interactor_id) values (?, ?)
Hibernate: insert into interactions (complex_id, interactor_id) values (?, ?)
Hibernate: select journal0_.nlmId as nlmId1_, journal0_.title as title1_ from Journal journal0_
PNAS[1]
Hibernate: select article0_.pmid as pmid0_, article0_.title as title0_, article0_.doi as doi0_, article0_.yearDate as yearDate0_, article0_.nlmId as nlmId0_ from Article article0_
Hibernate: select journal0_.nlmId as nlmId1_0_, journal0_.title as title1_0_ from Journal journal0_ where journal0_.nlmId=?
(1988)"Article title 1".PNAS
(1989)"Article title 2".PNAS

Hibernate: select interactor0_.id as id2_, interactor0_.name as name2_, case when interactor0_1_.interactorId is not null then 1 when interactor0_2_.interactorId is not null then 2 when interactor0_.id is not null then 0 else -1 end as clazz_ from Interactor interactor0_ left outer join Protein interactor0_1_ on interactor0_.id=interactor0_1_.interactorId left outer join Complex interactor0_2_ on interactor0_.id=interactor0_2_.interactorId


### Protein:prot1




### Protein:prot2




### Protein:prot3


Hibernate: select partners0_.complex_id as complex1_1_, partners0_.interactor_id as interactor2_1_, interactor1_.id as id2_0_, interactor1_.name as name2_0_, case when interactor1_1_.interactorId is not null then 1 when interactor1_2_.interactorId is not null then 2 when interactor1_.id is not null then 0 else -1 end as clazz_0_ from interactions partners0_ left outer join Interactor interactor1_ on partners0_.interactor_id=interactor1_.id left outer join Protein interactor1_1_ on interactor1_.id=interactor1_1_.interactorId left outer join Complex interactor1_2_ on interactor1_.id=interactor1_2_.interactorId where partners0_.complex_id=?


### Complex:cplx1. ID:4 interacts with prot2 prot1


Hibernate: select partners0_.complex_id as complex1_1_, partners0_.interactor_id as interactor2_1_, interactor1_.id as id2_0_, interactor1_.name as name2_0_, case when interactor1_1_.interactorId is not null then 1 when interactor1_2_.interactorId is not null then 2 when interactor1_.id is not null then 0 else -1 end as clazz_0_ from interactions partners0_ left outer join Interactor interactor1_ on partners0_.interactor_id=interactor1_.id left outer join Protein interactor1_1_ on interactor1_.id=interactor1_1_.interactorId left outer join Complex interactor1_2_ on interactor1_.id=interactor1_2_.interactorId where partners0_.complex_id=?


### Complex:cplx2. ID:5 interacts with cplx1 prot1


21:49:22,076 INFO SessionFactoryImpl:769 - closing
21:49:22,076 INFO DriverManagerConnectionProvider:147 - cleaning up connection pool: jdbc:derby:build/db/derby/hibernate;create=true
Done


That's it!
Pierre

Generating a C Pull Parser for dbSNP with XSLT



I've used the XSD schema describing dbSNP to generate a C "Pull parser" reading the content of the dbSNP XML files. To transform the schema into a C code I wrote the following XSLT stylesheet:. This stylesheet was specifically developed for dbSNP so it might not handle a more complicated schema (for example a schema that would use <xsd:elementType> ). Basically the C code generated is a scaffold for a Pull Parser using the libxml2 library. For example, here is a simplified snippet of code handling the tag <Assembly/>
/** A collection of genome sequence records (curated gene regions (NG's),
contigs (NWNT's) and chromosomes (NC/AC's) produced by a genome
sequence project. Structure is populated from ContigInfo tables. */

static int processAssembly(StatePtr state)
{
int returnValue=EXIT_SUCCESS;
int success;
int nodeType;
const int isEmptyElement= xmlTextReaderIsEmptyElement(state -> reader);

/** Name of the group(s) or organization(s) that generated the assembly */
xmlChar* assemblySourceAttr=NULL;

//(...) declare other attributes

assemblySourceAttr= xmlTextReaderGetAttribute(
state->reader,
BAD_CAST "assemblySource"
);

//(...) other attributes

if(!isEmptyElement)
{
success = xmlTextReaderRead( state -> reader );
if(!success)
{
fprintf( state->error,"In Assembly I/O Error. xmlTextReaderRead returned \n");
returnValue = EXIT_FAILURE;
goto cleanup;
}
nodeType = xmlTextReaderNodeType( state -> reader );


/* process childNode <Component/> */

while(nodeType == XML_READER_TYPE_ELEMENT)
{
if(xmlStrcmp(
xmlTextReaderConstName(state -> reader),
BAD_CAST "Component"
)!=0)
{
break;
}

if(processComponent(state)!=EXIT_SUCCESS)
{
returnValue = EXIT_FAILURE;
goto cleanup;
}

/* read next event */
success= xmlTextReaderRead(state->reader);
if(!success)
{
returnValue = EXIT_FAILURE;
fprintf( state->error,"In Assembly/Component I/O Error.\n");
goto cleanup;
}
nodeType=xmlTextReaderNodeType(state->reader);
}

/* process childNode <SnpStat/> */
(...)

}//end of if(!isEmptyElement)

cleanup:

//free attributes
if(assemblySourceAttr!=NULL)
{
xmlFree(assemblySourceAttr);
}
//(...) other attributes
return returnValue;
}
Using this prototype I was able to quickly write a fast parser , for example echoing a JSON description of the SNPs of the Human Mitochondrial Genome (time: 0.11user 0.01system 0:00.19elapsed 66%CPU).
[
{
"rsId":8896,
"seq5":"GGTGTTGGTTCTCTTAATCTTTAACTTAAAAGGTTAATGCTAAGTTAGCTTTACAGTGGGCTCTAGAGGGGG
TAGAGGGGGTG",
"observed":"C/T",
"seq3":"TATAGGGTAAATACGGGCCCTATTTCAAAGATTTTTAGGGGAATTAATTCTAGGACGATGGGCATGAAACTGTGGTTTGCTCCACAGATTTCAGAGCATT"
}
,
{
"rsId":8936,
"seq5":"ACTACGGCGGACTAATCTTCAACTCCTACATACTTCCCCCATTATTCCTAGAACCAGGCGACCTGCGACTCCTTGACGTTGACAATCGAGTAGTACTCCCGATTGAAGCCCCCATTCGTATAATAATTACATCACAAGACGTCTTGCACTCATGAGCTGTCCCCACATTAGGCTTAAAAACAGATGCAATTCCCGGACGT",
"observed":"A/C/T",
"seq3":"TAAACCAAACCACTTTCACCGCTACACGACCGGGGGTATACTACGGTCAATGCTCTGAAATCTGTGGAGCAAACCACAGTTTCATGCCCATCGTCCTAGAATTAATTCCCCTAAAAATCTTTGAAATAGGGCCCGTATTTACCCTATAGCACCCCCTCTACCCCCTCTAGAGCCCACTGTAAAGCTAACTTAGCATTAAC"
}
(...)
{
"rsId":72619366,
"seq5":"TGCTTACAAGCAAGTACAGCAATCAACCTTCAACTATCACACATCAACTGCAACTCCAAAGCCACCCCTCACCCACTAGGATACCAACAAACCTACCCAC",
"observed":"C/T",
"seq3":"CTTAACAGTACATAGTACATAAAGCCATTTACCGTACATAGCACATTACAGTCAAATCCCTTCTCGTCCCCATGGATGACCCCCCTCAGATAGGGGTCCC"
}
]



That's it
Pierre

01 September 2009

First steps with BerkeleyDB-XML. My notebook.

Berkeley DB XML is an embeddable XML database engine that provides support for XQuery access to documents stored in containers and indexed based on their content. Oracle Berkeley DB XML is built on top of Oracle Berkeley DB. Berkeley DB XML is available at : http://www.oracle.com/database/berkeley-db/xml/index.html. The distribution comes with a shell command.

pierre@linux-zfgk:.../dbxml-2.4.16> ./install/bin/dbxml

Creating a DataStore

dbxml> createContainer dbsnp.dbxml
Creating node storage container

Creating a set of XML documents describing some SNPs

dbxml> putDocument snp1 '<snp id="25">
<name>rs25</name>
<class>snp</class>
<het>0.5</het>
<observed>A/G</observed>
<mapping>
<location build="36_3" label="CRA_TCAGchr7v2" chrom="7" position="11637562"/>
<location build="36_3" label="Celera" chrom="7" position="11558958"/>
<location build="36_3" label="HuRef" chrom="7" position="11442496"/>
<location build="36_3" label="reference" chrom="7" position="11550666"/>
</mapping>
</snp>' s
Document added, name = snp1

dbxml> putDocument snp2 '<snp id="26">
<name>rs26</name>
<class>mixed</class>
<het>0</het>
<observed>-/A/G</observed>
<mapping>
<location build="36_3" label="CRA_TCAGchr7v2" chrom="7" position="11636891"/>
<location build="36_3" label="Celera" chrom="7" position="11558287"/>
<location build="36_3" label="HuRef" chrom="7" position="11441825"/>
<location build="36_3" label="reference" chrom="7" position="11549995"/>
</mapping>
</snp>' s
Document added, name = snp2

dbxml> putDocument snp3 '<snp id="27">
<name>rs27</name>
<class>snp</class>
<het>0.44</het>
<observed>C/G</observed>
<mapping>
<location build="36_3" label="CRA_TCAGchr7v2" chrom="7" position="11636645"/>
<location build="36_3" label="Celera" chrom="7" position="11558041"/>
<location build="36_3" label="HuRef" chrom="7" position="11441579"/>
<location build="36_3" label="reference" chrom="7" position="11549749"/>
</mapping>
</snp>' s
Document added, name = snp3


dbxml> putDocument snp4 '<snp id="300">
<name>rs300</name>
<class>snp</class>
<het>0.01</het>
<observed>A/G</observed>
<mapping>
<location build="36_3" label="Celera" chrom="8" position="18779978"/>
<location build="36_3" label="HuRef" chrom="8" position="18357119"/>
<location build="36_3" label="reference" chrom="8" position="19861166"/>
</mapping>
</snp>'
Document added, name = snp4

dbxml> putDocument snp5 '<snp id="600">
<name>rs600</name>
<class>snp</class>
<het>0.27</het>
<observed>C/G</observed>
<mapping>
<location build="36_3" label="Celera" chrom="X" position="148984992"/>
<location build="36_3" label="HuRef" chrom="X" position="137590170"/>
<location build="36_3" label="reference" chrom="X" position="148444179"/>
<location build="36_3" label="reference" chrom="X" position="148843642"/>
</mapping>
</snp>' s
Document added, name = snp5

dbxml> putDocument snp6 '<snp id="800">
<name>rs800</name>
<class>snp</class>
<het/>
<observed>C/G/T</observed>
<mapping>
<location build="36_3" label="Celera" chrom="22" position="18365004"/>
<location build="36_3" label="HuRef" chrom="22" position="17518747"/>
<location build="36_3" label="reference" chrom="22" position="32892297"/>
</mapping>
</snp>' s
Document added, name = snp6

Getting help

dbxml> help

Command Summary
---------------

# - Comment. Does nothing
abort - Aborts the current transaction
addAlias - Add an alias to the default container
addIndex - Add an index to the default container
append - Append to nodes specified in the query expression
commit - Commits the current transaction, and starts a new one
compactContainer - Compact a container to shrink it's size
contextQuery - Execute query expression using the last results as the context item
cquery - Execute an expression in the context of the default container
createContainer - Creates a new container, which becomes the default container
debug - Debug command -- internal use only
delIndex - Delete an index from the default container
echo - Echo to output
getDocuments - Gets document(s) by name from default container
getMetaData - Get a metadata item from the named document
help - Print help information. Use 'help commandName' for extended help
info - Get info on default container
insertAfter - Insert new content after nodes selected by the query expression
insertBefore - Insert new content before nodes selected by the query expression
listIndexes - List all indexes in the default container
lookupEdgeIndex - Performs an edge index lookup in the default container
lookupIndex - Performs an index lookup in the default container
lookupStats - Look up index statistics on the default container
openContainer - Opens a container, and uses it as the default container
preload - Pre-loads (opens) a container
prepare - Prepare the given query expression as the default pre-parsed query
print - Prints most recent results, optionally to a file
putDocument - Insert a document into the default container
query - Execute the given query expression, or the default pre-parsed query
queryPlan - Prints the query plan for the specified query expression
quit - Exit the program
reindexContainer - Reindex a container, optionally changing index type
removeAlias - Remove an alias from the default container
removeContainer - Removes a container
removeDocument - Remove a document from the default container
removeNodes - Remove content from documents specified by the query expression
renameNodes - Rename nodes specified by the query expression
run - Runs the given file as a script
setBaseUri - Set/get the base uri in the default context
setIgnore - Tell the shell to ignore script errors
setLazy - Sets lazy evaluation on or off in the default context
setMetaData - Set a metadata item on the named document
setNamespace - Create a prefix->namespace binding in the default context
setProjection - Enables or disables the use of the document projection optimization
setQueryTimeout - Set a query timeout in seconds in the default context
setReturnType - Sets the return type on the default context
setTypedVariable - Set a variable to the specified type in the default context
setVariable - Set a variable in the default context
setVerbose - Set the verbosity of this shell
sync - Sync current container to disk
time - Wrap a command in a wall-clock timer
transaction - Create a transaction for all subsequent operations to use
updateNodes - Update node content based on query expression and new content
upgradeContainer - Upgrade a container to the current container format

Printing the names of all the SNP

dbxml> query 'collection("dbsnp.dbxml")/snp/name/string()'
6 objects returned for eager expression 'collection("dbsnp.dbxml")/snp/name/string()'


dbxml> print
rs25
rs26
rs27
rs300
rs600
rs800

Finding the observed bases for the SNPs having het>0.3

dbxml> query 'collection("dbsnp.dbxml")/snp[number(het) > 0.3 ]/observed/string()'
2 objects returned for eager expression 'collection("dbsnp.dbxml")/snp[number(het) > 0.3 ]/observed/string()'


dbxml> print
A/G
C/G

Printing a HTML table of all the SNPs on chrom7

dbxml> query '<html><body><table><tr><th>Name</th><th>Chrom</th><th>Position</th></tr>
{ for $location in collection("dbsnp.dbxml")/snp/mapping/location[@chrom="7" and @label="reference"]
return
<tr><th>{$location/../../name/string()}</th><td>7</td><td>{$location/@position/string()}</td></tr>
}
</table></body></html>'
1 objects returned for eager expression '<html><body><table><tr><th>Name</th><th>Chrom</th><th>Position</th></tr> { for $location in collection("dbsnp.dbxml")/snp/mapping/location[@chrom="7" and @label="reference"]
return
<tr><th>{$location/../../name/string()}</th><td>7</td><td>{$location/@position/string()}</td></tr>}</table></body></html>'


dbxml> print
<html><body><table><tr><th>Name</th><th>Chrom</th><th>Position</th></tr><tr><th>rs25</th><td>7</td><td>11550666</td></tr><tr><th>rs26</th><td>7</td><td>11549995</td></tr><tr><th>rs27</th><td>7</td><td>11549749</td></tr></table></body></html>

Result viewed in a browser:
NameChromPosition
rs25711550666
rs26711549995
rs27711549749


That's it
Pierre