Showing posts with label sqlite. Show all posts
Showing posts with label sqlite. Show all posts

21 May 2016

Playing with the @ORCID_Org / @ncbi_pubmed graph. My notebook.

"ORCID provides a persistent digital identifier that distinguishes you from every other researcher and, through integration in key research workflows such as manuscript and grant submission, supports automated linkages between you and your professional activities ensuring that your work is recognized. "
I've recently discovered that pubmed now integrates ORCID identfiers.

And there are several minor problems, I found some articles where the ORCID id is malformed or where different people use the same ORCID-ID:







You can download the papers containing some orcid Identifiers using the entrez query http://www.ncbi.nlm.nih.gov/pubmed/?term=orcid[AUID].
I've used one of my tools pubmeddump to download the articles asXML and I wrote PubmedOrcidGraph to extract the author's orcid.
<?xml version="1.0" encoding="UTF-8"?>
<PubmedArticleSet>
  <!--Generated with PubmedOrcidGraph https://github.com/lindenb/jvarkit/wiki/PubmedOrcidGraph - Pierre Lindenbaum.-->
  <PubmedArticle pmid="27197243" doi="10.1101/gr.199760.115">
    <year>2016</year>
    <journal>Genome Res.</journal>
    <title>Improved definition of the mouse transcriptome via targeted RNA sequencing.</title>
    <Author orcid="0000-0002-4078-7413">
      <foreName>Giovanni</foreName>
      <lastName>Bussotti</lastName>
      <initials>G</initials>
      <affiliation>EMBL, European Bioinformatics Institute, Cambridge, CB10 1SD, United Kingdom;</affiliation>
    </Author>
    <Author orcid="0000-0002-4449-1863">
      <foreName>Tommaso</foreName>
      <lastName>Leonardi</lastName>
      <initials>T</initials>
      <affiliation>EMBL, European Bioinformatics Institute, Cambridge, CB10 1SD, United Kingdom;</affiliation>
    </Author>
    <Author orcid="0000-0002-6090-3100">
      <foreName>Anton J</foreName>
      <lastName>Enright</lastName>
      <initials>AJ</initials>
      <affiliation>EMBL, European Bioinformatics Institute, Cambridge, CB10 1SD, United Kingdom;</affiliation>
    </Author>
  </PubmedArticle>
  <PubmedArticle pmid="27197225" doi="10.1101/gr.204479.116">
    <year>2016</year>
    <journal>Genome Res.</journal>
(...)
Now, I want to insert those data into a sqlite3 database. I use the XSLT stylesheet below to convert the XML into some SQL statement.
<?xml version="1.0"?>
<xsl:stylesheet
 xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
 version="1.0"
    xmlns:xalan="http://xml.apache.org/xalan"
    xmlns:str="xalan://com.github.lindenb.xslt.strings.Strings"
    exclude-result-prefixes="xalan str"
 >
<xsl:output method="text"/>
<xsl:variable name="q">'</xsl:variable>

<xsl:template match="/">
create table author(orcid text unique,name text,affiliation text);
create table collab(orcid1 text,orcid2 text,unique(orcid1,orcid2));
begin transaction;
<xsl:apply-templates select="PubmedArticleSet/PubmedArticle"/>
commit;
</xsl:template>

<xsl:template match="PubmedArticle">
<xsl:for-each select="Author">
<xsl:variable name="o1" select="@orcid"/>insert or ignore into author(orcid,name,affiliation) values ('<xsl:value-of select="$o1"/>','<xsl:value-of select="translate(concat(lastName,' ',foreName),$q,' ')"/>','<xsl:value-of select="translate(affiliation,$q,' ')"/>');
<xsl:for-each select="following-sibling::Author">insert or ignore into collab(orcid1,orcid2) values(<xsl:variable name="o2" select="@orcid"/>
<xsl:choose>
 <xsl:when test="str:strcmp( $o1 , $o2) < 0">'<xsl:value-of select='$o1'/>','<xsl:value-of select='$o2'/>'</xsl:when>
 <xsl:otherwise>'<xsl:value-of select='$o2'/>','<xsl:value-of select='$o1'/>'</xsl:otherwise>
</xsl:choose>);
</xsl:for-each>
</xsl:for-each>
</xsl:template>
</xsl:stylesheet>

This stylesheet contains an extension 'strmcp' for the xslt processor xalan to compare two XML strings
This extension is just used to always be sure that the field "orcid1" in the table "collab" is always lower than "orcid2" to avoid duplicates pairs.
./src/xslt-sandbox/xalan/dist/xalan -XSL orcid2sqlite.xsl -IN orcid.xml

create table author(orcid text unique,name text,affiliation text);
create table collab(orcid1 text,orcid2 text,unique(orcid1,orcid2));
begin transaction;
insert or ignore into author(orcid,name,affiliation) values ('0000-0002-4078-7413','Bussotti Giovanni','EMBL, European Bioinformatics Institute, Cambridge, CB10 1SD, United Kingdom;');
insert or ignore into collab(orcid1,orcid2) values('0000-0002-4078-7413','0000-0002-4449-1863');
insert or ignore into collab(orcid1,orcid2) values('0000-0002-4078-7413','0000-0002-6090-3100');
insert or ignore into author(orcid,name,affiliation) values ('0000-0002-4449-1863','Leonardi Tommaso','EMBL, European Bioinformatics Institute, Cambridge, CB10 1SD, United Kingdom;');
insert or ignore into collab(orcid1,orcid2) values('0000-0002-4449-1863','0000-0002-6090-3100');
insert or ignore into author(orcid,name,affiliation) values ('0000-0002-6090-3100','Enright Anton J','EMBL, European Bioinformatics Institute, Cambridge, CB10 1SD, United Kingdom;');
(...)
and those sql statetements are loaded into sqlite3:
./src/xslt-sandbox/xalan/dist/xalan -XSL orcid2sqlite.xsl -IN orcid.xml |\
 sqlite3 orcid.sqlite

The next step is to produce a gexf+xml file to play with the orcid graph in gephi.
I use the following bash script to convert the sqlite3 database to gexf+xml.
DB=orcid.sqlite

cat << EOF
<?xml version="1.0" encoding="UTF-8"?>
<gexf xmlns="http://www.gexf.net/1.2draft" xmlns:viz="http://www.gexf.net/1.1draft/viz" version="1.2">
<meta>
<creator>Pierre Lindenbaum</creator>
<description>Orcid Graph</description>
</meta>
<graph defaultedgetype="undirected" mode="static">

<attributes class="node">
<attribute type="string" title="affiliation" id="0"/>
</attributes>
<nodes>
EOF

sqlite3 -separator ' ' -noheader  ${DB} 'select orcid,name,affiliation from author' |\
 sed  -e 's/&/&/g' -e "s/</\</g" -e "s/>/\>/g" -e "s/'/\'/g"  -e 's/"/\"/g' |\
 awk -F ' ' '{printf("<node id=\"%s\" label=\"%s\"><attvalue for=\"0\" value=\"%s\"/></node>\n",$1,$2,$3);}'

echo "</nodes><edges>"
sqlite3 -separator ' ' -noheader  ${DB} 'select orcid1,orcid2 from collab' |\
 awk -F ' ' '{printf("<edge source=\"%s\" target=\"%s\"/>\n",$1,$2);}'
echo "</edges></graph></gexf>"



The output is saved and then loaded into gephi.






That's it,

Pierre

01 May 2013

Inserting the result of a BLAST into a Database using XSLT.

Here is the XML output of a BLAST:

<?xml version="1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "NCBI_BlastOutput.dtd">
<BlastOutput>
  <BlastOutput_program>tblastn</BlastOutput_program>
  <BlastOutput_version>TBLASTN 2.2.27+</BlastOutput_version>
  <BlastOutput_reference>Stephen F. Altschul, Thomas L. Madden, Alejandro A. Sch&amp;auml;ffer, Jinghui Z
hang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), &quot;Gapped BLAST and PSI-BLAST: a new generation 
of protein database search programs&quot;, Nucleic Acids Res. 25:3389-3402.</BlastOutput_reference>
  <BlastOutput_db>nr</BlastOutput_db>
  <BlastOutput_query-ID>52385</BlastOutput_query-ID>
  <BlastOutput_query-def>myseq</BlastOutput_query-def>
  <BlastOutput_query-len>30</BlastOutput_query-len>
  <BlastOutput_param>
    <Parameters>
      <Parameters_matrix>BLOSUM62</Parameters_matrix>
      <Parameters_expect>10</Parameters_expect>
      <Parameters_gap-open>11</Parameters_gap-open>
      <Parameters_gap-extend>1</Parameters_gap-extend>
      <Parameters_filter>L;</Parameters_filter>
    </Parameters>
  </BlastOutput_param>
<BlastOutput_iterations>
<Iteration>
  <Iteration_iter-num>1</Iteration_iter-num>
  <Iteration_query-ID>52385</Iteration_query-ID>
  <Iteration_query-def>myseq</Iteration_query-def>
  <Iteration_query-len>30</Iteration_query-len>
<Iteration_hits>
<Hit>
  <Hit_num>1</Hit_num>
  <Hit_id>gi|110624327|dbj|AK225891.1|</Hit_id>
  <Hit_def>Homo sapiens mRNA for zinc finger CCCH-type containing 7B variant, clone: FCC121C02</Hit_def>
  <Hit_accession>AK225891</Hit_accession>
  <Hit_len>1829</Hit_len>
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>62.3882</Hsp_bit-score>
      <Hsp_score>150</Hsp_score>
      <Hsp_evalue>4.01658e-10</Hsp_evalue>
      <Hsp_query-from>1</Hsp_query-from>
      <Hsp_query-to>30</Hsp_query-to>
      <Hsp_hit-from>250</Hsp_hit-from>
      <Hsp_hit-to>339</Hsp_hit-to>
      <Hsp_query-frame>0</Hsp_query-frame>
      <Hsp_hit-frame>1</Hsp_hit-frame>
      <Hsp_identity>30</Hsp_identity>
      <Hsp_positive>30</Hsp_positive>
      <Hsp_gaps>0</Hsp_gaps>
      <Hsp_align-len>30</Hsp_align-len>
      <Hsp_qseq>MERQKRKADIEKGLQFIQSTLPLKQEEYEA</Hsp_qseq>
      <Hsp_hseq>MERQKRKADIEKGLQFIQSTLPLKQEEYEA</Hsp_hseq>
      <Hsp_midline>MERQKRKADIEKGLQFIQSTLPLKQEEYEA</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>
<Hit>
  <Hit_num>2</Hit_num>
  <Hit_id>gi|6176337|gb|AF188530.1|AF188530</Hit_id>
  <Hit_def>Homo sapiens ubiquitous tetratricopeptide containing protein RoXaN mRNA, partial cds</Hit_def&g
t;
  <Hit_accession>AF188530</Hit_accession>
  <Hit_len>2398</Hit_len>
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>62.3882</Hsp_bit-score>
      <Hsp_score>150</Hsp_score>
      <Hsp_evalue>4.12279e-10</Hsp_evalue>
      <Hsp_query-from>1</Hsp_query-from>
      <Hsp_query-to>30</Hsp_query-to>
      <Hsp_hit-from>105</Hsp_hit-from>
      <Hsp_hit-to>194</Hsp_hit-to>
      <Hsp_query-frame>0</Hsp_query-frame>
      <Hsp_hit-frame>3</Hsp_hit-frame>
      <Hsp_identity>30</Hsp_identity>
      <Hsp_positive>30</Hsp_positive>
      <Hsp_gaps>0</Hsp_gaps>
      <Hsp_align-len>30</Hsp_align-len>
      <Hsp_qseq>MERQKRKADIEKGLQFIQSTLPLKQEEYEA</Hsp_qseq>
      <Hsp_hseq>MERQKRKADIEKGLQFIQSTLPLKQEEYEA</Hsp_hseq>
      <Hsp_midline>MERQKRKADIEKGLQFIQSTLPLKQEEYEA</Hsp_midline>
    </Hsp>
  </Hit_hsps>
(...)
We want to insert that XML file into a database. I wrote the following XSLT stylesheet , it transforms the blast-xml into a set of SQL statements for sqlite3.
xsltproc --novalid blast2sqlite.xsl blast.xml 

create table BlastOutput(
(...)


BEGIN TRANSACTION;

insert into BlastOutput(
 program,
 version,
 reference,
 db,
 query_ID,
 query_def,
 query_len
 )
 values (
  'tblastn',
  'TBLASTN 2.2.27+',
  'Stephen F. Altschul, Thomas L. Madden, Alejandro A. Sch&auml;ffer, Jinghui Zhang,Zheng Zhang, Webb Miller, and David J. Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of protein database search programs", Nucleic Acids Res. 25:3389-3402.',
  'nr',
  '52385',
  'myseq',
  30
  );


insert into Parameters(
 blastOutput_id,
 expect,
 matrix,
 sc_match,
 sc_mismatch,
 gap_open,
 gap_extend,
 filter
 )
select MAX(id),
 10,
 'BLOSUM62',
 NULL,
 NULL,
 11,
 1,
 'L;'
 from BlastOutput;
 


insert into Iteration(
 blastOutput_id,
 iter_num,
 query_id,
 query_def,
 query_len
 )
select MAX(id),
 1,
 '52385',
 'myseq',
 30
from BlastOutput;

insert into Hit(iteration_id,num,hit_id,def,accession,len)
select MAX(id),
 1,
 'gi|110624327|dbj|AK225891.1|',
 'Homo sapiens mRNA for zinc finger CCCH-type containing 7B variant, clone: FCC121C02',
 'AK225891',
 1829
from Iteration;
(...)
All in one you can redirect the output to sqlite3.
xsltproc --novalid blast2sqlite.xsl blast.xml |\
 sqlite3 input.db
and query the database:
$ sqlite3 -header -line input.sqlite \
 'select * from Hsp,Hit where Hsp.hit_id=Hit.id limit 2'

          id = 1
      hit_id = 1
         num = 1
   bit_score = 62.3882
       score = 150.0
      evalue = 4.01658e-10
  query_from = 1
    query_to = 30
    hit_from = 250
      hit_to = 339
 query_frame = 0
   hit_frame = 1
    identity = 30
    positive = 30
        gaps = 0
   align_len = 30
        qseq = MERQKRKADIEKGLQFIQSTLPLKQEEYEA
        hseq = MERQKRKADIEKGLQFIQSTLPLKQEEYEA
     midline = MERQKRKADIEKGLQFIQSTLPLKQEEYEA
          id = 1
iteration_id = 1
         num = 1
      hit_id = gi|110624327|dbj|AK225891.1|
         def = Homo sapiens mRNA for zinc finger CCCH-type containing 7B variant, clone: FCC121C02
   accession = AK225891
         len = 1829

          id = 2
      hit_id = 2
         num = 1
   bit_score = 62.3882
       score = 150.0
      evalue = 4.12279e-10
  query_from = 1
    query_to = 30
    hit_from = 105
      hit_to = 194
 query_frame = 0
   hit_frame = 3
    identity = 30
    positive = 30
        gaps = 0
   align_len = 30
        qseq = MERQKRKADIEKGLQFIQSTLPLKQEEYEA
        hseq = MERQKRKADIEKGLQFIQSTLPLKQEEYEA
     midline = MERQKRKADIEKGLQFIQSTLPLKQEEYEA
          id = 2
iteration_id = 1
         num = 2
      hit_id = gi|6176337|gb|AF188530.1|AF188530
         def = Homo sapiens ubiquitous tetratricopeptide containing protein RoXaN mRNA, partial cds
   accession = AF188530
         len = 2398
That's it,

Pierre

02 November 2012

Saving your tweets in a database using sqlite, rhino, scribe, javascript

In the current post, I 'll describe a simple method to save your tweets in a sqlite database using Mozilla Rhino.

Prerequisites

  • sqlite
  • Apache Rhino. I think it should be de-facto available when the java developer toolkit (JDK) is installed
  • Scribe, the simple OAuth library for Java . It also requires Apache codec

The config.js file

Open an account on https://dev.twitter.com/ and create an App to receive an API-key and an API-secret.
Create the following file 'config.js' filled with the correct parameters.

The javascript

The following javascript file opens a Oauth connection, retrieves the tweets and stores them into sqlite. I've commented the code, I hope it is clear enough.

Running the script using Rhino

scribe.libs=/path/to/scribe-1.3.2.jar:/path/to/commons-codec.jar
rhino.libs=/usr/share/java/js.jar:/usr/share/java/jline.jar
sqlite.libs=/path/to/sqlitejdbc-v056.jar
CLASSPATH=${rhino.libs}:${scribe.libs}:${sqlite.libs}

java -cp ${CLASSPATH} org.mozilla.javascript.tools.shell.Main -f twitter2sqlite.js
At the first time, the user is asked to authorize the application to use the twitter API

The script runs forever (Ctrl-C to break), listening to the new tweets.

As a test, I wrote the following tweet:


... and the tweet was later inserted in the database...

Sleep...

Inserted ({created_at:"Fri Nov 02 20:29:04 +0000 2012", id:264464160664981500, id_str:"264464160664981504", text:"wrote a tool to save my tweets: This is a test . ( #rhino, #jdbc, #sqlite, #scribe #javascript )", source:"web", truncated:false, in_reply_to_status_id:null, in_reply_to_status_id_str:null, in_reply_to_user_id:null, in_reply_to_user_id_str:null, in_reply_to_screen_name:null, geo:null, coordinates:null, place:null, contributors:null, retweet_count:0, entities:{hashtags:[{text:"rhino", indices:[51, 57]}, {text:"jdbc", indices:[59, 64]}, {text:"sqlite", indices:[66, 73]}, {text:"scribe", indices:[75, 82]}, {text:"javascript", indices:[83, 94]}], urls:[], user_mentions:[]}, favorited:false, retweeted:false})

Sleep...
Sleep...
Sleep...

Later, the tweets can be extracted using the sqlite command line:

$  sqlite3 tweets.sqlite 'select * from tweet'

264464160664981504|({created_at:"Fri Nov 02 20:29:04 +0000 2012", id:264464160664981500, id_str:"264464160664981504", text:"wrote a tool to save my tweets: This
264421310841638913|({created_at:"Fri Nov 02 17:38:47 +0000 2012", id:264421310841638900, id_str:"264421310841638913", text:"The tools for recalibration have cha
264264932097400832|({created_at:"Fri Nov 02 07:17:24 +0000 2012", id:264264932097400830, id_str:"264264932097400832", text:"@warandpeace you're welcome. Your sh
264158323287416832|({created_at:"Fri Nov 02 00:13:46 +0000 2012", id:264158323287416830, id_str:"264158323287416832", text:"Drawing of the day November 1, 2012.
264142732174438400|({created_at:"Thu Nov 01 23:11:49 +0000 2012", id:264142732174438400, id_str:"264142732174438400", text:"[delicious] PLOS Collections: How th
264064117558624256|({created_at:"Thu Nov 01 17:59:26 +0000 2012", id:264064117558624260, id_str:"264064117558624256", text:"I've added a stupid basic dependency
264025607724204034|({created_at:"Thu Nov 01 15:26:24 +0000 2012", id:264025607724204030, id_str:"264025607724204034", text:"in the desert lab, checking my on-go
264013563704795136|({created_at:"Thu Nov 01 14:38:33 +0000 2012", id:264013563704795140, id_str:"264013563704795136", text:"Drawing of the day November 1, 2012.
263996436679630848|({created_at:"Thu Nov 01 13:30:29 +0000 2012", id:263996436679630850, id_str:"263996436679630848", text:"RT @RealistComics: he's tall, dark a
263966759210590208|({created_at:"Thu Nov 01 11:32:34 +0000 2012", id:263966759210590200, id_str:"263966759210590208", text:"RT @guermonprez: #Aubry Un avion nor
263946369847398402|({created_at:"Thu Nov 01 10:11:33 +0000 2012", id:263946369847398400, id_str:"263946369847398402", text:"[delicious] OVal: object validation 
263946366919790593|({created_at:"Thu Nov 01 10:11:32 +0000 2012", id:263946366919790600, id_str:"263946366919790593", text:"[delicious] MyBatis #tweet: a first 
263941020729896960|({created_at:"Thu Nov 01 09:50:17 +0000 2012", id:263941020729896960, id_str:"263941020729896960", text:"RT @josh_wills: I have never been pr
263938670187388928|({created_at:"Thu Nov 01 09:40:57 +0000 2012", id:263938670187388930, id_str:"263938670187388928", text:"RT @softmodeling @peterneubauer: Usi
263936362716200960|({created_at:"Thu Nov 01 09:31:47 +0000 2012", id:263936362716200960, id_str:"263936362716200960", text:"declined to review an article about 
263934528186351616|({created_at:"Thu Nov 01 09:24:29 +0000 2012", id:263934528186351600, id_str:"263934528186351616", text:"@figshare Thanks, ( was http://t.co/
263815846139412480|({created_at:"Thu Nov 01 01:32:53 +0000 2012", id:263815846139412480, id_str:"263815846139412480", text:"Drawing of the day October 30, 2012.
263731855919026176|({created_at:"Wed Oct 31 19:59:09 +0000 2012", id:263731855919026180, id_str:"263731855919026176", text:"[delicious] An integrated map of gen
263726281647067136|({created_at:"Wed Oct 31 19:36:59 +0000 2012", id:263726281647067140, id_str:"263726281647067136", text:"RT @bryan_howie: 1000 Genomes paper 
263695076516052992|({created_at:"Wed Oct 31 17:33:00 +0000 2012", id:263695076516053000, id_str:"263695076516052992", text:"\"Forget your Past\" ( abandoned Bul

That's it
Pierre

14 October 2012

Calculating time from submission to publication / Degree of burden in submitting a paper

After "404 not found": a database of non-functional resources in the NAR database collection, I've uploaded my second dataset on figshare:
Calculating time from submission to publication / Degree of burden in submitting a paper
.

Calculating time from submission to publication / Degree of burden in submitting a paper. Pierre Lindenbaum,  Ryan Delahanty.
figshare.
Retrieved 10:13, Oct 14, 2012 (GMT)
http://dx.doi.org/10.6084/m9.figshare.96403

This dataset was inspired by this post on biostar, initialy asked by Ryan Delahanty: I was wondering if it would be possible to calculate some kind of a metric for the speed-of-publication for each journal. I'm not sure submitted and accepted dates are available for all papers, but I noticed in XML data there are fields like the following:
<PubmedData>
        <History>
            <PubMedPubDate PubStatus="received">
                <Year>2011</Year>
                <Month>11</Month>
                <Day>29</Day>
                <Hour>6</Hour>
                <Minute>0</Minute>
            </PubMedPubDate>
            <PubMedPubDate PubStatus="accepted">
                <Year>2011</Year>
                <Month>12</Month>
                <Day>20</Day>
                <Hour>6</Hour>
                <Minute>0</Minute>
            </PubMedPubDate>
           (...)

In this dataset, the script 'pubmed.sh" downloads the the journals from http://www.ncbi.nlm.nih.gov/books/NBK3827/table/pubmedhelp.pubmedhelptable45/ , the 'eigenfactors' from http://www.eigenfactor.org.

For each journal , It scans pubmed (starting from year=2000) and get the difference between the date[@PubStatus='received'] and the date[@PubStatus='accepted'].

titleissneigenfactordays
"Acta biochimica Polonica"0001-527X0.003996119.770935960591
"Acta biomaterialia"1742-70610.02152129.682692307692
"Acta biotheoretica"0001-53420.000844161.897058823529
"Acta cirurgica brasileira / Sociedade Brasileira para Desenvolvimento Pesquisa em Cirurgia"0102-86500.00128122.038461538462
"Acta cytologica"0001-55470.00230565.3006134969325
"Acta diabetologica"0940-54290.001851299.6
"Acta haematologica"0001-57920.002825118.654676258993
"Acta histochemica"0065-12810.002162110.471204188482
"Acta histochemica et cytochemica"0044-59910.00067781.6455696202532
"Acta neurochirurgica"0001-62680.009685204.371830985916
"Acta neuropathologica"0001-63220.02347169.7277882797732
"Acta theriologica"0001-70510.000901147.0
"Acta tropica"0001-706X0.01011196.577777777778
"Acta veterinaria Scandinavica"0044-605X0.00161282.0
"Addictive behaviors"0306-46030.017915163.049731182796
"Advances in space research "0273-11770.021217205.0
Ambio0044-74470.007463181.878048780488
"American journal of human genetics"0002-92970.12015667.1898928024502
"American journal of hypertension"0895-70610.017359104.074576271186
(....)

Here is the kind of figure I got:

As far as I remember, "Cell" is the point having the highest eigenfactor.


Note: pubmed contains some errors: e.g. received > accepted (http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&id=20591334&retmode=xml) or some dates in the future: ( http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&id=12921703&retmode=xml )


That's it,

Pierre

08 May 2012

Downloading the XML data from the Exome Variant Server

From EVS: "The goal of the NHLBI GO Exome Sequencing Project (ESP) is to discover novel genes and mechanisms contributing to heart, lung and blood disorders by pioneering the application of next-generation sequencing of the protein coding regions of the human genome across diverse, richly-phenotyped populations and to share these datasets and findings with the scientific community to extend and enrich the diagnosis, management and treatment of heart, lung and blood disorders."

The NHLBI Exome Sequencing Project provides a download area but I wanted to build a local database for the richer XML data returned by their Web Services (previously described here on my blog ). The following java program sends some XML/SOAP requests to the EVS server for each chromosome using a genomic window of 150000 bp and parses the XML response.



it then dumps the results into a tab-delimited file:

Compilation:

javac DumpExomeVariantServerData.java

Execution:
java DumpExomeVariantServerData > input.evs.tsv

May 8, 2012 7:59:58 AM sandbox.DumpExomeVariantServerData fetchEvsData
INFO: 1:1-200011
May 8, 2012 8:00:02 AM sandbox.DumpExomeVariantServerData fetchEvsData
INFO: 1:200001-400011
May 8, 2012 8:00:03 AM sandbox.DumpExomeVariantServerData fetchEvsData
INFO: 1:400001-600011
May 8, 2012 8:00:04 AM sandbox.DumpExomeVariantServerData fetchEvsData
INFO: 1:600001-800011
May 8, 2012 8:00:05 AM sandbox.DumpExomeVariantServerData fetchEvsData
INFO: 1:800001-1000011
(...)

head input.evs.tsv
1 69116 <snpList><positionString>1:69116</positionString><chrPosition>69116</chr
1 69134 <snpList><positionString>1:69134</positionString><chrPosition>69134</chr
1 69270 <snpList><positionString>1:69270</positionString><chrPosition>69270</chr
1 69428 <snpList><positionString>1:69428</positionString><chrPosition>69428</chr
1 69453 <snpList><positionString>1:69453</positionString><chrPosition>69453</chr
1 69476 <snpList><positionString>1:69476</positionString><chrPosition>69476</chr
1 69496 <snpList><positionString>1:69496</positionString><chrPosition>69496</chr
1 69511 <snpList><positionString>1:69511</positionString><chrPosition>69511</chr
1 69552 <snpList><positionString>1:69552</positionString><chrPosition>69552</chr
1 69590 <snpList><positionString>1:69590</positionString><chrPosition>69590</chr

Inserting the EVS data in a sqlite database

We can now create a sqlite3 database to insert the data ...
$ sqlite3 evs.sqlite
sqlite> create table evsData(chrom TEXT NOT NULL,pos INT NOT NULL,xml TEXT NOT NULL);
sqlite> create index chrompos on evsData(chrom,pos);
sqlite> .separator "\t";
sqlite> .import "input.evs.tsv" evsData

... and query this database
$ sqlite3 evs.sqlite 'select xml from evsData where chrom="1" and pos=69552' |\
xmllint  --format -

<?xml version="1.0"?>
<snpList>
  <positionString>1:69552</positionString>
  <chrPosition>69552</chrPosition>
  <alleles>C/G</alleles>
  <uaAlleleCounts>C=4/G=4644</uaAlleleCounts>
  <aaAlleleCounts>C=0/G=2944</aaAlleleCounts>
  <totalAlleleCounts>C=4/G=7588</totalAlleleCounts>
  <uaMAF>0.0861</uaMAF>
  <aaMAF>0.0</aaMAF>
  <totalMAF>0.0527</totalMAF>
  <avgSampleReadDepth>143</avgSampleReadDepth>
  <geneList>OR4F5</geneList>
  <snpFunction>
    <chromosome>1</chromosome>
    <position>69552</position>
    <conservationScore>1.0</conservationScore>
    <conservationScoreGERP>-0.1</conservationScoreGERP>
    <snpFxnList>
      <mrnaAccession>NM_001005484.1</mrnaAccession>
      <fxnClassGVS>coding-synonymous</fxnClassGVS>
      <aminoAcids>none</aminoAcids>
      <proteinPos>154/306</proteinPos>
      <cdnaPos>462</cdnaPos>
      <pphPrediction>unknown</pphPrediction>
      <granthamScore>NA</granthamScore>
    </snpFxnList>
    <refAllele>G</refAllele>
    <ancestralAllele>C</ancestralAllele>
    <firstRsId>55874132</firstRsId>
    <secondRsId>0</secondRsId>
    <filters>SVM</filters>
    <clinicalLink>unknown</clinicalLink>
  </snpFunction>
  <conservationScore>1.0</conservationScore>
  <conservationScoreGERP>-0.1</conservationScoreGERP>
  <refAllele>G</refAllele>
  <altAlleles>C</altAlleles>
  <ancestralAllele>C</ancestralAllele>
  <chromosome>1</chromosome>
  <hasAtLeastOneAccession>true</hasAtLeastOneAccession>
  <rsIds>rs55874132</rsIds>
  <filters>SVM</filters>
  <clinicalLink>unknown</clinicalLink>
  <dbsnpVersion>dbSNP_129</dbsnpVersion>
  <uaGenotypeCounts>CC=0/CG=4/GG=2320</uaGenotypeCounts>
  <aaGenotypeCounts>CC=0/CG=0/GG=1472</aaGenotypeCounts>
  <totalGenotypeCounts>CC=0/CG=4/GG=3792</totalGenotypeCounts>
  <onExomeChip>false</onExomeChip>
  <gwasPubmedIds>unknown</gwasPubmedIds>
</snpList>

The Variation Toolkit

I also wrote a C++ program (that is part of my (always-beta) Variation Toolkit) to use this sqlite database to annotate some VCF-like files. See http://code.google.com/p/variationtoolkit/wiki/VcfEvs

Example 1

$ echo -e "#CHROM\tPOS\n1\t69511\n1\t69512\n1\t69552" |\
vcfevs -f evs.sqlite 

#CHROM  POS     evs.positionString      evs.chrPosition evs.alleles     evs.uaAlleleCounts      evs.aaAlleleCounts      evs.totalAlleleCounts   evs.uaMAF       evs.aaMAF       evs.totalMAF    evs.avgSampleReadDepth  evs.geneList    evs.conservationScore   evs.conservationScoreGERP       evs.refAllele   evs.altAllelesevs.ancestralAllele       evs.chromosome  evs.hasAtLeastOneAccession      evs.rsIds       evs.filters     evs.clinicalLink        evs.dbsnpVersion        evs.uaGenotypeCounts    evs.aaGenotypeCounts    evs.totalGenotypeCounts evs.onExomeChip evs.gwasPubmedIds
1       69511   1:69511 69511   G/A     G=4235/A=483    G=1707/A=1297   G=5942/A=1780   10.2374 43.1758 23.051  74      OR4F5   1.0     1.1     A       G       G1      true    rs75062661      PASS    unknown dbSNP_131       GG=1964/GA=307/AA=88    GG=703/GA=301/AA=498    GG=2667/GA=608/AA=586   false   unknown
1       69512   .       .       .       .       .       .       .       .       .       .       .       .       .       .       .       .       .       .       ..      .       .       .       .       .       .       .
1       69552   1:69552 69552   C/G     C=4/G=4644      C=0/G=2944      C=4/G=7588      0.0861  0.0     0.0527  143     OR4F5   1.0     -0.1    G       C       C1      true    rs55874132      SVM     unknown dbSNP_129       CC=0/CG=4/GG=2320       CC=0/CG=0/GG=1472       CC=0/CG=4/GG=3792       false   unknown

Example 2

$ echo -e "#CHROM\tPOS\n1\t69511\n1\t69512\n1\t69552" |\
vcfevs -f ~/WORK/20120506.evs.download/evs.sqlite -c uaMAF 
#CHROM  POS     evs.uaMAF
1       69511   10.2374
1       69512   .
1       69552   0.0861

Example 3

$ echo -e "#CHROM\tPOS\n1\t69511\n1\t69512\n1\t69552" |\
vcfevs -f evs.sqlite  -x -c _ | cut -c 1-200

#CHROM  POS     evs.xml
1       69511   <snpList><positionString>1:69511</positionString><chrPosition>69511</chrPosition><alleles>G/A</alleles><uaAlleleCounts>G=4235/A=483</uaAlleleCounts><aaAlleleCounts>G=1707/A=1297</aaAlleleCount
1       69512   .
1       69552   <snpList><positionString>1:69552</positionString><chrPosition>69552</chrPosition><alleles>C/G</alleles><uaAlleleCounts>C=4/G=4644</uaAlleleCounts><aaAlleleCounts>C=0/G=2944</aaAlleleCounts><to

That's it,

Pierre

28 January 2012

Insert your VCFs in a sqlite database.

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

Download

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

Dependencies

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

Compilation

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

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

Usage

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

Options

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

Schema


Example:

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

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

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

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

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

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

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

11 May 2010

Playing with the C API for sqlite: inserting the hapmap data , my notebook.



In this post, I'll insert a Hapmap genotype file into a SQLite database using the C API for SQLIte. Of course, this is just a test and I would not use SQLite to store this kind of information. Moreover the speed of insertion of SQLite was deadly slow here.

The Genotype file

The first line contains the name of the individuals, the other rows of the file contain one SNP per line and the genotype for each individual:
rs# alleles chrom pos strand assembly# center protLSID assayLSID panelLSID QCcode NA06984 NA06985 NA06986 NA06989 NA06991 NA06993 NA06994 NA06995 NA06997 NA07000 NA07014 NA07019 NA07022 NA07029 NA07031 NA07034 NA07037 NA07045 NA07048 NA07051 NA07055 NA07056 NA07345 NA07346 NA07347 NA07348 NA07349 NA07357 NA07435 NA10830 NA10831 NA10835 NA10836 NA10837 NA10838 NA10839 NA10840 NA10843 NA10845 NA10846 NA10847 NA10850 NA10851 NA10852 NA10853 NA10854 NA10855 NA10856 NA10857 NA10859 NA10860 NA10861 NA10863 NA10864 NA10865 NA11829 NA11830 NA11831 NA11832 NA11839 NA11840 NA11843 NA11881 NA11882 NA11891 NA11892 NA11893 NA11894 NA11917 NA11918 NA11919 NA11920 NA11930 NA11931 NA11992 NA11993 NA11994 NA11995 NA12003 NA12004 NA12005 NA12006 NA12043 NA12044 NA12045 NA12056 NA12057 NA12144 NA12145 NA12146 NA12154 NA12155 NA12156 NA12234 NA12236 NA12239 NA12248 NA12249 NA12264 NA12272 NA12273 NA12275 NA12282 NA12283 NA12286 NA12287 NA12335 NA12336 NA12340 NA12341 NA12342 NA12343 NA12344 NA12347 NA12348 NA12375 NA12376 NA12383 NA12386 NA12399 NA12400 NA12413 NA12489 NA12546 NA12707 NA12708 NA12716 NA12717 NA12718 NA12739 NA12740 NA12748 NA12749 NA12750 NA12751 NA12752 NA12753 NA12760 NA12761 NA12762 NA12763 NA12766 NA12767 NA12775 NA12776 NA12777 NA12778 NA12801 NA12802 NA12812 NA12813 NA12814 NA12815 NA12817 NA12818 NA12827 NA12828 NA12829 NA12830 NA12832 NA12842 NA12843 NA12864 NA12865 NA12872 NA12873 NA12874 NA12875 NA12877 NA12878 NA12889 NA12890 NA12891 NA12892
rs6423165 A/G chrY 109805 + ncbi_b36 broad urn:LSID:affymetrix.hapmap.org:Protocol:GenomeWideSNP_6.0:3 urn:LSID:broad.hapmap.org:Assay:SNP_A-8572888:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ GG AG GG AG AG AG AA AA GG AA AG NN AG AA AG NN AA AG NN GG GG GG GG AA GG GG AG AG AA AG GG AG AG AG AA AG AG AG GG GG AA AA NN GG AG NN AA AG NN GG NN GG GG GG AA AG AA AA AG GG AG AG GG AG AA AG GG AG GG GG AG GG AG GG GG AG GG GG AA NN AG GG AG GG GG GG AG GG GG AA AA GG GG GG NN AG AG AG GG AG AG AA AA GG AA GG AA AG AA AG GG AG GG GG AG AG GG GG AG AG AG NN AG AG AG AG AG NN AA AA AG AG AG AG AG AA AA AG AG AA AG GG AG GG GG AG AG AG AG AG AG AA AG AG AG AG AA AA AG GG GG GG GG AG AG GG GG AG AA GG AG AG GG AG
rs6608381 C/T chrY 109922 + ncbi_b36 broad urn:LSID:affymetrix.hapmap.org:Protocol:GenomeWideSNP_6.0:3 urn:LSID:broad.hapmap.org:Assay:SNP_A-8528859:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ CC CT CC CT CT CT TT TT CC TT CT NN CT TT CT NN TT CC NN CC CC CC CC TT CC CC CT CT TT CT CC CT CT CT TT CT CT CT CC CC TT TT NN CC CT NN TT CT NN CC NN CC CC CC TT CT TT TT CT CC CT CT CC CT TT CT CC CT CC CC CT CC CT CC CC CT CC CC TT NN CT CC CT CC CC CC CT CC CC TT TT CC CC CC NN CT CT CT CC CT CT TT TT CC TT CC TT CT TT CT CC CT CC CC CT CT CC CC CT CT CT NN CT CT CT CT CT NN TT TT CT CT CT CT CT TT TT CT CT TT CT CC CT CC CC CT CT CC CT CT CC TT CT CT CT CT TT TT CT CC CC CC CC CT CT CC CC CT TT CC CT CT CC CT
rs6644972 A/G chrY 118624 + ncbi_b36 sanger urn:LSID:illumina.hapmap.org:Protocol:Human_1M_BeadChip:3 urn:LSID:sanger.hapmap.org:Assay:H1Mrs6644972:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ GG GG AG GG AG AG AG GG AG GG AG NN GG GG AG NN GG GG NN AG AG GG GG GG GG GG GG GG GG AG AA GG GG GG AG GG GG GG GG GG GG GG NN GG AG GG GG GG NN AG NN AG AG AA GG GG GG GG GG GG GG GG GG AG GG GG AA AG GG AG GG GG GG GG GG GG GG AG AG NN AG GG AG GG GG GG AG GG GG GG GG AA AA AG NN GG GG GG GG GG GG GG AG GG GG GG GG AG GG GG AG GG GG GG GG GG GG GG GG GG GG AG GG GG GG GG GG NN GG GG AG GG GG AG GG AG GG AG AA GG GG GG GG AG GG GG GG AG GG AG GG GG GG AG GG GG AG GG AG GG GG GG GG GG GG AG GG GG GG GG GG GG GG GG
rs6655397 A/G chrY 141935 + ncbi_b36 sanger urn:LSID:illumina.hapmap.org:Protocol:Human_1M_BeadChip:3 urn:LSID:sanger.hapmap.org:Assay:H1Mrs6655397:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ GG GG GG GG AG AG GG GG GG GG GG NN GG GG NN NN GG GG NN GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG NN AG GG GG GG GG NN AG NN GG GG GG NN AG GG GG GG AG GG GG AG GG GG GG GG AG GG AG GG GG GG GG AG AG GG GG GG NN AG GG GG GG GG GG GG GG GG GG GG GG GG GG NN GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG NN GG GG AG GG AG NN GG GG GG GG GG GG GG GG AG GG GG AG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG GG AG GG GG GG GG GG GG
rs6603172 C/T chrY 142439 + ncbi_b36 broad urn:LSID:affymetrix.hapmap.org:Protocol:GenomeWideSNP_6.0:3 urn:LSID:broad.hapmap.org:Assay:SNP_A-8527413:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ TT TT TT TT CT CT TT TT TT TT TT NN TT TT TT NN TT TT NN TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT NN CT TT TT TT TT NN CT NN TT TT TT TT CT TT TT TT CT TT TT CT TT TT TT TT CT TT CT TT TT TT TT CT CT TT TT TT NN CT TT TT TT TT TT TT TT TT TT TT TT TT TT NN TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT CT TT CT NN TT TT TT TT TT TT TT TT CT TT TT CT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT CT TT TT TT TT TT TT
rs6644970 A/G chrY 142664 + ncbi_b36 broad urn:LSID:affymetrix.hapmap.org:Protocol:GenomeWideSNP_6.0:3 urn:LSID:broad.hapmap.org:Assay:SNP_A-4207883:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ AG AA AA AG AG GG AG AG AA AA GG NN AG AA GG NN AG AG NN AG AG AG GG AG AG AG AA AG GG GG AA AG AG AG AG AG AA AA AG GG AG GG NN GG AG NN AG AG NN AG NN AA AG AA AA AG GG GG AG AG AG AG AG AG AA AG AA AG GG GG AG AG AA GG AG GG AA AG AG NN AG AG GG AG GG GG AG GG GG AA AG AA AG GG NN GG AA GG AG AG GG GG AA GG AA AG AG GG GG AG AG GG AG AG AG GG AG AG AG AG AG NN AG GG AG AG GG NN AA AG AA AA GG AG AG AG GG AG AG AG GG AG AA AG GG AA AA AG GG AA AG AG GG AA GG AA AA AG AG AG AG AG GG AG GG GG AG AG AG AA AA GG AG AG
rs11019 A/G chrY 159978 + ncbi_B36 affymetrix urn:LSID:affymetrix.hapmap.org:Protocol:GenomeWideSNP_6.0:2 urn:LSID:affymetrix.hapmap.org:Assay:SNP_A-4207841:2 urn:lsid:dcc.hapmap.org:Panel:CEPH-30-trios:1 QC+ NN AG NN NN AG GG AA NN NN AG NN AG GG AA NN AG NN NN GG NN GG AG GG NN NN GG NN GG NN GG AG GG NN NN GG AG NN NN NN AG AG NN NN NN NN AA AG GG GG AA AG GG GG NN NN GG GG AG GG AG AG NN AA AA NN NN NN NN NN NN NN NN NN NN GG AG GG AG AG GG AG AG GG GG NN GG GG AG AG GG GG AG AG AG AG AG GG GG 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 AG NN AA GG NN NN GG NN NN GG GG AG AA AG GG AA AA NN NN NN NN NN NN AG GG GG AG GG AG NN NN NN NN NN NN NN NN NN GG GG AG AG GG GG NN GG NN NN GG AG
rs1132353 C/T chrY 160116 + ncbi_b36 sanger urn:LSID:illumina.hapmap.org:Protocol:Human_1M_BeadChip:3 urn:LSID:sanger.hapmap.org:Assay:H1Mrs1132353:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ CC CT TT CT CT TT NN CT TT CT CC NN TT CC CT NN TT CT NN CT TT CT CT TT TT CC TT CT CC TT CT TT TT TT TT CT TT TT TT CT CT TT NN TT CT NN CT TT NN CC NN TT TT CC NN TT TT CT TT CT CT CC CC CC CT CT CC CT CC CT CT TT TT TT TT CT TT CT CT NN CT CT TT TT TT TT TT CC CT TT TT CT CT CT NN CT TT TT TT TT TT TT TT TT TT TT CT CT CC TT CT TT CT CT CT CT CT CC CC CC CC TT CC CT CT CT CC NN TT CT TT TT CT TT CT CT CC CT TT CC CC TT TT TT CT TT TT CT TT TT CT TT CT TT TT CT TT CT TT CT CT CT TT TT CT CT TT TT CT TT TT CC TT CT
rs35047434 C/T chrY 169118 + ncbi_b36 sanger urn:LSID:illumina.hapmap.org:Protocol:Human_1M_BeadChip:3 urn:LSID:sanger.hapmap.org:Assay:H1Mrs35047434:3 urn:lsid:dcc.hapmap.org:Panel:CEPH-60-trios:3 QC+ TT TT TT TT CT CT TT TT TT TT TT NN TT TT TT NN TT TT NN CT TT TT TT TT TT TT TT TT TT TT TT TT TT CT TT TT TT CT TT TT CT TT NN CT TT TT TT TT NN TT NN TT TT TT TT TT TT TT TT TT CT TT TT TT TT TT TT CT TT TT TT CT TT CT CT TT TT TT TT NN TT TT TT TT TT TT TT TT TT CT TT TT TT TT NN TT CT TT TT TT CT CT TT TT TT CT CT TT TT CT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT NN TT CT TT CT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT TT CT TT TT TT TT CT TT TT TT CT TT TT TT TT CT TT TT TT CC TT TT TT TT TT TT
(...)

Opening the SQLIte Database

if(sqlite3_open(env->fileout,&(env->connection))!=SQLITE_OK)
{
fprintf(stderr,"Cannot open sqlite file %s.\n",env->fileout);
exit(EXIT_FAILURE);
}

Closing the SQLIte Database

sqlite3_close(env->connection);

Creating a table

if(sqlite3_exec(env->connection,
"create table markers(id INTEGER PRIMARY KEY,name varchar(20) unique,chrom varchar(10),position integer)",
NULL,NULL,&error
)!=SQLITE_OK)
{
fprintf(stderr,"Cannot create table markers: %s\n",error);
exit(EXIT_FAILURE);
}

Filling a prepared statement and inserting some data

if(sqlite3_prepare(env->connection,
"insert into markers(name,chrom,position) values(?,?,?)",
-1,&pstmt_insert_marker,NULL
)!=SQLITE_OK)
{
fprintf(stderr,"Cannot compile insert into markers.\n");
exit(EXIT_FAILURE);
}
(...)

if(sqlite3_bind_text(
pstmt_insert_marker,1,
name,
strlen(name),
NULL)!=SQLITE_OK)
{
fprintf(stderr,"Cannot bind markers's name.\n");
exit(EXIT_FAILURE);
}
(
if(sqlite3_bind_text(
pstmt_insert_marker,2,
chrom,
strlen(chrom),
NULL)!=SQLITE_OK)
{
fprintf(stderr,"Cannot bind markers's chrom.\n");
exit(EXIT_FAILURE);
}
if( sqlite3_bind_int(
pstmt_insert_marker,3,
genomic_position
)!=SQLITE_OK)
{
fprintf(stderr,"Cannot bind markers's position.\n");
exit(EXIT_FAILURE);
}

if (sqlite3_step(pstmt_insert_marker) != SQLITE_DONE) {
fprintf(stderr,"Could not step insert marker.\n");
exit(EXIT_FAILURE);
}

marker_id = sqlite3_last_insert_rowid(env->connection);

sqlite3_reset(pstmt_insert_marker);

Creating a custom function for sqlite

Here I implemented a custom function named homozygous it returns true if a genotype is a string of two identical characters (but not 'N').
void is_homozygous(
sqlite3_context* context, /* context */
int argc, /* number of arguments */
sqlite3_value** argv /* arguments */
)
{
/* one text arg */
if(argc==1 &&
sqlite3_value_type(argv[0])==SQLITE_TEXT)
{
const char *alleles=(const char*)sqlite3_value_text(argv[0]);
if( alleles!=NULL &&
strlen(alleles)==2 &&
alleles[0]==alleles[1] &&
alleles[0]!='N')
{
/* set result to TRUE */
sqlite3_result_int(context,1);
return;
}
}
/* set result to FALSE */
sqlite3_result_int(context,0);
}

Telling sqlite about the new custom function

sqlite3_create_function(env->connection,
"homozygous",/* function name */
1,/* number of args */
SQLITE_UTF8,/* encoding */
NULL,
is_homozygous,/* the implementation callback */
NULL,
NULL
);

Implementing a callback for a simple 'SELECT' query

The following callback prints the results to stdout.
static int callback_select_homozygous(
void* notUsed,
int argc,/* number of args */
char** argv,/* arguments as string */
char** columns /* labels for the columns */
)
{
int i;
/* prints all the columns to stdout*/
for(i=0; i<argc; i++)
{
printf( "\"%s\": \"%s\"\n",
columns[i],
argv[i]!=NULL? argv[i] : "NULL"
);
}
printf("\n");
return 0;
}

Executing a SELECT query

Here, we select all the homoygous genotypes:
if(sqlite3_exec(
env->connection, /* An open database */
"select individuals.name,"
"markers.name,"
"markers.chrom,"
"markers.position,"
"genotypes.alleles "
"from "
"genotypes,"
"markers,"
"individuals "
"where "
"markers.id=genotypes.marker_id and "
"individuals.id=genotypes.individual_id and "
"homozygous(genotypes.alleles)=1"
, /* SQL to be evaluated */
callback_select_homozygous, /* Callback function */
NULL, /* 1st argument to callback */
&errormsg /* Error msg written here */
)!=SQLITE_OK)
{
fprintf(stderr,"Cannot select.\n");
exit(EXIT_FAILURE);
}

The full source code



Excuting the C program

[LOG]inserted marker ID.1
[LOG]inserted marker ID.2
[LOG]inserted marker ID.3
[LOG]inserted marker ID.4
[LOG]inserted marker ID.5
(... WAIT !!! )
"name": "NA12878"
"name": "rs6603251"
"chrom": "chrY"
"position": "240580"
"alleles": "CC"

"name": "NA12890"
"name": "rs6603251"
"chrom": "chrY"
"position": "240580"
"alleles": "CC"

"name": "NA12892"
"name": "rs6603251"
"chrom": "chrY"
"position": "240580"
"alleles": "CC"


That's it

Pierre

04 July 2007

Systems-Biology using GoogleGears: my notebook


Google gears is an open source browser extension that enables web applications to provide offline functionality. The data are stored locally in a fully-searchable relational database using the sqlite engine.


My Biological Network is a tool I created as a test to play with Google gears: it is used to build a network of protein-protein interactions. It uses Google Gears to record your entries on the local disk, so Gears needs to be installed on your computer. Programming with gears with JAVASCRIPT is really cool as you don't have to implement the storage of the data on the server side and you're using some standard SQL statements to handle the data.




Screenshots


My Biological Network


Tutorial


Open the tab Organism (fig. 4): add one or more organism. (Homo Sapiens already inserted by default)

Open the tab Protein (fig. 1): add one or more protein.

Open the tab Paper (fig. 3): add one or more article that will be used as an evidence for an interaction.

Open the tab Technology (fig. 2): add one or more technology that was used to characterize an interaction.

Open the tab Component: add one or more cellular component using Gene Ontology (GO:0005575 \"cellular component\" was inserted by default)

Open the tab Interaction (fig. 5):


  • Name and describe this interaction

  • Select one or more protein and/or one or more previously defined proteic complex. You Cannot describe self interactions with this tool.

  • (optional) choose one or more paper/technology/component...



Open the RDF table (fig. 6): I choose to display the content of the database using RDF. Such format can then be validated and visualized using the W3C RDF validator, or transformed using XSLT, etc.... I also used the life science identifier (LSID) as an URI for my resources.


On my computer, the database is stored in /env/islande/home/lindenb/.mozilla/firefox/<profile-id>/Google Gears for Firefox/islande/<host>/mynetwork#database. The database can be manualy accessed using sqlite3:

sqlite3 mynetwork#database
SQLite version 3.4.0
Enter '.help' for instructions
sqlite> .tables
component interactionhash paper technology
interaction organism prote
sqlite> .schema organism
CREATE TABLE organism(id integer primary key ,name varchar(50) not null unique);
sqlite> select * from organism;
9606|Homo Sapiens
sqlite>


Internals


We the page is loaded, we check that gears was installed



if (!window.google || !google.gears) {
debug("NOTE: You must install Google Gears first.")

We then create the database if does not exist. The file is created in firefox in ${HOME}/.mozilla/firefox/<profile-id>/Google Gears for Firefox/<server>/mynetwork#database


connection = google.gears.factory.create("beta.database","1.0");

I create the tables just by invoking some standards SQL 'CREATE TABLE' statements. I also insert some default values (e.g. human organism)



connection.execute("create table if not exists organism(id integer primary key ,name varchar(50) not null unique)");
connection.execute("insert or ignore into organism(id,name) values(9606,\"Homo Sapiens\")");
connection.execute("create table if not exists protein(id integerprimary key autoincrement,name varchar(50) not null,taxId int not null,acn varchar(50) not null unique)");
connection.execute("create table if not exists paper(pmid integerprimary key ,title varchar(255) not null,citation varchar(255) not null,firstAuthor varchar(50) not null)");
connection.execute("create table if not exists component(id integer primary key autoincrement,go varchar(50) not null unique, name varchar(50) not null unique)");

connection.execute("insert or ignore into component(go,name) values(\"GO:0005575\",\"cellular component\")");
connection.execute("insert or ignore into component(go,name) values(\"GO:0008372\",\"cellular component unknown\")");

connection.execute("create table if not exists technology(id integer primary key autoincrement,name varchar(50) not null unique, description varchar(255) not null)");

connection.execute("insert or ignore into technology(name,description) values(\"Y2H\",\"Yeast Two Hybrid System\")");
connection.execute("insert or ignore into technology(name,description) values(\"CoIP\",\"Co-Immuno Precipitation\")");


connection.execute("create table if not exists interaction(id integer primary key autoincrement, name varchar(50) not null unique,description varchar(255) not null)");
connection.execute("create table if not exists interactionhash(id integer primary key autoincrement,LINK_interaction int ,type varchar(20) not null,child int not null)");

When a data is about to be inserted we check all the fields and we insert them using SQL: INSERT INTO


var id= getById("organism-input-id");
if(!isInteger(id.value))
{
debug("TaxId not a Number");
return;
}
var name=getById("organism-input-name");
if(trim(name.value).length==0)
{
debug("Taxon Name empty");
return;
}

try
{
connection.execute("insert into organism(id,name) values("+sqlescape(trim(id.value))+","+sqlquote(trim(name.value))+")");
id.value="";
name.value="";
}
catch(err)
{
debug(err.message);
return;
}

a simple SELECT is used to retrieve the data and insert them in a HTML table



var rs= connection.execute("select id,name from organism order by name");
while (rs.isValidRow())
{
var tr= ce("tr");
table.appendChild(tr);
var td= ce("td");
tr.appendChild(td);
td.appendChild(ct(rs.field(0)));

td= ce("td");
tr.appendChild(td);
var a= ce("a");
a.setAttribute("title","Open in NCBI");
a.setAttribute("target","tax"+rs.field(0));
a.setAttribute("href","http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id="+rs.field(0));
td.appendChild(a);
a.appendChild(ct(rs.field(1)));
rs.next();
}
rs.close();



That's it !

Pierre

updated 2010-08-12: source code

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<script type="text/javascript" src="gears_init.js"></script>
<script type="text/javascript" src="network.js"></script>
<link rel="stylesheet" type="text/css" href="./network.css" />
<title>My Biological Network</title>
</head>
<body onload="init()">
<h1>My Biological Network</h1>
<p>Pierre Lindenbaum PhD <a href="mailto:plindenbaum@yahoo.fr">plindenbaum@yahoo.fr</a><br/><a href="http://plindenbaum.blogspot.com">http://plindenbaum.blogspot.com</a><br/><address>Bioinformatics department<br/><a href="http://www.integragen.com">Integragen S.A.</a><br/>Evry, France</address></p>
<p/>
<div>
<button onclick="javascript:showCard('home-pane');">Home</button>
<button onclick="showOrganismPane()">Organisms</button>
<button onclick="showProteinPane()">Proteins</button>
<button onclick="showPaperPane()">Papers</button>
<button onclick="showTechnologyPane()">Technology</button>
<button onclick="showComponentPane()">Component</button>
<button onclick="showInteractionPane()">Interactions</button>
<button onclick="showRDFPane()">RDF</button>
</div>
<div style="color:red;" id="stderr"></div>
<p/>

<!-- ====================================== ORGANISM ====================================== -->
<div style="display:none;" id="organism-pane">
<table>
<caption>Add an Organism</caption>
<tr><th>NCBI Taxon ID <i>(e.g. 10912)</i></th><td><input id="organism-input-id" length="10"/></td></tr>
<tr><th>NCBI Taxon Name <i>(e.g. Rotavirus)</i></th><td><input id="organism-input-name" length="10"/></td></tr>
<tr><th/><td><button onclick="addOrganism()">Add</button></td></tr>
</table>

<hr/>

<table width="80%">
<caption>All Organisms</caption>
<thead>
<tr><th>Taxon ID</th><th>Taxon Name</th></tr></tr>
</thead>
<tbody id="organism-table">
</tbody>
</table>

</div>

<!-- ====================================== COMPONENT ====================================== -->
<div style="display:none;" id="component-pane">
<table>
<caption>Add a Component</caption>
<tr><th>Name</th><td><input id="component-input-name" length="10"/></td></tr>
<tr><th>GO</th><td><input id="component-input-go" length="10"/></td></tr>
<tr><th/><td><button onclick="addComponent()">Add</button></td></tr>
</table>

<hr/>

<table width="80%">
<caption>All Components</caption>
<thead>
<tr><th>Name</th><th>GO</th></tr>
</thead>
<tbody id="component-table">
</tbody>
</table>

</div>

<!-- ====================================== TECHNOLOGY ====================================== -->
<div style="display:none;" id="technology-pane">
<table>
<caption>Add a Technology</caption>
<tr><th>Name</th><td><input id="technology-input-name" length="50"/></td></tr>
<tr><th>Description</th><td><input id="technology-input-desc" length="50"/></td></tr>
<tr><th/><td><button onclick="addTechnology()">Add</button></td></tr>
</table>

<hr/>

<table width="80%">
<caption>All Technologies</caption>
<thead>
<tr><th>Name</th><th>Description</th></tr></tr>
</thead>
<tbody id="technology-table">
</tbody>
</table>

</div>


<!-- ====================================== PROTEIN ====================================== -->

<div style="display:none;" id="protein-pane">
<table>
<caption>Add a Protein</caption>
<tr><th>Uniprot accession number <i>(e.g. Q3T8J2)</i></th><td><input id="protein-input-acn" length="10"/></td></tr>
<tr><th>Uniprot Name <i>(e.g. Replicase polyprotein 1ab)</i></th><td><input id="protein-input-name" length="10"/></td></tr>
<tr><th>Organism</th><td><select id="protein-input-taxon" length="10"><option>A</option></select></td></tr>
<tr><th/><td><button onclick="addProtein()">Add</button></td></tr>
</table>

<hr/>

<table width="80%">
<caption>All Proteins</caption>
<thead>
<tr><th>Primary accession</th><th>Name</th><th>Taxon</th></tr></tr>
</thead>
<tbody id="protein-table">
</tbody>
</table>

</div>

<!-- ====================================== PAPER ====================================== -->
<div style="display:none;" id="paper-pane">
<table>
<caption>Add a Paper</caption>
<tr><th>PMID</th><td><input id="paper-input-pmid" length="10"/></td></tr>
<tr><th>Title</th><td><input id="paper-input-title" length="50"/></td></tr>
<tr><th>Citation</th><td><input id="paper-input-citation" length="50"/></td></tr>
<tr><th>First Author</th><td><input id="paper-input-author" length="50"/></td></tr>
<tr><th/><td><button onclick="addPaper()">Add</button></td></tr>
</table>

<hr/>

<table width="80%">
<caption>All Papers</caption>
<thead>
<tr><th>PMID</th><th>Citation</th><th>First Author</th><th>Title</th></tr></tr>
</thead>
<tbody id="paper-table">
</tbody>
</table>

</div>

<!-- ====================================== INTERACTION ====================================== -->


<div style="display:none;" id="interaction-pane">

<table>
<caption>Add an Interaction</caption>
<tr><th>Name</th><td colspan="4"><input id="interaction-input-name" length="50"/></td></tr>
<tr><th>Description</th><td colspan="4"><input id="interaction-input-desc" length="50"/></td></tr>
<tr>
<th>Protein</th>
<th>Interactors</th>
<th>Methods</th>
<th>Evidences</th>
<th>Components</th></tr>
<tr>
<td><select id="interactors-input-proteins" size="5" multiple="true"/></td>
<td><select id="interactors-input-interactors" size="5" multiple="true"></td>
<td><select id="interactors-input-technologies" size="5" multiple="true"></td>
<td><select id="interactors-input-evidences" size="5" multiple="true"></td>
<td><select id="interactors-input-components" size="5" multiple="true"></td>
</tr>
<tr><th colspan="4"/><td><button onclick="addInteraction()">Add</button></td></tr>
</table>

<hr/>

<table width="80%">
<caption>All Interactions</caption>
<thead>
<tr><th>Name</th><th>Description</th></tr></tr>
</thead>
<tbody id="interaction-table">
</tbody>
</table>

</div>

<!-- ====================================== RDF ====================================== -->
<div style="display:none;" id="rdf-pane">
<h2>RDF Pane</h2>
<textarea wrap="off" id="rdf-area" rows="20" cols="80"></textarea>

</div>

<!-- ====================================== HOME ====================================== -->
<div style="display:none;" id="home-pane">
<h3>About My Biological Network</h3>
<p><a href="http://gears.google.com/">Google gears</a> is an open source browser extension that enables web applications to provide offline functionality. The data are stored locally in a fully-searchable relational database using the <a href="http://www.sqlite.org/">sqlite engine</a>.</p>
<p><b>My Biological Network</b> is a tool I created as a test to play with Google gears: it is used to build a network of protein-protein interactions. It uses Google Gears to record your entries on the <u>local disk</u>, so Gears needs to be installed on your computer. </p>

<p>
Open the tab <b>Organism</b>: add one or more organism. (Homo Sapiens already inserted by default)<br/>
Open the tab <b>Protein</b>: add one or more protein.<br/>
Open the tab <b>Paper</b>: add one or more article that will be used as an evidence for an interaction.<br/>
Open the tab <b>Technology</b>: add one or more technology that was used to characterize an interaction.<br/>
Open the tab <b>Component</b>: add one or more cellular component using Gene Ontology (GO:0005575 \"cellular component\" was inserted by default)<br/>
Open the tab <b>Interaction</b>:<ul>
<li>Name and describe this interaction</li>
<li>Select one or more protein and/or one or more previously defined proteic complex. You <i>Cannot</i> describe self interactions with this tool.<li>
<li>(optional) choose one or more paper/technology/component...</li>
</ul><br/>
Open the <b>RDF table</b>: I choose to display the content of the database using <a href="http://www.w3.org/RDF/">RDF</a>. Such format can then be validated and visualized using the <a href="http://www.w3.org/RDF/Validator/">W3C RDF validator</a>, or transformed using <a href="http://www.w3.org/TR/xslt">XSLT</a>, etc.... I also used the <a href="http://lsid.sourceforge.net/">life science identifier (LSID)</a> as an URI for my resources.<br/>

</p>

<p>On my computer, the database is stored in <code>$HOME/.mozilla/firefox/&lt;profile-id&gt;/Google Gears for Firefox/&lt;host&gt;/mynetwork#database</code>. The database can be manualy accessed using <a href="http://www.sqlite.org/">sqlite3</a>:<pre style='color:black;border:1pt solid;background:lightgray;'>sqlite3 mynetwork#database
SQLite version 3.4.0
Enter &apos;.help&apos; for instructions
sqlite&gt; .tables
component interactionhash paper technology
interaction organism prote
sqlite&gt; .schema organism
CREATE TABLE organism(id integer primary key ,name varchar(50) not null unique);
sqlite&gt; select * from organism;
9606|Homo Sapiens
sqlite&gt;</pre>

</p>

</div>


<!-- google analytics -->

<script src="http://www.google-analytics.com/urchin.js"
type="text/javascript">
</script>
<script type="text/javascript">
_uacct = "XXXXXX";
urchinTracker();
</script>

<!-- google analytics -->


</body>
</html>