YOKOFAKUN

20 July 2011

Drawing a SVG timeline with the http://data.bnf.fr data

The National French Library has recently started to release
its data as RDF/XML. Here, I've played with the biographies of the famous French writers to create a simple timeline.

Unfortunately, this timeline is too large to be displayed here. :-)

However, here is the java code I used to generate this map as a SVG document:

See also: Freebase and the History of Sciences

That's it,

Pierre

18 July 2011

An interactive 'dialog' generating a SAM Flag

The following shell script uses the dialog utility to generate a numeric Sam flag.

$ sh selectsam.sh 
      ┌───────────────────────SAM FLAGS──────────────────────────┐
      │ SELECT FLAGS                                             │  
      │ ┌──────────────────────────────────────────────────────┐ │  
      │ │ [X] 1     Read Paired                                │ │  
      │ │ [ ] 2     Read mapped in proper pair                 │ │  
      │ │ [X] 4     Read unmapped                              │ │  
      │ │ [ ] 8     Mate unmapped                              │ │  
      │ │ [ ] 16    Read reverse strand                        │ │  
      │ │ [X] 32    Mate reverse strand                        │ │  
      │ │ [ ] 64    First in pair                              │ │  
      │ │ [X] 128   Second in pair                             │ │  
      │ │ [ ] 256   Not primary alignment                      │ │  
      │ │ [ ] 512   Read fails platform/vendor quality checks  │ │  
      │ │ [ ] 1024  Read is PCR or optical duplicate           │ │  
      │ └─────↓(+)─────────────────────────────────────────────┘ │  
      │                                                          │  
      │                                                          │  
      ├──────────────────────────────────────────────────────────┤  
      │                         <  OK  >                         │  
      └──────────────────────────────────────────────────────────┘  
        
$ 165

Source available at: https://github.com/lindenb/samtools-utilities/blob/master/script/selectflag.sh

That's it,

Pierre

15 July 2011

I can't save you, but I can help you: reverse-complementing a DNA in OpenOffice

People in my lab are using "Microsoft Word" to edit and annotate their sequences. Just like Neil, I can't save them but I can help them: here is an OpenOffice / LibreOffice macro reverse-complementing the current selected text.

open libreoffice/openoffice
Open menu "Tools / Macros / Organize macros / Basic...
Select "My Macros/Standard" and click on "New"
Copy+paste my macro here
Click on the top button "compile"
Close the dialog
Menu "Tools/Customize" , "Keyboard" tab, "Toolbar content" . Button "Add... " / Category "Libre Office Macros" . Select "My Macro/Standard/Modulexxx"/ReverseComplementSelection". Click on "Modify" to assign an icon to this new button
There is now a new button in the toolbar. Clicking on that button reverse-complements the selected text.

That's it,

Pierre

PS: Hum ? what did you say, they use Microsoft Word ? Not OpenOffice ? pfff....

14 July 2011

A text alignment viewer using the samtools API

I wrote a modified version of samtools tview. The original code uses the curses API to display an interactive sequence viewer.

I wanted to generate this kind of screen for a large number of positions and I wanted to be able to redirect the output to a unix pipeline so I changed the original code in order to only handle an extensible 2-dimensional array of characters. The code I wrote is available on github at:

https://github.com/lindenb/samtools-utilities/blob/master/src/ttview.c

Compilation

cd samtools/
make #compile samtools
gcc -o bamttview -g -Wall -O2 -DSTANDALONE_VERSION  -I. -Lbcftools  bam_ttview.c  bam2bcf.o   errmod.o  bam_color.o libbam.a -lbcf  -lm -lz

Example

print a position

$ ./bamttview  -g "ref:5" examples/toy.bam examples/toy.fa |\
cat -n
     1       11              21        31         41        51        61          
     2 TGTTAGATAA****GATA**GCTGTGCTAGTAGGCAG*TCAGCGCCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
     3   ........    ....  ......K.K......K. ..........                                
     4   ........AGAG....***...      ,,,,,    ,,,,,,,,,                                
     5     ......GG**....AA                                                            
     6     ..C...**** ...**...>>>>>>>>>>>>>>T.....

print a list of positions:

$ cat positions.txt 
ref2:10-100
ref:10-15
ref2:11

$ ./bamttview  -f positions.txt examples/toy.bam examples/toy.fa |\
cat -n | head -n 25
     1 
     2 
     3 > ref2:10-100
     4 
     5  11            21        31        41        51        61        71             
     6 aaaac****aattaagtctacagagcaactaNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
     7 ....Y    ..W...................                                                 
     8 .....****..A...                                                                 
     9 .....****..A...T.                                                               
    10 .....AAAT.............                                                          
    11 C...T****....................                                                   
    12   ..T****.....................                                                  
    13     T****......................                                                 
    14 
    15 
    16 
    17 > ref:10-15
    18 
    19  11              21        31         41        51        61        71          
    20 GATAA****GATA**GCTGTGCTAGTAGGCAG*TCAGCGCCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
    21 .....    ....  ......K.K......K. ..........                                     
    22 .....AGAG....***...      ,,,,,    ,,,,,,,,,                                     
    23 .....GG**....AA                                                                 
    24 .C...**** ...**...>>>>>>>>>>>>>>T.....                                          
    25 

./bamttview  -f positions.txt -d examples/toy.bam examples/toy.fa |\
cat -n | head -n 25
     1 
     2 
     3 > ref2:10-100
     4 
     5  11            21        31        41        51        61        71             
     6 aaaac****aattaagtctacagagcaactaNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
     7 ....Y    ..W...................                                                 
     8 AAAAC****AAATAA                                                                 
     9 AAAAC****AAATAATT                                                               
    10 AAAACAAATAATTAAGTCTACA                                                          
    11 CAAAT****AATTAAGTCTACAGAGCAAC                                                   
    12   AAT****AATTAAGTCTACAGAGCAACT                                                  
    13     T****AATTAAGTCTACAGAGCAACTA                                                 
    14 
    15 
    16 
    17 > ref:10-15
    18 
    19  11              21        31         41        51        61        71          
    20 GATAA****GATA**GCTGTGCTAGTAGGCAG*TCAGCGCCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
    21 .....    ....  ......K.K......K. ..........                                     
    22 GATAAAGAGGATA***CTG      taggc    cagcgccat                                     
    23 GATAAGG**GATAAA                                                                 
    24 GCTAA**** ATA**GCT>>>>>>>>>>>>>>TTCAGC                                          
    25

That's it,

Pierre

09 July 2011

Storing SNPs in a HDF5 file: my notebook

I discovered the benefits of using twitter in 2008, the day Deepak Singh replied to one of my tweets related to the storage of a large number of genotypes.

Since that day, I've tried to use the HDF5 library, without any success (there's a large disheartening documentation/API on the HDF5 site and the API seems to be only used by a small number of geeks). Furthermore, HDF5 is a technology used by the IGV genome browser.

In that post, I'll describe a C program loading a set of SNPs defined by the following C structure:

typedef struct structSnp
 {
 char rsId[RS_LENGTH];//rs##
 char chrom[CHROM_LENGTH];//chromosome name
 int chromStart;//genomic start index
 int chromEnd;//genomic end index
 }Snp;

The SNPs will be read from stdin. Four columns are expected: rs-id, chrom, chromStart and chromEnd. Here is the command line I used to get a small input file:

curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp132.txt.gz" |\
 gunzip -c |\
 cut -d ' ' -f 2-5 |\
 head -n 10000 > sample.tsv

Source code

The C program is described below. I hope my comments will make the code readable.

The Makefile

Compile and run

gcc -Wall -o a.out -I  /path/to/hdf5/include -L /path/to/hdf5/lib dbsnp2hdf5.c  -lhdf5
./a.out < sample.tsv

Dump the data

At the en, the program creates a structured binary file containing our SNPs. The HDF5 utility h5dump can be used to display the data as text:

$ h5dump storage.h5

HDF5 "storage.h5" {
GROUP "/" {
   GROUP "variations" {
      DATASET "dbSNP" {
         DATATYPE  H5T_COMPOUND {
            H5T_STRING {
               STRSIZE 13;
               STRPAD H5T_STR_NULLTERM;
               CSET H5T_CSET_ASCII;
               CTYPE H5T_C_S1;
            } "rs";
            H5T_STRING {
               STRSIZE 7;
               STRPAD H5T_STR_NULLTERM;
               CSET H5T_CSET_ASCII;
               CTYPE H5T_C_S1;
            } "chrom";
            H5T_STD_I32LE "chromStart";
            H5T_STD_I32LE "chromEnd";
         }
         DATASPACE  SIMPLE { ( 10000 ) / ( H5S_UNLIMITED ) }
         DATA {
         (0): {
               "rs112750067",
               "chr1",
               10326,
               10327
            },
         (1): {
               "rs56289060",
               "chr1",
               10433,
               10433
            },
         (2): {
               "rs112155239",
               "chr1",
               10439,
               10440
            },
         (3): {
               "rs112766696",
               "chr1",
               10439,
               10440
            },
          (...)
               450425,
               450426
            }
         }
         ATTRIBUTE "rdfs:comment" {
            DATATYPE  H5T_STRING {
                  STRSIZE 15;
                  STRPAD H5T_STR_NULLTERM;
                  CSET H5T_CSET_ASCII;
                  CTYPE H5T_C_S1;
               }
            DATASPACE  SCALAR
            DATA {
            (0): "My list of SNPs"
            }
         }
      }
   }
}
}

That's it !
Pierre

04 July 2011

The Neo4j REST API. My notebook

Neo4j is a open-source graph engine implemented in Java. This post is my notebook for the Neo4J-server, a server combining a REST API and a webadmin application into a single stand-alone server.

Download & extract Neo4J

Download Neo4J from here

$ tar xfz neo4j-community-1.4.M04-unix.tar.gz
$ rm neo4j-community-1.4.M04-unix.tar.gz
$ cd neo4j-community-1.4.M04/

Start the Neo4J Server

#edit the config file if needed
$ nano conf/neo4j-server.properties
#start the server
$ ./bin/neo4j start
Starting Neo4j Server...
Waiting for Neo4j Server.....
7/4/11 7:07:13 PM org.neo4j.server.database.Database INFO: Using database at NEO4J/neo4j-community-1.4.M04/data/graph.db
7/4/11 7:07:13 PM org.neo4j.server.modules.DiscoveryModule INFO: Mounted discovery module at [/]
Adding JAXRS packages [org.neo4j.server.rest.discovery] at [/]
Adding JAXRS packages [org.neo4j.server.rest.web] at [/db/data]
Adding JAXRS packages [org.neo4j.server.webadmin.rest] at [/db/manage]
7/4/11 7:07:13 PM org.neo4j.server.modules.RESTApiModule INFO: Mounted REST API at [/db/data/]
7/4/11 7:07:13 PM org.neo4j.server.modules.ManagementApiModule INFO: Mounted management API at [/db/manage/]
7/4/11 7:07:13 PM org.neo4j.server.modules.WebAdminModule INFO: Mounted webadmin at [/webadmin]
7/4/11 7:07:13 PM org.neo4j.server.NeoServerWithEmbeddedWebServer INFO: Starting Neo Server on port [7474]
7/4/11 7:07:13 PM org.neo4j.server.web.Jetty6WebServer INFO: Mounting static content at [/webadmin] from [webadmin-html]
7/4/11 7:07:15 PM org.neo4j.server.NeoServerWithEmbeddedWebServer INFO: Server started on [http://okazaki:7474/]
running: PID:2816

Get the server root

$ curl -D - -H Accept:application/json "http://localhost:7474/db/data/"
HTTP/1.1 200 OK
Content-Length: 410
Content-Encoding: UTF-8
Content-Type: application/json
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

{
  "relationship_index" : "http://localhost:7474/db/data/index/relationship",
  "node" : "http://localhost:7474/db/data/node",
  "relationship_types" : "http://localhost:7474/db/data/relationship/types",
  "extensions_info" : "http://localhost:7474/db/data/ext",
  "node_index" : "http://localhost:7474/db/data/index/node",
  "reference_node" : "http://localhost:7474/db/data/node/0",
  "extensions" : {
  }

Create an empty node

$ curl -D - -H Accept:application/json -X POST http://localhost:7474/db/data/node
HTTP/1.1 201 Created
Content-Length: 968
Location: http://localhost:7474/db/data/node/2
Content-Encoding: UTF-8
Content-Type: application/json
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

{
  "outgoing_relationships" : "http://localhost:7474/db/data/node/2/relationships/out",
(...)

Set properties for this node

$ curl -D - -H Content-Type:application/json -X PUT \
   -d '{"name":"Charles Darwin","birthDate":"1809-02-12","deathDate":"1882-04-19","knownFor":["Voyage of the Beagle","On the Origin of Species evolution by natural selection"]}' \
   http://localhost:7474/db/data/node/2/properties
HTTP/1.1 204 No Content
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

View this node

$ curl -D - -H Accept:application/json  http://localhost:7474/db/data/node/2HTTP/1.1 200 OK
Content-Length: 1166
Content-Encoding: UTF-8
Content-Type: application/json
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

{
  "outgoing_relationships" : "http://localhost:7474/db/data/node/2/relationships/out",
  "data" : {
    "knownFor" : [ "Voyage of the Beagle", "On the Origin of Species evolution by natural selection" ],
    "name" : "Charles Darwin",
    "birthDate" : "1809-02-12",
    "deathDate" : "1882-04-19"
  },
  "traverse" : "http://localhost:7474/db/data/node/2/traverse/{returnType}",
  "all_typed_relationships" : "http://localhost:7474/db/data/node/2/relationships/all/{-list|&|types}",
  "property" : "http://localhost:7474/db/data/node/2/properties/{key}",
  "self" : "http://localhost:7474/db/data/node/2",
  "properties" : "http://localhost:7474/db/data/node/2/properties",
  "outgoing_typed_relationships" : "http://localhost:7474/db/data/node/2/relationships/out/{-list|&|types}",
  "incoming_relationships" : "http://localhost:7474/db/data/node/2/relationships/in",
  "extensions" : {
  },
  "create_relationship" : "http://localhost:7474/db/data/node/2/relationships",
  "all_relationships" : "http://localhost:7474/db/data/node/2/relationships/all",
  "incoming_typed_relationships" : "http://localhost:7474/db/data/node/2/relationships/in/{-list|&|types}"
}

Only show the properties for that node:

$ curl -H Accept:application/json  http://localhost:7474/db/data/node/2/properties
{
  "knownFor" : [ "Voyage of the Beagle", "On the Origin of Species evolution by natural selection" ],
  "name" : "Charles Darwin",
  "birthDate" : "1809-02-12",
  "deathDate" : "1882-04-19"
}

Create a new node with some properties

$ curl -D - -H Accept:application/json -H Content-Type:application/json -X POST -d '{"name":"Alfred Russel Wallace","birthDate":"1823-01-08","deathDate":"1913-11-07","knownFor":["natural selection","biogeography"]}' "http://localhost:7474/db/data/node"
HTTP/1.1 201 Created
Content-Length: 1127
Location: http://localhost:7474/db/data/node/3
Content-Encoding: UTF-8
Content-Type: application/json
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

{
  "outgoing_relationships" : "http://localhost:7474/db/data/node/3/relationships/out",
  "data" : {
    "knownFor" : [ "natural selection", "biogeography" ],
    "name" : "Alfred Russel Wallace",
    "birthDate" : "1823-01-08",
    "deathDate" : "1913-11-07"
  },
  "traverse" : "http://localhost:7474/db/data/node/3/traverse/{returnType}",
  (...)
}

Set one property

$ curl D - -H Accept:application/json -H Content-Type:application/json -X PUT \
     -d '"United Kingdom"' \
     "http://localhost:7474/db/data/node/3/properties/citizenship"
HTTP/1.1 204 No Content
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

$ curl -H Accept:application/json  http://localhost:7474/db/data/node/3/properties{
  "knownFor" : [ "natural selection", "biogeography" ],
  "name" : "Alfred Russel Wallace",
  "citizenship" : "United Kingdom",
  "birthDate" : "1823-01-08",
  "deathDate" : "1913-11-07"
}

Remove a node

$ curl -D - -H Accept:application/json -X POST http://localhost:7474/db/data/node
HTTP/1.1 201 Created
Content-Length: 968
Location: http://localhost:7474/db/data/node/4
Content-Encoding: UTF-8
Content-Type: application/json
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

$ curl -D - -X DELETE http://localhost:7474/db/data/node/4
HTTP/1.1 204 No Content
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

Create a relationship from Darwin to Wallace

$ curl -D - -H Accept:application/json -H Content-Type:application/json -X POST -d '{"type":"KNOWS","to":"http://localhost:7474/db/data/node/3","data":{"ref":"http://en.wikipedia.org/wiki/Charles_Darwin"}}' "http://localhost:7474/db/data/node/2/relationships"
HTTP/1.1 201 Created
Content-Length: 439
Location: http://localhost:7474/db/data/relationship/0
Content-Encoding: UTF-8
Content-Type: application/json
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

{
  "start" : "http://localhost:7474/db/data/node/2",
  "data" : {
    "ref" : "http://en.wikipedia.org/wiki/Charles_Darwin"
  },
  "self" : "http://localhost:7474/db/data/relationship/0",
  "property" : "http://localhost:7474/db/data/relationship/0/properties/{key}",
  "properties" : "http://localhost:7474/db/data/relationship/0/properties",
  "type" : "KNOWS",
  "extensions" : {
  },
  "end" : "http://localhost:7474/db/data/node/3"
}

#view properties for this relationship
$ curl -H Content-Type:application/json "http://localhost:7474/db/data/relationship/0/properties" 
{
  "ref" : "http://en.wikipedia.org/wiki/Charles_Darwin"
}

Add another property to the relationship

$ curl -D - -H Content-Type:application/json  -X PUT -d '"Darwin received a letter from Wallace asking if the book would examine human origins"' "http://localhost:7474/db/data/relationship/0/properties/comment"
HTTP/1.1 204 No Content
Access-Control-Allow-Origin: *
Server: Jetty(6.1.25)

$ curl -H Content-Type:application/json "http://localhost:7474/db/data/relationship/0/properties"
{
  "ref" : "http://en.wikipedia.org/wiki/Charles_Darwin",
  "comment" : "Darwin received a letter from Wallace asking if the book would examine human origins"
}

List Types of Relationship

$ curl -H Content-Type:application/json "http://localhost:7474/db/data/relationship/types"
["KNOWS"]

List Relationships

#from Darwin
$ curl -H Content-Type:application/json "http://localhost:7474/db/data/node/2/relationships/out/KNOWS"
[ {
  "start" : "http://localhost:7474/db/data/node/2",
  "data" : {
    "ref" : "http://en.wikipedia.org/wiki/Charles_Darwin",
    "comment" : "Darwin received a letter from Wallace asking if the book would examine human origins"
  },
  "self" : "http://localhost:7474/db/data/relationship/0",
  "property" : "http://localhost:7474/db/data/relationship/0/properties/{key}",
  "properties" : "http://localhost:7474/db/data/relationship/0/properties",
  "type" : "KNOWS",
  "extensions" : {
  },
  "end" : "http://localhost:7474/db/data/node/3"
} ]
#in to Darwin
$ curl -H Content-Type:application/json "http://localhost:7474/db/data/node/2/relationships/in/KNOWS"
[ ]
#out from wallace
$ curl -H Content-Type:application/json "http://localhost:7474/db/data/node/3/relationships/out/KNOWS"
[ ]
#all from/to wallace
$ curl -H Content-Type:application/json "http://localhost:7474/db/data/node/3/relationships/all/KNOWS"
[ {
  "start" : "http://localhost:7474/db/data/node/2",
  "data" : {
    "ref" : "http://en.wikipedia.org/wiki/Charles_Darwin",
    "comment" : "Darwin received a letter from Wallace asking if the book would examine human origins"
  },
  "self" : "http://localhost:7474/db/data/relationship/0",
  "property" : "http://localhost:7474/db/data/relationship/0/properties/{key}",
  "properties" : "http://localhost:7474/db/data/relationship/0/properties",
  "type" : "KNOWS",
  "extensions" : {
  },
  "end" : "http://localhost:7474/db/data/node/3"
} ]

Stop the Neo4J Server

$ ./bin/neo4j stop
Stopping Neo4j Server...
7/4/11 7:09:30 PM org.neo4j.server.NeoServerBootstrapper INFO: Neo4j Server shutdown initiated by kill signal
7/4/11 7:09:30 PM org.neo4j.server.NeoServerWithEmbeddedWebServer INFO: Successfully shutdown Neo Server on port [7474]
Waiting for Neo4j Server to exit...
Stopped Neo4j Server.

That's all for today, next time I'll dive into the indexes.