'R' = dna.translate("AGG") . A custom C function for R, My notebook.

In the following post, I will show how I've implemented a custom C function for R. This C function will translate a DNA to a protein. I'm very new to 'R' so feel free to make any comment about the code.

C code

The data in 'R' are stored in an opaque structure named 'SEXP'. A custom C function receives some SEXP arguments and returns another SEXP. So, the declaration of my function translate_dna is:

SEXP translate_dna(SEXP sequences)
{
(...)
}

We check if the argument 'sequences' has a type='character'.

if(!isString(sequences))
error("argument is not a string[]");

The input consists of one or more DNA, so our return value will be an SEXP array of type 'character' (STRSXP) with a size equals to the number of DNAs. This array is allocated:

SEXP array_of_peptides = allocVector(STRSXP,length(sequences));

We loop over each DNA sequence

for(in i=0;i<length(sequences);++i) ...

and we allocate some memory for the peptide:

const char *dna=CHAR(STRING_ELT(sequences, i));
int dna_length=strlen(dna);
char* peptide =malloc(dna_length/3+1);

The peptide is then filled with its amino acids:

for(j=0;j+2 < dna_length;j+=3)
{
peptide[n++]=_translate(dna[j],dna[j+1],dna[j+2]);
}
peptide[n]=0;

And we put this new string/peptide in the returned value

SET_STRING_ELT(array_of_peptides,i,Rf_mkChar(peptide));

Here, I am not sure how 'R' handles its memory. Should I care about the way R runs its garbage manager ? how should I use the macros PROTECT and UNPROTECT ?

Compilation

The code 'translate.c' is compiled as a dynamic library ' libtranslate.so':

gcc -fPIC -I -g -c -Wall -I ${R_HOME}/include translate.c
gcc -shared -Wl,-soname,libtranslate.so.1 -o libtranslate.so translate.o

R code

on the R side , the previous dynamic library is loaded:

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

A function dna.translate is declared: it forces the array to be an array of string and it invokes the C function 'translate_dna'

storage.mode(dna) <- "character"
.Call("translate_dna", dna)

We can now invoke the R function dna.translate:

peptides <- dna.translate(
c( "ATGGAGAGGCAGAAACGGAAGGCGGACATCGAGAAAGGGCTGCAGTTCATTCAGTCGACACTAC",
NULL,
"CCCAAAAGCAAGAAGAATATGAGGCCTTTCTGCT",
"CAAACTGGTGCAGAATCTGTTTGCTGAGGGCAATGA",
NULL,
"GGCAGATCAGGGAACTTCTAATGGATTGGGGTCC",
"GGATAACTGCACCTTCGCCTACCATCAGGAGGA",
"GGGTCCCAGGCAGCGCTGCCTGGGGGCTGGGGAG",
"TTGCGACAGGCTCCAGAAGGGCAAAGCCTGCCCAGAT",
"GCACCCCCCTCCTCTCCACCCTACCTTCCATCAACCA",
"AGG"
))

print(peptides)

Result:

[1] "MERQKRKADIEKGLQFIQSTL" "PKSKKNMRPFC" "QTGAESVC*GQ*"
[4] "GRSGNF*WIGV" "G*LHLRLPSGG" "GSQAALPGGWG"
[7] "LRQAPEGQSLPR" "APPSSPPYLPST" "R"

Full source code

:
C code:

#include <stdio.h>
#include <ctype.h>
#include <R.h>
#include <Rinternals.h>
/* the genetic code */
static const char* STANDARD_GENETIC_CODE="FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG";

static char _translate(char a,char b,char c);


/** translates DNA to protein
* @arg sequences one or more DNA sequence
* @return an array of peptides
*/
SEXP translate_dna(SEXP sequences)
{
int i;
SEXP array_of_peptides;
//check input has type= characters
if(!isString(sequences))
error("argument is not a string[]");
//prepare a new list for the translated sequence
array_of_peptides = allocVector(STRSXP,length(sequences));

//loop over the input
for(i=0;i< length(sequences);++i)
{
int j=0,n=0;
//transform the sequence into a ptr* of char
const char *dna=CHAR(STRING_ELT(sequences, i));
//get the length of this sequence
int dna_length=strlen(dna);
//alloc the protein sequence
char* peptide =malloc(dna_length/3+1);
if(peptide==NULL) error("out of memory");
//loop over the codons
for(j=0;j+2 < dna_length;j+=3)
{
//set the amino acid at 'n'
peptide[n++]=_translate(dna[j],dna[j+1],dna[j+2]);
}
//add EOS
peptide[n]=0;
//put a copy of this peptide in the vector
SET_STRING_ELT(array_of_peptides,i,Rf_mkChar(peptide));
//free our copy
free(peptide);
}
//return the array of petptide
return array_of_peptides;
}


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

static char _translate(char a,char b,char c)
{
int base1= base2index(a);
int base2= base2index(b);
int base3= base2index(c);
if(base1==-1 || base2==-1 || base3==-1)
{
return '?';
}
else
{
return STANDARD_GENETIC_CODE[base1*16+base2*4+base3];
}

}

Makefile:

run:translate.so
${R_HOME}/bin/R --no-save < translate.R

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

R code:

#load the dynamic library translate
dyn.load(paste("libtranslate", .Platform$dynlib.ext, sep=""))
#declare the function dna.translate
dna.translate <- function(dna)
{
#force dna to be type=character
storage.mode(dna) <- "character"
#call the C function
.Call("translate_dna", dna)
}

peptides <- dna.translate(
c( "ATGGAGAGGCAGAAACGGAAGGCGGACATCGAGAAAGGGCTGCAGTTCATTCAGTCGACACTAC",
NULL,
"CCCAAAAGCAAGAAGAATATGAGGCCTTTCTGCT",
"CAAACTGGTGCAGAATCTGTTTGCTGAGGGCAATGA",
NULL,
"GGCAGATCAGGGAACTTCTAATGGATTGGGGTCC",
"GGATAACTGCACCTTCGCCTACCATCAGGAGGA",
"GGGTCCCAGGCAGCGCTGCCTGGGGGCTGGGGAG",
"TTGCGACAGGCTCCAGAAGGGCAAAGCCTGCCCAGAT",
"GCACCCCCCTCCTCTCCACCCTACCTTCCATCAACCA",
"AGG"
))

print(peptides)

That's it

Pierre

Playing with Biomart. My notebook.

Thanks to Neil, I've discovered that Biomart was not just a web interace, it also contains a powerful API: Citing http://www.biomart.org/martservice.html:
To submit a query using our webservices generate an XML document conforming to our Query XML syntax. This can be achieved simply by building up your query using MartView and hitting the XML button . This XML should be posted to http://(...)/martservice attached to a single parameter of query.

The base URL for Biomart is http://www.biomart.org/biomart.

Retrieving the Registry information for this Biomart installation: http://www.biomart.org/biomart/martservice?type=registry

<MartRegistry>
<MartURLLocation database="ensembl_mart_57" default="1" displayName="ENSEMBL GENES 57 (SANGER UK)" host="www.biomart.org" includeDatasets="" martUser="" name="ensembl" path="/biomart/martservice" port="80" serverVirtualSchema="default" visible="1" />
<MartURLLocation database="hapmart_27_36" default="1" displayName="HAPMAP 27 (NCBI US)" host="hapmap.ncbi.nlm.nih.gov" includeDatasets="" name="HapMap_rel27" path="/biomart/martservice" port="80" redirect="1" serverVirtualSchema="rel27_NCBI_Build36" visible="1" />
<MartURLLocation database="vega_mart_57" default="0" displayName="VEGA 37 (SANGER UK)" host="www.biomart.org" includeDatasets="" martUser="" name="vega" path="/biomart/martservice" port="80" serverVirtualSchema="default" visible="1" />
<MartURLLocation database="genomic_features_mart_57" default="0" displayName="ENSEMBL GENOMIC FEATURES 57 (SANGER UK)" host="www.biomart.org" includeDatasets="" martUser="" name="genomic_features" path="/biomart/martservice" port="80" serverVirtualSchema="default" visible="0" />
(...)
</MartRegistry>

Retrieving the datasets available for the registry HapMap_rel27: http://www.biomart.org/biomart/martservice?type=datasets&mart=HapMap_rel27 :

TableSet hm27_variation_yri HapMap Population YRI - Yoruba in Ibadan, Nigeria (West Africa) 1 20050000 default 2009-10-28 11:19:52
TableSet hm27_variation_tsi HapMap Population TSI - Toscans in Italy 1 200 50000 default 2009-10-28 11:13:07
TableSet hm27_gene hm27_gene 0 200 50000 default 2009-09-17 12:04:20
TableSet hm27_variation_chd HapMap Population CHD - Chinese in Metropolitan Denver, Colorado 1 200 50000 default 2009-10-28 11:07:01
TableSet hm27_variation All Populations 1 200 50000 default 2009-10-28 11:22:21
TableSet hm27_variation_chb HapMap Population CHB - Han Chinese in Beijing, China 1 200 50000 default 2009-10-28 11:05:57
TableSet hm27_variation_asw HapMap Population ASW - African ancestry in Southwest USA 1 20050000 default 2009-10-28 08:53:43
TableSet hm27_variation_ceu HapMap Population CEU - Utah residents with Northern and Western European ancestry from the CEPH collection 1 200 50000 default 2009-10-28 11:04:39
TableSet hm27_variation_mkk HapMap Population MKK - Maasai in Kinyawa, Kenya 1 200 50000 default 2009-10-28 11:12:12
(...)

Where can I find a description of those columns ?


Retrieving the attributes for the dataset hm27_variation_ceu: http://www.biomart.org/biomart/martservice?type=attributes&dataset=hm27_variation_ceu :

chrom chromosome naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main chrom
start position naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main start
strand strand naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main strand
marker1 marker id naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main marker1
alleles alleles naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main alleles
refallele reference allele naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main refallele
center genotyping center naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__maincenter
platform genotyping platform naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main platform
ceu_id CEU genotype naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__dm genotype
ref_allele reference allele naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main ref_allele
ref_allele_freq reference allele frequency naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main ref_allele_freq
ref_allele_count reference allele count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main ref_allele_count
other_allele other allele naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__mainother_allele
other_allele_freq other allele frequency naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main other_allele_freq
other_allele_count other allele count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main other_allele_count
total_allele_count total allele count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main total_allele_count
refhom_gt reference homozygote genotype naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main refhom_gt
het_gt heterozygote genotype naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__mainhet_gt
otherhom_gt other homozygote genotype naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main otherhom_gt
refhom_gtfreq reference homozygote genotype frequency naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main refhom_gtfreq
het_gtfreq heterozygote genotype frequency naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main het_gtfreq
otherhom_gtfreq other homozygote genotype frequency naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main otherhom_gtfreq
refhom_gtcount reference homozygote genotype count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main refhom_gtcount
het_gtcount heterozygote genotype count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main het_gtcount
otherhom_gtcount other homozygote genotype count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main otherhom_gtcount
total_gtcount total genotype count naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__variation__main total_gtcount
prot_lsid protocol LSID naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__assay__dm prot_lsid
assay_id assay LSID naive_attributes html,txt,csv,tsv,xls hm27_variation_ceu__assay__dm assay_lsid

(Where can I find a description of those columns ?). Each attribute is a specific column in the final result. Example of attributes available for the dataset hm27_variation_ceu:

• chrom
• start
• marker1
• alleles
• (...)

Retrieving the filters for the dataset hm27_variation_ceu: http://www.biomart.org/biomart/martservice?type=filters&dataset=hm27_variation_ceu :

allele_freq Minor Allele Frequency [>=] [0.0,0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5] naive_filters list >= hm27_variation_ceu__variation__main avg_dprime
center Genotyping Center [affymetrix,bcm,illumina,imsut-riken,mcgill-gqic,perlegen,ucsf-wu,broad,sanger] naive_filters list = hm27_variation_ceu__variation__main center
platform Genotyping Platform [AFFY_100K,AFFY_500K,BeadArray,FP-TDI,Infinium_H1,Infinium_H300,Invader,MIP,Perlegen,AFFY_6.0,Illumina_1M] naive_filters list = hm27_variation_ceu__variation__main platform
monomorphic Monomorphic SNPs [,0] naive_filters boolean_num =only,excluded hm27_variation_ceu__variation__main is_reference_bool
marker_name List of HapMap SNPs [] naive_filters = hm27_variation_ceu__variation__mainmarker1
utr3 3' UTR [] naive_filters boolean_list only,excluded hm27_variation_ceu__variation__main is_utr_bool
utr5 5' UTR [] naive_filters boolean_list only,excluded hm27_variation_ceu__variation__main is_5utr_bool
coding_synon Exons - synonymous coding SNPs [] naive_filters boolean_list only,excluded hm27_variation_ceu__variation__main is_coding_synon_bool
coding_nonsynon Exons - non-synonymous coding SNPs [] naive_filters boolean_list only,excluded hm27_variation_ceu__variation__main is_coding_nonsynon_bool
intronic introns [] naive_filters boolean_list only,excluded hm27_variation_ceu__variation__mainis_intronic_bool
chrom chromosome [chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr20,chr21,chr22,chrX,chrY] naive_filters list = hm27_variation_ceu__variation__main chrom
start From position [] naive_filters text >= hm27_variation_ceu__variation__main start
stop To position [] naive_filters text <= hm27_variation_ceu__variation__main stop
gene_list List of Genes [] naive_filters =,in pointer dataset gene_primary_symbol
encode_region ENCODE region [ENm010: 7:26924045..27424045,ENm013: 7:89621624..90736048,ENm014: 7:125865891..127029088,ENr112: 2:51512208..52012208,ENr113: 4:118466103..118966103,ENr123: 12:38626476:39126476,ENr131: 2:234156563..234656627,ENr213: 18:23719231..24219231,ENr232: 9:130725122..131225122,ENr321: 8:118882220..119382220] naive_filters list = pointer dataset glook_encode_region

(Where can I find a description of those columns ?) Here we can find some filters for hm27_variation_ceu (e.g.: allele_freq)

Retrieving the configurations for the dataset hm27_variation_ceu: http://www.biomart.org/biomart/martservice?type=configuration&dataset=hm27_variation_ceu :

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE DatasetConfig>
<DatasetConfig dataset="hm27_variation_ceu" datasetID="2" displayName="HapMap Population CEU - Utah residents with Northern and Western European ancestry from the CEPH collection" interfaces="default" internalName="default" martUsers="default" modified="2009-10-28 11:04:39" softwareVersion="0.6" template="hm27_variation_ceu" type="TableSet" visible="1">
<MainTable>hm27_variation_ceu__variation__main</MainTable>
<Key>marker_id_key</Key>
<Importable filters="chrom,start,stop" internalName="encode_filter" linkName="encode_filter" name="encode_filter" type="link"/>
<Importable displayName="gene_filter_35" filters="chrom,start,stop" internalName="gene_filter_35" linkName="gene_filter_35" name="gene_filter_35" type="link"/>
<FilterPage displayName="FILTERS" internalName="naive_filters">
<FilterGroup displayName="FILTERS" internalName="filters">
<FilterCollection displayName="MINOR ALLELE FREQUENCY [&gt;=]" internalName="pop_filter">
<FilterDescription displayName="Minor Allele Frequency [&gt;=]" displayType="list" field="avg_dprime" internalName="allele_freq" key="marker_id_key" legal_qualifiers="&gt;=" qualifier="&gt;=" style="menu" tableConstraint="main" type="list">
<Option displayName="0.0" internalName="0.0" isSelectable="true" value="0.0"/>
<Option displayName="0.01" internalName="0.01" isSelectable="true" value="0.01"/>
<Option displayName="0.05" internalName="0.05" isSelectable="true" value="0.05"/>
<Option displayName="0.1" internalName="0.1" isSelectable="true" value="0.1"/>
<Option displayName="0.15" internalName="0.15" isSelectable="true" value="0.15"/>
<Option displayName="0.2" internalName="0.2" isSelectable="true" value="0.2"/>
<Option displayName="0.25" internalName="0.25" isSelectable="true" value="0.25"/>
<Option displayName="0.3" internalName="0.3" isSelectable="true" value="0.3"/>
<Option displayName="0.35" internalName="0.35" isSelectable="true" value="0.35"/>
(...)

As far as I understand , this file is used to 'group' the filters.

At the end, here is Query sent to Biomart:

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE Query>
<Query virtualSchemaName="default" formatter="TSV" header="0" uniqueRows="0" count="" datasetConfigVersion="0.6">
<Dataset name="hm27_variation_ceu" interface="default">
<Filter name="stop" value="1000000"/>
<Filter name="allele_freq" value="0.3"/>
<Filter name="chrom" value="chr1"/>
<Filter name="start" value="0"/>
<Attribute name="chrom"/>
<Attribute name="start"/>
<Attribute name="strand"/>
<Attribute name="marker1"/>
<Attribute name="alleles"/>
<Attribute name="ceu_id"/>
</Dataset>
</Query>

and we can send it with curl:

curl --data 'query=<?xml version="1.0" encoding="UTF-8"?><!DOCTYPE Query><Query virtualSchemaName="default" formatter="TSV" header="0" uniqueRows="0" count="" datasetConfigVersion="0.6" ><Dataset name="hm27_variation_ceu" interface="default" ><Filter name="stop" value="1000000"/><Filter name="allele_freq" value="0.3"/><Filter name="chrom" value="chr1"/><Filter name="start" value="0"/><Attribute name="chrom" /><Attribute name="start" /><Attribute name="strand" /><Attribute name="marker1" /><Attribute name="alleles" /><Attribute name="ceu_id" /></Dataset></Query>' http://www.biomart.org/biomart/martservice

Response from the Biomart server:

chr1 843817 + rs1806509 A/C AC AA AC AC AC CC AC CC AC CC AC NN AC CC AA NN AC AC NN CC AC CC AC AC AA AC AC AC CC AA AC AC AC AC AA CC AC AC AA AA AA AA NN AA AA AA AC AA NN AA NN AC AC CC AA AA AC AC AC AC AA AA AA AC AA AA CC CC AC AA AC AC AC AA CC AA AC CC AA NN CC CC AA AC AA AC AA AA AC AC AC AC AA AA NN AA AA CC AC AC AA AA AA AC AC AC AA AC AA AC AC AA AA AC AA CC AC AC AC AC AA AA AA CC AC CC AC NN AC AC AA AA CC AA AC AC AC AA AC AA AC AC AC AA AC CC AA AA CC AA AA AC AC AC AC AA CC AA AC CC CC AC AC AA CC AA AC AC AC AA AA CC AA AA
chr1 881808 + rs13303106 A/G AG GG AG GG GG AG AG AG GG GG GG NN AA AG GG NN GG GG NN GG AG AG GG AG GG GG AG GG AA GG AG GG AG AA GG AG GG AG GG GG GG GG NN GG GG GG AA AG NN GG NN AG AG AA GG AG AG AA AG AG GG GG GG AG GG GG AG AA GG GG AG GG AG GG AG GG AG AG GG NN AA GG GG AG GG AG AG AG AG AG GG AG GG GG NN GG GG AG AG AG AG GG AG AG GG GG GG AG GG AG AG GG AG AG AG AA GG AG AG AG GG GG GG AG AG AA AG NN AG AG GG GG AA GG GG GG AA GG GG AG AG AG AG GG AG AA GG GG AA GG GG AG AG AG GG GG AA GG GG AA AA AG AG AG AA GG AA AG AG GG AG AA GG GG
chr1 906697 + rs6694632 A/G AG AA AA AG AA AG AG AG AA AA AG NN GG AA GG NN AA AA NN AA AG AG AA AG AA AA AG AA GG GG GG AG AA GG AA AA AA AG AG AG AG AA NN AA AA AA AG AG NN AA NN AG AG AG AG AG AG GG AG AG AA AA AA AG AG AG AA GG AA AA AA GG GG AA AG GG AA AG AA NN AG AA AG AG AA AG GG GG AA GG AG AG AG AG NN AA AG AG AG AG AG AA AA GG AA AA AA AA AA AG AG AA AG AG AG AA AG AA AG AG AA AA GG AG AG GG AG NN AG AG GG AG AA AG AG AA AG AA AA AA AG AA AG AA AA GG AA AA GG AA AA AG AG AG AA AA GG AA AA GG GG AG AG AG AG AG GG AA AG AG AG AG AG AG
chr1 908247 + rs13303118 G/T GT TT TT GT TT GT GT GT TT TT GT NN GG TT GG NN TT TT NN TT GT GT TT GT TT TT GT TT GG GG GG GT TT GG TT TT TT GT GT GT GT TT NN TT TT TT GT GT NN TT NN NN GT GT GT GT GT GG GT GT TT TT TT GT GT GT TT GG TT TT TT GG GG TT GT GG TT GT TT NN GT TT GT GT TT GT GG GG TT GG GT GT GT GT NN TT GT GT GT GT GT TT TT GG TT TT TT TT TT GT GT TT GT GT GT TT GT TT GT GT TT TT GG GT GT GG GT NN GT GT GG GT TT GT GT TT GT TT TT TT GT TT GT TT TT GG TT TT GG TT TT GT GT GT TT TT GG TT TT GG GG GT GT GT GT GT GG TT GT GT GT GT GT GT
chr1 908436 + rs2341354 A/G AG GG GG AG GG AG AG AG GG GG AG NN AA GG AA NN GG GG NN GG AG AG GG AG GG GG AG GG AA AA AA AG GG AA GG GG GG AG AG AG AG GG NN GG GG GG AG AG NN GG NN AG AG AG AG AG AG AA AG AG GG GG GG AG AG AG AG AA GG GG GG AA AA GG AG AA GG AG GG NN AG GG AG AG GG AG AA AA GG AA AG AG AG AG NN GG AG AG AG AG AG GG GG AA GG GG GG GG GG AG AG GG AG AG AG GG AG GG AG AG GG GG AA AG AG AA AG NN AG AG AA AG GG AG AG GG AG GG GG GG AG GG AG GG GG AA GG GG AA GG GG AG AG AG GG GG AA GG GG AA AA AG AG AG AG AG AA GG AG AG AG AG AG AG
chr1 934427 + rs3128117 C/T TT TT CT CC TT CT CC CT TT TT CT NN CT CT CC NN TT TT NN TT TT CT TT CC TT TT CT TT CC CC TT TT TT CT TT TT CT CC CT CT CT TT NN TT CT TT CT CT NN TT NN CT CC TT CT CT CT CC CT CT TT CT TT CT CT CT TT CT TT TT CT CC CT TT CT CC CT CT TT NN CT TT CT CT CT CT CC CC TT CC CT TT TT CC NN TT TT CT CC CT TT TT TT CT CT TT TT TT TT TT TT TT CT CT CT CT CT CT CC CT CT TT CC CT CT CT CT NN CC CC CC CT CT CT CC TT CT TT TT TT CT TT CC CT TT CC CT TT CC TT TT CT CC CT TT TT CT TT TT CT CC TT CT TT CT CT CT CT TT CT TT TT CT CT
chr1 940106 + rs1891906 A/C AA AA AC CC AA AC AC AC AA AA AC NN AC AC CC NN AA AA NN AA AA AC AA CC AA AA AC AA CC CC AA AA AA AC AA AA AC CC AC AC AC AA NN AA AC AA AC AC NN AA NN AC CC AA AC AC AC CC AC AC AA AC AA AC AC AC AA AC AA AA AC CC AC AA AC AC AC AC AA NN AC AA AC AC AC AC CC CC AA NN AC AA AA CC NN AA AA AC CC AC AA AA AA AC AC AA AA AA AA AA AA AA AC AC AC AC AC AC CC AC AC AA CC AC AC AC AC NN CC CC AC AC AC AC CC AA AC AA AA AA AC AA CC AC AA CC AC AA CC AA AA AC CC AC AA AA AC AA AA AC CC AA AC AA AC AC AC AC AA AC AA AA AC AC
chr1 949705 + rs2710888 C/T CT CC CT CT CC CT CT CT CC CC CT NN CT CT TT NN CC CC NN CC CC CC NN CT CC CC CT CC TT CT CC CC CT CC CC CC CT TT CT CT CT CC NN CC CC CT TT CT NN CC NN CT TT CC CT CC CT TT CT CT CT CC CC CT CC CT CC CT CC CC CT CT CT CT CT CC CT CT CC NN CT CC CT CT CT CT CT TT CT TT CC CC CC TT NN CC CC CT TT CC CC CC CT CT CT CC CC CT CC CC CT CC CT CT CT TT CC CT TT CT CT CC CT CT CT CT CC NN CC TT CT CT CT CT CT CC CT CC CC CC CT CT TT CT TT TT CT CC TT CC CC CT TT CT CC CC CT CC CC CT TT CC CT CC CT CT CT CT CT CC CT CT CC CT
chr1 952073 + rs3128126 A/G NN AA NN NN AA AG AG NN NN AA NN AG AG AG NN AG NN NN AG NN AG AA AG NN NN AA NN AA NN AG AA AA NN NN AA AG NN NN NN AG AG NN AA NN NN AA AG AG AA AA AG AG GG NN NN AA AG GG AA AG AA NN AG AG NN NN NN NN NN NN NN NN NN NN AG AG AG AG AA AG GG AA AG AG NN AG AG GG AG GG AA AA AA GG GG AA AG AG GG NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN NN GG NN AG GG NN NN AG NN NN AG AG AA NN AA AA AA AG NN NN NN NN NN NN AA GG AA AA AG GG NN NN NN NN NN NN NN NN NN GG AA AG GG AG AG NN AA NN NN AA AG
chr1 952469 + rs4970393 A/G AG AA AG GG AA AG AG AG AA AA AG NN AG AG AG NN AA AA NN AA AG AA AG AG AA AA AG AA GG GG AA AA AA GG AA AA AG GG AG AG AG AA NN AA AG AA AG AG NN AA NN AG GG AG AG AG AG GG AA AG AA AG AA AG AG AG AA GG AA AA AG AG AG AG AG AG AG AG AA NN AG AA AG AG AG AG GG NN AG GG AG AA AA GG NN AA AA AG GG AG AG AA AA AG AG AA AA AG AA AA AG AA AG AG AG AG AG AG GG AG AG AG NN AG AG AG AA NN AG NN GG AG AG AG GG AA AG AA AA AA AG AG GG AG GG GG AG AA GG AA AA AG GG AG AA AA AG AA AA GG AG AG AG AG AG AG NN AG AA AG AA AA AG AG

Note: I've also tried to use the SOAP/WSDL Biomart API, but as it was discussed on Biostar, the API does not work.


That's it

Pierre

Learning DesignPatterns: The Factory Method Pattern

"The Factory method pattern defines a separate method for creating the objects, which subclasses can then override to specify the derived type of product that will be created"

Example

/** Abstract class for melting temperature calculation */
abstract class TmCalculator
{
public abstract double calcTm(CharSequence s);
}

/** the factory creates two kinds of TmCalculator
1) 2AT4GC
2) Nearest-Neighbor
*/
public class TmCalculatorFactory
{
public TmCalculator create2AT4GC()
{
return new TmCalculator()
{
@Override
public double calcTm(CharSequence s)
{
(...implements ... )
}
};
}

public TmCalculator createNearestNeighbor()
{
return new TmCalculator()
{
@Override
public double calcTm(CharSequence s)
{
(...implements ... )
}
};
}
}

That's it

Pierre

Biostar.com: A StackExchange engine for Bioinformatics

Just In Case You Don't Know:

Biostar is a site running the StackExchange engine (the same engine used by http://stackoverflow.com). The site's focus is bioinformatics, computational genomics and biological data analysis(...)No question is too trivial or too "newbie". Please look around to see if your question has already been asked (and maybe even answered!) before you ask. But if you end up asking a question that has been asked before, that is OK and deliberately allowed. Other users will hopefully edit in links to related or similar questions to help future visitors find their way..

Learning DesignPatterns: The Builder Pattern

largely inspired by Wikipedia:"Builder Pattern is a software design pattern. The intention is to abstract steps of construction of objects so that different implementations of these steps can construct different representations of objects."

Example

/**
* Taxon, the class to be constructed
*/
class Taxon
{
private int id=-1;
private int parent_id=-1;
private String name=null;

public void setId(int id) { this.id=id; }
public void setParentId(int parent_id) { this.parent_id=parent_id; }
public void setName(String name) { this.name=name; }
public int getId() { return this.id;}
public int getParentId() { return this.parent_id;}
public String getName() { return this.name;}
}

/** Abstract interface for creating Taxons */
abstract class TaxonBuilder
{
protected Taxon taxon=new Taxon();

public void createNewTaxon()
{
this.taxon = new Taxon();
}

public Taxon getTaxon()
{
return this.taxon;
}
public abstract void buildId();
public abstract void buildParentId();
public abstract void buildName();
}
/** concrete TaxonBuilder builder for HomoSapiens */
class HomoSapiensBuilder extends TaxonBuilder
{
public void buildId() { super.taxon.setId(9606);}
public void buildParentId() { super.taxon.setParentId(9605);}
public void buildName() { super.taxon.setName("Homo Sapiens");}
}

/** concrete TaxonBuilder builder for Platypus */
class PlatypusBuilder extends TaxonBuilder
{
public void buildId() { super.taxon.setId(9258);}
public void buildParentId() { super.taxon.setParentId(9257);}
public void buildName() { super.taxon.setName("Platypus");}
}

/** constructs a Taxon object by calling its TaxonBuilder */
class TaxonAssembler
{
private TaxonBuilder builder;

public void setTaxonBuilder(TaxonBuilder builder)
{
this.builder=builder;
}

public void assemble()
{
this.builder.createNewTaxon();
this.builder.buildId();
this.builder.buildParentId();
this.builder.buildName();
}
}

public class Test
{
public static void main(String[] args) {
{
TaxonBuilder builder=new PlatypusBuilder();
TaxonAssembler assembler=new TaxonAssembler();
assembler.setBuilder(builder);
assembler.assemble();
Taxon taxon=builder.getTaxon();
}
}

That's it.

Pierre

Learning DesignPatterns: AbstractFactoryPattern

In the very next posts, I'll post my notes about the "Design Patterns". A Design Pattern is a "general reusable solution to a commonly occurring problem in software design". Here, I'll follow the patterns described at http://c2.com/cgi/wiki?DesignPatternsBook and will try to apply them to Bioinformatics. Today I'm starting with the AbstractFactoryPattern.

via Wikipedia:AbstractFactoryPattern provides a way to encapsulate a group of individual factories that have a common theme. In normal usage, the client software creates a concrete implementation of the abstract factory and then uses the generic interfaces to create the concrete objects that are part of the theme.
Use of this pattern makes it possible to interchange concrete classes without changing the code that uses them, even at runtime.

Example

/** the result of a pairwise alignment */
interface Alignment
{
public void print(PrintWriter out);
public double getScore();
}

/** interface for a pairwise alignment */
interface PairwiseAlignment
{
public Alignment align(CharSequence seq1,CharSequence seq2);
}

/** implementation of a global PairwiseAlignment */
class NeedlemanWunsch
implements PairwiseAlignment
{
public Alignment align(CharSequence seq1,CharSequence seq2)
{
(...)
}
}

/** implementation of a local PairwiseAlignment */
class SmithWaterman
implements PairwiseAlignment
{
public Alignment align(CharSequence seq1,CharSequence seq2)
{
(...)
}
}

/** abstract PairwiseAlignmentFactory
* contains an abstract method 'createPairwiseAligner'
* the static function 'newInstance' returns a new 'PairwiseAlignmentFactoryImpl'
*/
abstract class PairwiseAlignmentFactory
{
protected boolean global=true;

protected PairwiseAlignmentFactory()
{
}

public void setGlobal(boolean global)
{
this.global=global;
}

public boolean isGlobal()
{
return this.global;
}

/** abstract. to be implemented */
public abstract PairwiseAlignment createPairwiseAligner();

/** creates a new concrete instance of PairwiseAlignmentFactory */
static public PairwiseAlignmentFactory newInstance()
{
//described later...
return new PairwiseAlignmentFactoryImpl();
}
}

/** a concrete PairwiseAlignmentFactory
* implements createPairwiseAligner()
* instance created by PairwiseAlignment.newInstance()
*/
class PairwiseAlignmentFactoryImpl
extends
{
@Override
public PairwiseAlignment createPairwiseAligner()
{
return isGlobal() ?
new NeedlemanWunsch():
: new SmithWaterman();
}
}

That's it

Pierre

A Mysql user defined function (UDF) for Gene Ontology (GO)

In this post I'll create a Mysql Defined function (UDF) answering if a defined GO term is a descendant of another term. This post is not much different from two previous posts I wrote here:

• MYSQL UDF, trees of data, hierarchy
• Mysql user defined function (UDF) for Bioinformatics.

Here I built a binary file containing an array of pairs(term,parent_term) of GO identifiers from the XML/RDF file o_daily-termdb.rdf-xml.gz. This file is then used by the mysql UDF to find all the ancestors of a given term. The source code is available here:

http://code.google.com/p/lindenb/source/browse/trunk/proj/go_udf/

Building the database

The libxml API is used to load go_daily-termdb.rdf-xml.gz. For each go:term, the parent terms are searched and each pair(term,parent-term) is saved into the array.

static void scanterm(xmlNode *term)
{
xmlAttrPtr rsrc;
xmlNode *n=NULL;
xmlChar * acn=NULL;
for(n=term->children; n!=NULL; n=n->next)
{
if (n->type != XML_ELEMENT_NODE) continue;
if(strcmp(n->name,"accession")==0 &&
n->ns!=NULL && n->ns->href!=NULL &&
strcmp(n->ns->href,GO_NS)==0)
{
acn= xmlNodeGetContent(n);
break;
}
}
if(acn==NULL) return;

for(n=term->children; n!=NULL; n=n->next)
{
if (n->type != XML_ELEMENT_NODE) continue;
if(strcmp(n->name,"is_a")==0 &&
n->ns!=NULL && n->ns->href!=NULL &&
strcmp(n->ns->href,GO_NS)==0
)
{
xmlChar* is_a=xmlGetNsProp(n,"resource",RDF_NS);
if(is_a!=NULL)
{
char* p=strstr(is_a,"#GO:");
if(p!=NULL)
{
++p;
termdb.terms=(TermPtr)realloc(termdb.terms,sizeof(Term)*(termdb.n_terms+1));
if(termdb.terms==NULL)
{
fprintf(stderr,"out of memory\n");
exit(EXIT_FAILURE);
}
strncpy(termdb.terms[termdb.n_terms].parent,p,MAX_TERM_LENGTH);
strncpy(termdb.terms[termdb.n_terms].child,acn,MAX_TERM_LENGTH);
++termdb.n_terms;
}
xmlFree(is_a);
}
}
}
xmlFree(acn);
}

When the RDF file has been read, the content of the array is sorted by term and saved to a file

qsort(termdb.terms,termdb.n_terms,sizeof(Term),cmp_term);
out= fopen(argv[1],"w");
if(out==NULL)
{
fprintf(stderr,"Cannot open %s\n",argv[1]);
return -1;
}
fwrite(&(termdb.n_terms),sizeof(int),1,out);
fwrite(termdb.terms,sizeof(Term),termdb.n_terms,out);
fflush(out);
fclose(out);

Initializing the mysql UDF

When the mysql UDF is inited the binary file is loaded into memory (error handling ignored)

my_bool go_isa_init(
UDF_INIT *initid,
UDF_ARGS *args,
char *message
)
{
TermDBPtr termdb;
FILE* in=NULL;

initid->maybe_null=1;
initid->ptr= NULL;

termdb=(TermDBPtr)malloc(sizeof(TermDB));


in=fopen(GO_PATH,"r");
fread(&(termdb->n_terms),sizeof(int),1,in);
termdb->terms=malloc(sizeof(Term)*(termdb->n_terms));
fread(termdb->terms,sizeof(Term),termdb->n_terms,in);
fclose(in);
initid->ptr=(void*)termdb;
return 0;
}

Disposing the mysql UDF

When the UDF is disposed, we just free the memory

/* The deinitialization function */
void go_isa_deinit(UDF_INIT *initid)
{
/* free the memory **/
if(initid->ptr!=NULL)
{
TermDBPtr termdb=(TermDBPtr)initid->ptr;
if(termdb->terms!=NULL) free(termdb->terms);
free(termdb);
initid->ptr=NULL;
}
}

invoking the UDF

The UDF itself will scan the GO directed tree and will get all the ancestors of a given term:

long long go_isa(UDF_INIT *initid, UDF_ARGS *args,
char *is_null, char *error)
{
(...)
index=termdb_findIndexByName(termdb,term);
if(index==-1)
{
return 0;
}
return recursive_search(termdb,index,parent);
}

static int recursive_search(const TermDBPtr db,int index, const char* parent)
{
int rez=0;
int start=index;
int parent_idx=0;

if(start<0 || start>=db->n_terms) return 0;
if(strcmp(db->terms[index].child,parent)==0) return 1;
while(index < db->n_terms)
{
if(strcmp(db->terms[index].child,db->terms[start].child)!=0) break;
if(strcmp(db->terms[index].parent,parent)==0) return 1;
parent_idx= termdb_findIndexByName(db,db->terms[index].parent);
rez= recursive_search(db,parent_idx,parent);
if(rez==1) return 1;
++index;
}
return 0;
}

As the terms have been ordered by their names, a binary_search is used to find the index of a given term:

static int lower_bound(const TermDBPtr termsdb, const char* name)
{
int low = 0;
int len= termsdb->n_terms;

while(len>0)
{
int half=len/2;
int mid=low+half;
if( strncmp(termsdb->terms[mid].child,name,MAX_TERM_LENGTH)<0)
{
low=mid;
++low;
len=len-half-1;
}
else
{
len=half;
}
}
return low;
}


static int termdb_findIndexByName(const TermDBPtr termsdb,const char* name)
{
int i=0;
if(name==NULL || termsdb==NULL || termsdb->terms==NULL || termsdb->n_terms==0) return -1;
i= lower_bound(termsdb,name);
if(i<0 || i >= termsdb->n_terms || strcmp(termsdb->terms[i].child,name)!=0) return -1;
return i;
}

Compiling

On my laptop, the UDF was compiled using

gcc -shared -DGO_PATH='"/tmp/terms.bin"' -I /usr/include -I /usr/local/include -I /usr/include/mysql -o /usr/lib/go_udf.so src/go_udf.c

Creating the function

mysql> create function go_isa RETURNS INTEGER SONAME 'go_udf.so';

Invoking the UDF function

Finf all the GO terms that are a descendant of GO:0016859 (cis-trans isomerase activity):

mysql> select LEFT(T.name,20) as name,T.acc from go_latest.term as T where go_isa(T.acc,"GO:0016859");
+----------------------+------------+
| name | acc |
+----------------------+------------+
| peptidyl-prolyl cis- | GO:0003755 |
| retinal isomerase ac | GO:0004744 |
| maleylacetoacetate i | GO:0016034 |
| cis-trans isomerase | GO:0016859 |
| cis-4-[2-(3-hydroxy) | GO:0018839 |
| trans-geranyl-CoA is | GO:0034872 |
| carotenoid isomerase | GO:0046608 |
| 2-chloro-4-carboxyme | GO:0047466 |
| 4-hydroxyphenylaceta | GO:0047467 |
| farnesol 2-isomerase | GO:0047885 |
| furylfuramide isomer | GO:0047907 |
| linoleate isomerase | GO:0050058 |
| maleate isomerase ac | GO:0050076 |
| maleylpyruvate isome | GO:0050077 |
| retinol isomerase ac | GO:0050251 |
+----------------------+------------+
15 rows in set (1.27 sec)

Disposing the UDF function

drop function go_isa;

That's it
Pierre

eXist: The Open Source Native XML Database : My notebook

In a previous post, I've played with Oracle's BerkeleyDB-XML. Here, I used with eXist-db, an open source database management system built using XML technology. It stores XML data according to the XML data model and features efficient, index-based XQuery processing.

Download & Install

wget http://downloads.sourceforge.net/project/exist/Stable/1.4/eXist-setup-1.4.0-rev10440.jar
java -jar eXist-setup-1.4.0-rev10440.jar
export EXIST_HOME=PATH_TO_EXIST/eXist

And tha'ts it: it was far more easy than installing (compiling...) BerkeleyDB-XML.

Starting the Server

eXist/bin/startup.sh

Using locale: en_US.UTF-8
18 Feb 2010 15:31:32,579 [main] INFO (JettyStart.java [run]:90) - Configuring eXist from EXIST/eXist/conf.xml
18 Feb 2010 15:31:32,580 [main] INFO (JettyStart.java [run]:91) -
18 Feb 2010 15:31:32,580 [main] INFO (JettyStart.java [run]:92) - Running with Java 1.6.0_07 [Sun Microsystems Inc. (Java HotSpot(TM) Server VM) in /usr/local/package/jdk1.6.0_07/jre]
(...)

Inserting the data

First using the web console, I've created a 'collection' named '/db/dbsnp'.
I've then downloaded about 1000 XML documents from dbsnp:

for S in `mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e 'select right(name,length(name)-2) from snp130 limit 1000' -N`
do
wget -O rs${S}.xml "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=snp&id=${S}&retmode=xml"
done

Those document have been inserted in the XML database:

EXIST/eXist/bin/client.sh -u admin --password=mypassword -s -m 'db/dbsnp' -p rs*.xml
(...)
parsed 7704 bytes in 33ms.
storing document rs10.xml (1 of 1) ...done.
parsing 16306 bytes took 36ms.

parsed 16306 bytes in 36ms.

Using XQUERY

The following XQuery search for all the SNPs in the database '/db/dbsnp' having a heterozygosity greater than 0.49. For each SNP it prints its name, its sequence and the position on the reference genome.

xquery version "1.0";

declare namespace s="http://www.ncbi.nlm.nih.gov/SNP/docsum";

<MyListOfSnp>
{
for $x in collection("/db/dbsnp")/s:ExchangeSet/s:Rs
where data($x/s:Het/@value)>0.49
return <SNP>
<name>rs{data($x/@rsId)}</name>
<sequence>{data($x/s:Sequence/s:Seq5)}[{data($x/s:Sequence/s:Observed)}]{data($x/s:Sequence/s:Seq3)}</sequence>

{
for $as in $x/s:Assembly
where $as/@groupLabel="reference" return

for $comp in $as/s:Component return

for $maploc in $comp/s:MapLoc
return
<map>
<chromosome>chr{data($comp/@chromosome)}</chromosome>
<position>{data($maploc/@physMapInt)}</position>
</map>
}
</SNP>
}
</MyListOfSnp>

Executing the query:

EXIST/eXist/bin/client.sh -u admin --password=mypassword -F input.xquery

Result:

<MyListOfSnp>
<SNP>
<name>rs10000300</name>
<sequence>ATCAAATACCCAAGCAAAGATTTACATTCAAATCTGTTTACTGAAGTTCTATTTATAATACAATGCAATGAACATAATAGTATATATTTACACGTAATGTAATAAACACAAATATTCAATGGTATAAAAATGGTCAATAAATCGTGGCATAGCCACAGCTTAGAGTACCTGTTTAATGTTCTCAGCTATTTTAACTTTGCTAAATAATATTTAAAGATATGcggtagtcccccttcatctgaggaggacctgttccaagacccccagtggatgcctgaaacctctgatagtaatgaaccctatatatactgttttttcctatacatatttacatatataatacatacctatgattaagtttaatttataaattaggcacagtaagagattaacaacaacaataataaaatgtaacaattatagcaatactctaataataaagttatgtgagtgtggtctctctctctgtctcaaaatatcatactgtatgcctctatttt[G/T]ggaatacagttgacaacgggtaactgaaaccgagaaaagtgaaactgcagatgggggctgactactgTATATGAAAATTAAACAATCagccaggcatggtggctcacgcctgtaataccagcactttgggaggccgaggcgggaggatcacgaggtcaggagatcgagaccacggtgaaaccccgtctctattaaaaatacaaaaaaaaaattagccgggtacagtggcaggcacctgtagtcccagctactcgggaggctgaggcaggagaatggcgtgaacccgggaggcagagcttgcagtgagccgagatcgcgccactgcactccagcctgggcgacagagcaagactctgtctcaaaaaaaaaaaaaaaaaaaaaaaaaaGGAAAAGAAAATTAAACAACCAAACAAAATCAGAGTAAATACACCATGTTAATTCTGGTTATATTTGGATTGTGGGCTTATGGGTAGATTTTGTTACATTTTTCTATAATTTCC</sequence>
<map>
<chromosome>chr4</chromosome>
<position>40161303</position>
</map>
</SNP>
<SNP>
<name>rs10000307</name>
<sequence>GTTCAAAGACTCCTGATTAGAGTGTCCTTTCTATAACCAATCTTGTTCCTTAAAACATCTTGAATGATTTGATCTCAGATCCCCTGAAGGGACTGCTGAGATCATCTGCACCAATCCTAAAAAAAAAAATCTTTCATGCCCAAACCCTTAGCAAAGCTAGTTTCTTGTGGGACTCTTAATCCCTCTTATCCTGCTTGACACAGAGGTGCTCACCTGCTGTGCATCAGAAACACTATGGATACTTCTTGAAAGTGCCTGCAACAGAGATACTGATGCATCTGTTGTGGGGTGGGCCCTAGGTATCAGTAATTTAAAAAGTTTTTTAAAATACG[C/T]CCCAGATAATTCTGATTGCTTGTAAATGGCAAAGGTTGAGAAGCACTGCTGGAAGCTTTTGAGCTCCTGTTGGGTAAGTTCAAGCGACAGGAGAATCTCATAGTGATCATAAAACAGCACTCTGAATTCTTGGAGAAACCCAGACTCATCTTATGTGACTAATTTCCTTAATGTGTACCCCAAAACTATCCTAGCGCGTTCACAGGTACACCAGGTAATGCTATTCTGATTGAGCACCCAAGAGTCTC</sequence>
<map>
<chromosome>chr4</chromosome>
<position>188340930</position>
</map>
</SNP>
<SNP>
<name>rs1000031</name>
<sequence>CTTTGAGGATCTCGATGAAAAATCTGCACCTCTCCCAGAAAAATGCACCTCTGCACAGGTTCACAGATGTCTGCATACAATTTCAGGGTTCTCAGACCCTGAAGGCCACCAAGGGACCCAAGTACATGAGCCTTACACAGCACAACCTAAATCGTCAATGGCAATGTCTCAGGAGTGTAGGACAGTGACTGCCTCTGTAAGACCATCAGCACAGCCATGGCCACACATGTTGTCTGGAGGATCAGGTGGCCTTTTTCTGTGGCTTTTGAGGTTGAGGCTGGGTACCCTTGTGGCTAATGCATAATGCCAGGATGGCCAATAAAGACACCATAAAAATTCCCTGCCGTGTGCCTGACACTGGACAGATTTAATCTCCAGGTCTTCTGGGAACCCCGCAgaggcaggggctgttttctcattttactgatggaaactgaggctcaaggaagtgaaggaatttgtttcaagtcccaggcagtacca[C/T]gaacatgggatttgaaatcacgcaagtctgacACGCAAACCTTGGTTCTTTCCTTTTTCCCTTCTCACAGAGGGTGCTTTTCGCTTCCCGGAAGCTGGCAGGGAGTTCCTCTAAAGCGCAGGTTGGAGTGGTCAGAAGGGAGCGAACTGACAGCACGAGGAAGGCTCAGCGCATGCCAGCTCCACTCACGGGAAATGACTCACTGCAGCCCTGCTGCTCTCGGGCTCCGGGGGACACATCCACATTTCCTGTATCTCGGCTAGAGCCTTGGGCAGTGTGAGCTGGCAGGGCAGATCGCTGAAGGCGGCTAGAGATAGAAAACCACCCAGCTCTGCATCCTGAGACAAAGAAGCCTTTCCCTGGGCTCATATGATAGAGGTACGTTGCctctgggcctcagttttgccatctgtaaaatagggTGAAGGTCAGATTAGATTGGGCATATTCAGTGTG</sequence>
<map>
<chromosome>chr18</chromosome>
<position>44615438</position>
</map>
</SNP>
(...)
</MyListOfSnp>

That's it
Pierre

The path from EgonWillighagen to Jandot : Neo4j , a graph API for java: my notebook.

neo4j "Neo4j is a graph database. It is an embedded, disk-based, fully transactional Java persistence engine that stores data structured in graphs rather than in tables. A graph (mathematical lingo for a network) is a flexible data structure that allows a more agile and rapid style of development.".
In the current post, I'll use the neo4j API to load a set of pubmed entries and find the shortest path between two authors.
First, a few constants are defined:

private static final String NODE_TYPE="node-type";
private static final String TYPE_ARTICLE="article";
private static final String TYPE_AUTHOR="author";
private static final String PROPERTY_PMID="pmid";
private static final String PROPERTY_NAME="name";
private static final String PROPERTY_TITLE="title";

We also need a Neo4J embedded graph as well as two indexes that will be used to quickly find if an article or an author have already been inserted in the graph:

//this graph will be stored in /tmp/neo4j
private GraphDatabaseService graphDB =new EmbeddedGraphDatabase("/tmp/neo4j");
private IndexService findByPmid= new LuceneIndexService( graphDB );
private IndexService findByName= new LuceneIndexService( graphDB );

when the program exists, all those objects have to be closed:

private void close()
{
if(findByPmid!=null)
{
findByPmid.shutdown();
findByPmid=null;
}
if(findByName!=null)
{
findByName.shutdown();
findByName=null;
}
if(graphDB!=null)
{
graphDB.shutdown();
}
graphDB=null;
}

A Stax parser is used to read the XML pubmed entries. For this test, I've loaded a bunch of articles from 'Nature' and from the Biogang. Each time a new pubmed identifier (PMID) is found, the Node is searched in the findByPmid index, else a new node is created:

Transaction txn=graphDB.beginTx();
(...)
if(name.equals("PMID") && article==null)
{
Integer pmid=new Integer(reader.getElementText());
article= findByPmid.getSingleNode(PROPERTY_PMID, pmid);

//article was not already in the graph
if(article==null)
{
article= this.graphDB.createNode();
article.setProperty(NODE_TYPE, TYPE_ARTICLE);
article.setProperty(PROPERTY_PMID, pmid);
findByPmid.index(article, PROPERTY_PMID, pmid);
}
}
(..)
txn.success();

Likewise, each time a new pubmed author (PMID) is found, the Node is searched in the findByName index, else a new node is created.

while(!(evt=reader.nextEvent()).isEndDocument())
{
if(evt.isEndElement())
{
if(initials.isEmpty() && !firstName.isEmpty()) initials=""+firstName.charAt(0);
String s= initials+" "+middle+" "+lastName+" "+suffix;
String name= s.replaceAll("[ ]+", " ").trim();
if(name.isEmpty()) return null;
Node author=findByName.getSingleNode(PROPERTY_NAME, name);
if(author==null)
{
author= this.graphDB.createNode();
author.setProperty(NODE_TYPE, TYPE_AUTHOR);
author.setProperty(PROPERTY_NAME, name);
findByName.index(author, PROPERTY_NAME, name);
}
return author;
}
if(!evt.isStartElement()) continue;
String tag= evt.asStartElement().getName().getLocalPart();
String content= reader.getElementText().trim();
if(tag.equals("LastName"))
{
lastName= content;
}
else if(tag.equals("FirstName") || tag.equals("ForeName"))
{
firstName= content;
}
else if(tag.equals("Initials"))
{
initials= content;
}
else if(tag.equals("MiddleName"))
{
middle= content;
}
else if(tag.equals("CollectiveName"))
{
return null;
}
else if(tag.equals("Suffix"))
{
suffix= content;
}

For each author in an Article, a Relationship is created:

author.createRelationshipTo(article, PubmedRelations.IS_AUTHOR_OF);

Now we can find the Path between Jan Aerts and Egon Willighagen (ok, those author could be some homonyms, we need a unique identifier for the Authors ! ).

txn=this.graphDB.beginTx();
Node startAuthor= findByName.getSingleNode(PROPERTY_NAME, "J Aerts");
Node endAuthor= findByName.getSingleNode(PROPERTY_NAME, "EL Willighagen");
SingleSourceShortestPat<Integer> pathBFS;
pathBFS = new SingleSourceShortestPathBFS(
startAuthor,
Direction.BOTH,
PubmedRelations.IS_AUTHOR_OF
);

for(Node n: pathBFS.getPathAsNodes(endAuthor))
{
echo(n);
}

txn.finish();

Result:

1. J Aerts
2. "A physical map of the chicken genome."
3. RK Wilson
4. "Initial sequencing and comparative analysis of the mouse genome"
5. RB Weiss
6. "Independent evolution of bitter-taste sensitivity in humans and chimpanzees"
7. MT Howard
8. "The Blue Obelisk-interoperability in chemical informatics"
9. EL Willighagen

Source code

import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.util.Iterator;
import java.util.List;

import javax.xml.stream.XMLEventReader;
import javax.xml.stream.XMLInputFactory;
import javax.xml.stream.XMLStreamException;
import javax.xml.stream.events.EndElement;
import javax.xml.stream.events.StartElement;
import javax.xml.stream.events.XMLEvent;

import org.neo4j.graphalgo.shortestpath.SingleSourceShortestPath;
import org.neo4j.graphalgo.shortestpath.SingleSourceShortestPathBFS;
import org.neo4j.graphdb.Direction;
import org.neo4j.graphdb.GraphDatabaseService;
import org.neo4j.graphdb.Node;
import org.neo4j.graphdb.RelationshipType;
import org.neo4j.graphdb.Transaction;
import org.neo4j.index.IndexService;
import org.neo4j.index.lucene.LuceneIndexService;
import org.neo4j.kernel.EmbeddedGraphDatabase;

public class PubmedGraph
{
private GraphDatabaseService graphDB;
private IndexService findByPmid;
private IndexService findByName;

private static final String NODE_TYPE="node-type";
private static final String TYPE_ARTICLE="article";
private static final String TYPE_AUTHOR="author";
private static final String PROPERTY_PMID="pmid";
private static final String PROPERTY_NAME="name";
private static final String PROPERTY_TITLE="title";

private static enum PubmedRelations
implements RelationshipType
{
IS_AUTHOR_OF
}



private PubmedGraph()
{

}
private void open()
{
close();
graphDB=new EmbeddedGraphDatabase("/tmp/neo4j");
findByPmid= new LuceneIndexService( graphDB );
findByName= new LuceneIndexService( graphDB );
}
private void close()
{
if(findByPmid!=null)
{
findByPmid.shutdown();
findByPmid=null;
}
if(findByName!=null)
{
findByName.shutdown();
findByName=null;
}
if(graphDB!=null)
{
graphDB.shutdown();
}
graphDB=null;
}

private boolean read(InputStream in)
{
XMLInputFactory xmlInputFactory = XMLInputFactory.newInstance();
xmlInputFactory.setProperty(XMLInputFactory.IS_NAMESPACE_AWARE, Boolean.FALSE);
xmlInputFactory.setProperty(XMLInputFactory.IS_COALESCING, Boolean.TRUE);
xmlInputFactory.setProperty(XMLInputFactory.IS_REPLACING_ENTITY_REFERENCES, Boolean.TRUE);
XMLEventReader reader=null;

Transaction txn=null;
try {
reader= xmlInputFactory.createXMLEventReader(in);
txn = graphDB.beginTx();

Node article=null;
while(reader.hasNext())
{
XMLEvent event = reader.nextEvent();
if(event.isStartElement())
{
StartElement e=event.asStartElement();
String name=e.getName().getLocalPart();
if(name.equals("PubmedArticle"))
{
article=null;
}
else if(name.equals("ArticleTitle") && article!=null && !article.hasProperty(PROPERTY_TITLE))
{
article.setProperty(PROPERTY_TITLE, reader.getElementText());
}
else if(name.equals("PMID") && article==null)
{
Integer pmid=new Integer(reader.getElementText());

article= findByPmid.getSingleNode(PROPERTY_PMID, pmid);


if(article==null)
{
article= this.graphDB.createNode();
article.setProperty(NODE_TYPE, TYPE_ARTICLE);
article.setProperty(PROPERTY_PMID, pmid);
findByPmid.index(article, PROPERTY_PMID, pmid);
}
}
else if(name.equals("Author") && article!=null)
{
Node author= author(reader);
if(author!=null)
{
author.createRelationshipTo(article, PubmedRelations.IS_AUTHOR_OF);
}
}
}
else if(event.isEndElement())
{
EndElement e=event.asEndElement();
String name=e.getName().getLocalPart();
if(name.equals("PubmedArticle"))
{
article=null;
}
}
}
txn.success();
return true;
}
catch (Exception e)
{
e.printStackTrace();
return false;
}
finally
{
if(txn!=null) txn.finish();
if(reader!=null) try {reader.close();}catch(Exception e2){}
}

}
private Node author(XMLEventReader reader) throws IOException,XMLStreamException
{
String lastName="";
String firstName="";
String initials="";
String middle="";
String suffix="";
XMLEvent evt;
while(!(evt=reader.nextEvent()).isEndDocument())
{
if(evt.isEndElement())
{
if(initials.isEmpty() && !firstName.isEmpty()) initials=""+firstName.charAt(0);
String s= initials+" "+middle+" "+lastName+" "+suffix;
String name= s.replaceAll("[ ]+", " ").trim();
if(name.isEmpty()) return null;
Node author=findByName.getSingleNode(PROPERTY_NAME, name);
if(author==null)
{
author= this.graphDB.createNode();
author.setProperty(NODE_TYPE, TYPE_AUTHOR);
author.setProperty(PROPERTY_NAME, name);
findByName.index(author, PROPERTY_NAME, name);
}
return author;
}
if(!evt.isStartElement()) continue;
String tag= evt.asStartElement().getName().getLocalPart();
String content= reader.getElementText().trim();
if(tag.equals("LastName"))
{
lastName= content;
}
else if(tag.equals("FirstName") || tag.equals("ForeName"))
{
firstName= content;
}
else if(tag.equals("Initials"))
{
initials= content;
}
else if(tag.equals("MiddleName"))
{
middle= content;
}
else if(tag.equals("CollectiveName"))
{
return null;
}
else if(tag.equals("Suffix"))
{
suffix= content;
}

else
{
//"###ignoring "+tag+"="+content);
}
}
throw new IOException("Cannot parse Author");
}

public void echo(Node n)
{
System.out.println("nodeId."+n.getId());
for(String p: n.getPropertyKeys())
{
System.out.println(" "+p+":\t"+n.getProperty(p));
}
System.out.println();
}

public void dump()
{
System.out.println("Dump:");
Transaction txn=null;
try {
txn=this.graphDB.beginTx();
for(Node n: this.graphDB.getAllNodes())
{
echo(n);
}
} finally
{
txn.finish();
}

}



private void search()
{
Transaction txn=null;
try {
txn=this.graphDB.beginTx();
Node startAuthor= findByName.getSingleNode(PROPERTY_NAME, "J Aerts");
if(startAuthor==null) { System.err.println("No found;");return;}
Node endAuthor= findByName.getSingleNode(PROPERTY_NAME, "EL Willighagen");

if(endAuthor==null) return;

System.err.println("Waking paths");
SingleSourceShortestPath<Integer> pathBFS;
pathBFS = new SingleSourceShortestPathBFS(
startAuthor,
Direction.BOTH,
PubmedRelations.IS_AUTHOR_OF
);


for(Node n: pathBFS.getPathAsNodes(endAuthor))
{
echo(n);
}

System.err.println("Done");


} finally
{
txn.finish();
}

}


public static void main(String[] args)
{

PubmedGraph app=new PubmedGraph();
try {
int optind=0;

app.open();


while(optind< args.length)
{
InputStream in= new FileInputStream(args[optind++]);
app.read(in);
in.close();
}


app.search();
} catch (Exception e) {
e.printStackTrace();
}
finally
{
if(app!=null) app.close();
}
}
}

That's it
Pierre

Semantic Web Services with the SADI Framework: my notebook.

At Biohackathon2010 , Mark Wilkinson and Luke McCarthy introduced The SADI Framework. From sadiframework.org:SADI is a framework for discovery of, and interoperability between, distributed data and analytical resources. It combines simple, stateless, GET/POST-based Web Services with standards from the W3C Semantic Web initiative. The objective of SADI is to make it easy for data and analytical tool providers to quickly make their resources available on the Semantic Web with minimal disruption to their usual practices.(...)

• SADI Services consume and provide data via simple HTTP POST and GET
• SADI Services consume and produce data in RDF format. This allows SADI Services to exploit existing OWL reasoners and SPARQL query engines to enhance interoperability between Services and the interpretation of the data being passed between them
• Service interfaces (i.e., Inputs and Outputs) are defined in terms of OWL-DL classes; the property restrictions on these OWL classes define what specific data elements are required by the Service and what data will be provided by the Service, respectively
• Input RDF data - data that is compliant with the Input OWL Class - is “decorated” or “annotated” by the service provider to include new properties. These properties will (of course) be a function of the lookup/analytical operations performed by the Web Service.
• Importantly, discovery of SADI Services can include searches for the properties the user wants to add to their data. This contrasts with other Semantic Web Service standards which attempt only to define the computational process by which input data is analysed, rather than the properties that process generates between the input and output data. This is KEY to the semantic behaviours of SADI.
  
• SADI Web Services are stateless and atomic.

In the current post, I just want to understand how SADI invokes the services.

The classical Web Services are described using WSDL (Web Services Description Language) and the messages are transported with SOAP.

Your browser does not support the <CANVAS> element !

In Sadi, as far as I understand it, the description of the services, the operations, the inputs, the ouputs and the transport of the messages use a RDF/OWL format.

Your browser does not support the <CANVAS> element !

http://sadiframework.org/registry/ contains all the services handled by the SADI framework. For example calling http://sadiframework.org/services/getPubMedReferencesForPDB returns the following XML document describing the service:

<rdf:RDF
xmlns="http://www.w3.org/2002/07/owl#"
xmlns:a="http://www.mygrid.org.uk/mygrid-moby-service#"
xmlns:b="http://protege.stanford.edu/plugins/owl/dc/protege-dc.owl#"
xml:base="http://bioinfo.icapture.ubc.ca/SADI"
xmlns:rdfs="http://www.w3.org/2000/01/rdf-schema#"
xmlns:databases="http://sadiframework.org/ontologies/Databases.owl#"
xmlns:misc="http://sadiframework.org/ontologies/miscellaneousObjects.owl#"
xmlns:owl="http://www.w3.org/2002/07/owl#"
xmlns:xsd="http://www.w3.org/2001/XMLSchema#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
<rdf:Description rdf:about="http://sadiframework.org/services/getPubMedReferencesForPDB">
<rdf:type rdf:resource="http://www.mygrid.org.uk/mygrid-moby-service#serviceDescription"/>
<b:format>sadi</b:format>
<b:identifier>urn:lsid:myservices:getPubMedReferencesForPDB</b:identifier>
<a:locationURI>http://sadiframework.org/services/getPubMedReferencesForPDB</a:locationURI>
<a:hasServiceDescriptionText>A implementation of the 'getPubMedReferencesForPDB' service</a:hasServiceDescriptionText>
<a:hasServiceDescriptionLocation>http://sadiframework.org/services/getPubMedReferencesForPDB</a:hasServiceDescriptionLocation>
<a:hasServiceNameText>getPubMedReferencesForPDB</a:hasServiceNameText>
<a:providedBy>
<rdf:Description rdf:about="getPubMedReferencesForPDB_mark.ubic.ca_0">
<a:authoritative>0</a:authoritative>
<b:creator>markw@illuminae.com</b:creator>
<b:publisher>mark.ubic.ca</b:publisher>
<rdf:type rdf:resource="http://www.mygrid.org.uk/mygrid-moby-service#organisation"/>
</rdf:Description>
</a:providedBy>
<a:hasOperation>
<rdf:Description rdf:about="getPubMedReferencesForPDB_mark.ubic.ca_1">
<a:hasOperationNameText>getPubMedReferencesForPDB</a:hasOperationNameText>
<rdf:type rdf:resource="http://www.mygrid.org.uk/mygrid-moby-service#operation"/>
<a:performsTask>
<rdf:Description rdf:about="getPubMedReferencesForPDB_mark.ubic.ca_2">
<rdf:type rdf:resource="http://www.mygrid.org.uk/mygrid-moby-service#operationTask"/>
<rdf:type rdf:resource="http://mygrid.org.uk/sometype"/>
</rdf:Description>
</a:performsTask>
<a:inputParameter>
<rdf:Description rdf:about="getPubMedReferencesForPDB_mark.ubic.ca_3">
<rdf:type rdf:resource="http://www.mygrid.org.uk/mygrid-moby-service#parameter"/>
<a:objectType>
<rdf:Description rdf:about="http://purl.oclc.org/SADI/LSRN/PDB_Thing"/>
</a:objectType>
</rdf:Description>
</a:inputParameter>
<a:outputParameter>
<rdf:Description rdf:about="getPubMedReferencesForPDB_mark.ubic.ca_4">
<rdf:type rdf:resource="http://www.mygrid.org.uk/mygrid-moby-service#parameter"/>
<a:objectType>
<rdf:Description rdf:about="http://sadiframework.org/ontologies/pdb_2_pmid.owl#getPubMedReferencesForPDB_Output"/>
</a:objectType>
</rdf:Description>
</a:outputParameter>
</rdf:Description>
</a:hasOperation>
</rdf:Description>
</rdf:RDF>

A part of this service can be visualized as a graph:

Your browser does not support the <CANVAS> element !

So, the service getPubMedReferencesForPDB takes as input a parameter of type http://purl.oclc.org/SADI/LSRN/PDB_Thing and returns an object of type http://sadiframework.org/ontologies/pdb_2_pmid.owl#getPubMedReferencesForPDB_Output. Ok, let's try this service: I'm going to invoke the service getPubMedReferencesForPDB with two structures of type http://purl.oclc.org/SADI/LSRN/PDB_Thing (Many thanks to Mark for helping me via Twitter ...:-P ): 1KNZ and 1LJ2.

File: input.rdf

<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:lsrn="http://purl.oclc.org/SADI/LSRN/">
<lsrn:PDB_Thing rdf:about="http://lsrn.org/PDB:1KNZ"/>
<lsrn:PDB_Thing rdf:about="http://lsrn.org/PDB:1LJ2"/>
</rdf:RDF>

Calling curl with this file for the service getPubMedReferencesForPDB:

curl -d @input.rdf http://sadiframework.org/services/getPubMedReferencesForPDB

Result:

<rdf:RDF>
<rdf:Description rdf:about="http://lsrn.org/PDB:1LJ2">
<a:hasReference>
<rdf:Description rdf:about="http://lsrn.org/PMID:12086624">
<rdf:type rdf:resource="http://purl.oclc.org/SADI/LSRN/PMID_Thing"/>
</rdf:Description>
</a:hasReference>
<rdf:type rdf:resource="http://sadiframework.org/ontologies/pdb_2_pmid.owl#getPubMedReferencesForPDB_Output"/>
</rdf:Description>
<rdf:Description rdf:about="http://lsrn.org/PDB:1KNZ">
<a:hasReference>
<rdf:Description rdf:about="http://lsrn.org/PMID:11792322">
<rdf:type rdf:resource="http://purl.oclc.org/SADI/LSRN/PMID_Thing"/>
</rdf:Description>
</a:hasReference>
<rdf:type rdf:resource="http://sadiframework.org/ontologies/pdb_2_pmid.owl#getPubMedReferencesForPDB_Output"/>
</rdf:Description>
</rdf:RDF>

Here, for each PDB entry (input), we were returned a Pubmed Identifier (output).


That's it !
Pierre