06 September 2011

Merging two XML files using the libxml streaming API.

I just wrote a program appending one or more XML files to another XML file using the libxml streaming API. This program named mergexml is available on github at: mergexml.

Usage

$ mergexml -h
mergexml: Pierre Lindenbaum PHD. 2011.
Compilation: Sep 6 2011 at 19:54:22.
Usage:
mergexml [options] -i database.xml xml1, xml2, ....xmln
Options:
-i (required)
-o (default:stdout)
--replace replace xml-in
-r ignore root of XML docs

Example

Say you have a RDF-based collection of quotes in a large file database.rdf:
<?xml version="1.0" encoding="UTF-8"?>
<rdf:RDF
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:foaf="http://xmlns.com/foaf/0.1/"
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:dcterms="http://purl.org/dc/terms/"
xmlns:q="http://quote.lindenb.org/ontology"
>
<!-- (...) -->
<foaf:Person rdf:about="http://en.wikipedia.org/wiki/Louis_Pasteur">
<foaf:name>Louis Pasteur</foaf:name>
</foaf:Person>

<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Louis_Pasteur"/>
<q:quote xml:lang="en">Science knows no country, because knowledge belongs to humanity, and is the torch which illuminates the world.</q:quote>
</q:Quote>

<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Louis_Pasteur"/>
<q:quote xml:lang="en">I am utterly convinced that Science and Peace will triumph over Ignorance and War, that nations will eventually unite not to destroy but to edify, and that the future will belong to those who have done the most for the sake of suffering humanity.</q:quote>
</q:Quote>

</rdf:RDF>
And you want to append another file 'chunck.rdf' at the end to 'database.rdf':
<?xml version="1.0" encoding="UTF-8"?>
<rdf:RDF
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:foaf="http://xmlns.com/foaf/0.1/"
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:dcterms="http://purl.org/dc/terms/"
xmlns:q="http://quote.lindenb.org/ontology"
>
<foaf:Person rdf:about="http://en.wikipedia.org/wiki/Isaac_Asimov">
<foaf:name>Isaac Asimov</foaf:name>
</foaf:Person>

<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="en">Violence is the last refuge of the incompetent.</q:quote>
</q:Quote>

<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="en">Never let your sense of morals prevent you from doing what is right.</q:quote>
</q:Quote>

<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="en">Creationists make it sound as though a "theory" is something you dreamt up after being drunk all night</q:quote>
</q:Quote>


<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="fr">Si le savoir peut créer des problèmes, ce n'est pas l'ignorance qui les résoudra.</q:quote>
</q:Quote>

</rdf:RDF>


Append the file:
mergexml --replace -r -i database.rdf chunk.rdf

Result:
$ xmllint --format database.rdf
<?xml version="1.0" encoding="UTF-8"?>
<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:foaf="http://xmlns.com/foaf/0.1/" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:q="http://quote.lindenb.org/ontology">
<foaf:Person rdf:about="http://en.wikipedia.org/wiki/Louis_Pasteur">
<foaf:name>Louis Pasteur</foaf:name>
</foaf:Person>
<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Louis_Pasteur"/>
<q:quote xml:lang="en">Science knows no country, because knowledge belongs to humanity, and is the torch which illuminates the world.</q:quote>
</q:Quote>
<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Louis_Pasteur"/>
<q:quote xml:lang="en">I am utterly convinced that Science and Peace will triumph over Ignorance and War, that nations will eventually unite not to destroy but to edify, and that the future will belong to those who have done the most for the sake of suffering humanity.</q:quote>
</q:Quote>
<foaf:Person rdf:about="http://en.wikipedia.org/wiki/Isaac_Asimov">
<foaf:name>Isaac Asimov</foaf:name>
</foaf:Person>
<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="en">Violence is the last refuge of the incompetent.</q:quote>
</q:Quote>
<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="en">Never let your sense of morals prevent you from doing what is right.</q:quote>
</q:Quote>
<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="en">Creationists make it sound as though a "theory" is something you dreamt up after being drunk all night</q:quote>
</q:Quote>
<q:Quote>
<q:author rdf:resource="http://en.wikipedia.org/wiki/Isaac_Asimov"/>
<q:quote xml:lang="fr">Si le savoir peut créer des problèmes, ce n'est pas l'ignorance qui les résoudra.</q:quote>
</q:Quote>
</rdf:RDF>


That's it

Pierre

13 August 2011

A FUSE-based filesystem reproducing the NCBI Taxonomy hierarchy.

In this post I will show how I've used FUSE to create a new loadable filesystem for Linux, reproducing the tree of the NCBI Taxonomy.

From the Fuse Homepage:With FUSE it is possible to implement a fully functional filesystem in a userspace program. Features include:

  • Simple library API
  • Simple installation (no need to patch or recompile the kernel)
  • Secure implementation
  • Userspace - kernel interface is very efficient
  • Usable by non privileged users
  • Runs on Linux kernels 2.4.X and 2.6.X
  • Has proven very stable over time
.

Downloading and installing Fuse

Download Fuse from sourceforge:http://sourceforge.net/projects/fuse/files/fuse-2.X/.

$ tar xvfz fuse-2.8.5.tar.gz

$ cd fuse-2.8.5/
$ ./configure
$ make
$ make install


The classes


Taxon

A Taxon is a node in the NCBI Taxonomy:

class Taxon
{
public:
/* node id */
int id;
/* parent id */
int parent_id;
/* number of children */
int count_children;
/* name for this node */
char* name;
/* node children */
Taxon** children;
Taxon();
~Taxon();
/* compare by id */
static int comparator(const void* p1,const void* p2);
/* compare by name*/
static int compareByName(const void* p1,const void* p2);
const Taxon* getChildrenAt(int i) const;
/* find children by its name*/
const Taxon* findChildByName(const char* name) const;
/* recursive. find child from a unix path */
const Taxon* findChildByPath(const char* path) const;
/* stupid representation of a node as a XML file */
string xml() const;
};

FuseTaxonomy

The NCBI Taxonomy:
class FuseTaxonomy

{
private:
/* number of nodes */
size_t nTaxons;
/** taxons ordered by ids */
Taxon** taxons;
/** taxons ordered by names */
Taxon** names;
public:
/* global instance */
static FuseTaxonomy* INSTANCE;
FuseTaxonomy();
~FuseTaxonomy();
/* find taxon by ID */
const Taxon* findTaxonById(int id)const;
/* find taxon by name */
const Taxon* findTaxonByName(const char* name) const;
/* root of taxonomy */
const Taxon* getRoot() const;
/* read file nodes.dmp */
void readNodes(const char* nodes);
/* read file names.dmp */
void readNames(const char* namesfile);
/* find taxon node from unix path */
const Taxon* findByPath(const char* path) const;
/** FUSE CALLBACK: This function returns metadata concerning a file specified by path in a special stat structure. */
static int getattr(const char *path, struct stat *stbuf);
/* FUSE CALLBACK: used to read directory contents */
static int readdir(const char *path, void *buf, fuse_fill_dir_t filler, off_t offset, struct fuse_file_info *fi);
/* FUSE CALLBACK: checks whatever user is permitted to open the /hello file with flags given in the fuse_file_info structure. */
static int open(const char *path, struct fuse_file_info *fi);
/* FUSE CALLBACK: used to feed the user with data from the file. */
static int read(const char *path, char *buf, size_t size, off_t offset, struct fuse_file_info *fi);
};

The static functions 'getattr', 'readdir', 'open' and 'read' are the callbacks called by the FUSE API to explore the new filesystem and must be initialized in the 'main' method.

Test


Load the NCBI taxonomy
$ wget -O taxdump.tar.gz "ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz"

$ tar xvfz taxdump.tar.gz names.dmp
$ tar xvfz taxdump.tar.gz nodes.dmp

Create a temporary folder for our new filesystem:
$ mkdir -p tmp_fuse

load the NCBI taxonomy and install the new filesystem
$ ./fusetaxonomy nodes.dmp names.dmp  tmp_fuse

What's is in this new filesystem ?
$ls tmp_fuse

root

And what are the nodes under "root" ?
$ls tmp_fuse

cellular_organisms other_sequences unclassified_sequences Viroids Viruses

And what's under 'Eukaryota' ?
ls -la tmp_fuse/root/cellular_organisms/Eukaryota
total 0
drwxr-xr-x 22 root root 0 1970-01-01 01:00 .
drwxr-xr-x 3 root root 0 1970-01-01 01:00 ..
drwxr-xr-x 9 root root 0 1970-01-01 01:00 Alveolata
drwxr-xr-x 11 root root 0 1970-01-01 01:00 Amoebozoa
drwxr-xr-x 4 root root 0 1970-01-01 01:00 Apusozoa
drwxr-xr-x 5 root root 0 1970-01-01 01:00 Centroheliozoa
drwxr-xr-x 4 root root 0 1970-01-01 01:00 Cryptophyta
drwxr-xr-x 38 root root 0 1970-01-01 01:00 environmental_samples
drwxr-xr-x 5 root root 0 1970-01-01 01:00 Euglenozoa
drwxr-xr-x 7 root root 0 1970-01-01 01:00 Fornicata
drwxr-xr-x 5 root root 0 1970-01-01 01:00 Glaucocystophyceae
drwxr-xr-x 11 root root 0 1970-01-01 01:00 Haptophyceae
drwxr-xr-x 4 root root 0 1970-01-01 01:00 Heterolobosea
drwxr-xr-x 4 root root 0 1970-01-01 01:00 Jakobida
drwxr-xr-x 3 root root 0 1970-01-01 01:00 Katablepharidophyta
drwxr-xr-x 1 root root 0 1970-01-01 01:00 Malawimonadidae
drwxr-xr-x 6 root root 0 1970-01-01 01:00 Opisthokonta
drwxr-xr-x 8 root root 0 1970-01-01 01:00 Oxymonadida
drwxr-xr-x 9 root root 0 1970-01-01 01:00 Parabasalia
drwxr-xr-x 9 root root 0 1970-01-01 01:00 Rhizaria
drwxr-xr-x 7 root root 0 1970-01-01 01:00 Rhodophyta
drwxr-xr-x 25 root root 0 1970-01-01 01:00 stramenopiles
drwxr-xr-x 36 root root 0 1970-01-01 01:00 unclassified_eukaryotes
drwxr-xr-x 3 root root 0 1970-01-01 01:00 Viridiplantae

Let's use find !
find tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/ | head

tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha/Astrophorida
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha/Astrophorida/Geodiidae
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha/Astrophorida/Geodiidae/Geodia
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha/Astrophorida/Geodiidae/Geodia/Geodia_cydonium
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha/Astrophorida/Geodiidae/Geodia/Geodia_neptuni
tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Porifera/Demospongiae/Tetractinomorpha/Astrophorida/Geodiidae/Geodia/Geodia_papyracea

What is the content of Homo_sapiens_neanderthalensis ?:
$ cat tmp_fuse/root/cellular_organisms/Eukaryota/Opisthokonta/Metazoa/Eumetazoa/Bilateria/Coelomata/Deuterostomia/Chordata/Craniata/Vertebrata/Gnathostomata/Teleostomi/Euteleostomi/Sarcopterygii/Tetrapoda/Amniota/Mammalia/Theria/Eutheria/Euarchontoglires/Primates/Haplorrhini/Simiiformes/Catarrhini/Hominoidea/Hominidae/Homininae/Homo/Homo_sapiens/Homo_sapiens_neanderthalensis

<?xml version="1.0"?>
<Taxon-Id>63221</Taxon-Id>

unmount the NCBI filesystem
sudo ${FUSEDIR}/bin/fusermount -u  tmp_fuse

The Full source code.

The full source code is available as a 'gist' on github.com.
The code is also on github here

That's it !

Pierre

10 August 2011

Wikipedia OpenSearch and semantic bash completion

I wrote a tool that use the opensearch API for mediawiki to display the results of a query in wikipedia:

$ opensearch Charles Darw
Charles_Darwin Charles Robert Darwin FRS (12 February 1809 – 19 April 1882) was an English...
Charles_Darwin_University Charles Darwin University (CDU) is an Australian public university wi...
Charles_Darwin%27s_health For much of his adult life, Charles Darwin's health was repeatedly co...
Charles_Darwin_Foundation The Charles Darwin Foundation was founded in 1959, under the auspices...
Charles_Darwin_National_Park Charles Darwin National Park is in the Northern Territory of Austr...
Charles_Darwin_Research_Station The Charles Darwin Research Station (CDRS) is a biological rese...
Charles_Darwin%27s_education Charles Darwin's education gave him a foundation in the doctrine o...
Charles_Darwin_%281758%E2%80%931778%29 Charles Darwin (3 September 1758–15 May 1778) was the ...
Charles_Darwin%27s_religious_views Charles Darwin's views on religion have been the subject of ...
Charles_Darwin_and_the_Tree_of_Life Charles Darwin and the Tree of Life is a 2009 television do...
Charles_Darwin_School Charles Darwin School is the only secondary school in the Biggin Hill area.
Charles_Darwin_%28aviator%29 Major Charles John Wharton Darwin was a First World War flying ace...
Charles_Darwin_Reserve Charles Darwin Reserve is a 686 km2 nature reserve in Western Australia.
Charles_Darwin_bibliography This is a partial list of the writings of Charles Darwin, including...
Darwin Darwin may refer to:


This tool is available on github at:

Bash completion

Now I want to use this tool for bash completion. Say, I wrote a tool named 'foo' that expects the title of a wikipedia article as an argument. In ~/.bash_completion.d/, I wrote a file named opensearch:
_opensearch()

{
local cur choice reply
COMPREPLY=()
cur="${COMP_WORDS[COMP_CWORD]}"
choice=$(opensearch -d -c ${cur})
reply=$(compgen -W "${choice}" -- ${cur})
COMPREPLY=( $reply )
return 0
} &&

complete -F _opensearch foo

This script is activated when a new terminal is loaded in ${HOME}/.bashrc:
(...)
if [ -f ${HOME}/.bash_completion.d/opensearch ]; then
. ${HOME}/.bash_completion.d/opensearch
fi
Now, typing 'foo', a keyword and <TAB> activates this new completion:
$ foo Charles_D <TAB> <TAB>

Charles_Dance Charles_Darwin_University Charles_Doolittle_Walcott Charles_Dutoit
Charles_Darwin Charles_Dickens Charles_Durning
$ foo Charles_D


That's it,

Pierre

09 August 2011

Quick tip: bash completion for Bioinformatics

The default behavior for the <tab> completion can be extended by creating a new file:


${HOME}/.bash_completion
in your home.
For example, you can write your own Bash Completion for samtools, bwa or your favourite tool by adding the following line in "${HOME}/.bash_completion":

complete -f -X '!*.@(bam|sam|fasta|fa|fa.gz|fastq.gz|fasta.gz)' samtools bwa

Open a new bash and type:

$ samtools index <tab>
myfile.fastq.gz sequence.fasta


Pierre

01 August 2011

Using the Freebase and the Bioportal Widgets to create a semantic object.


The following HTML code uses:
After completion, it generates a JSON object describing a semantic object:
{
"subject":"http://www.freebase.com/view/en/nsp3",
"predicate":"http://purl.org/obo/owl/MI#MI_0407",
"value":"http://www.freebase.com/view/en/roxan",
"pmid":15047801
}



Source


You can download the source and test it:


That's it,

Pierre

Memory-Mapping the Human Genome with 'mmap': my notebook

In this post, I've explored how to use a memory-mapped file to handle the 3Go of the Human Genome sequence as if it was entirely loaded in memory.
Via wikipedia: "A memory-mapped file is a segment of virtual memory which has been assigned a direct byte-for-byte correlation with some portion of a file or file-like resource. This resource is typically a file that is physically present on-disk... In computing, mmap is a POSIX-compliant Unix system call that maps files or devices into memory. It is a method of memory-mapped file I/O. It naturally implements demand paging, because initially file contents are not entirely read from disk and do not use physical RAM at all. The actual reads from disk are performed in "lazy" manner, after a specific location is accessed."
Using a C++ program, I'm going to memory-map the fasta sequence of the human genome (indexed with samtools faidx) and search the position of some short-reads using a BoyerMoore algorithm:

Required members:


/* maps a chromosome to its samtools faidx index */
map<string,faidx1_t> name2index;
/* used to get the size of the file */
struct stat buf;
/* genome fasta file file descriptor */
int fd;
/* the mmap (memory mapped) pointer */
char *mapptr;

Opening the mmap

string faidx(fastaFile);
string line;
faidx.append(".fai");
/* open *.fai file */
ifstream in(faidx.c_str(),ios::in);
/* read indexes in the .fai file that was created with samtools */
while(getline(in,line,'\n'))
{
faidx1_t index;
//parse the faidx line...
//...
name2index.insert(make_pair(chrom,index));
}
/* close index file */
in.close();

/* get the whole size of the fasta file */
stat(fasta, &buf);
/* open the fasta file */
fd = open(fastaFile, O_RDONLY);
/* open a memory mapped file associated to this fasta file descriptor */
mapptr = (char*)mmap(0, buf.st_size, PROT_READ, MAP_SHARED, fd, 0);

It's a kind of MAGIC: Getting the base at index 'position-th' of chromosome 'chrom'


std::map<std::string,faidx1_t>::iterator r=name2index.find(chrom);
faidx1_t& index=r->second;
char base=at(&index,position);
(...)
/* returns the base at position 'index' for the chromosome indexed by faidx */
char at(const FaidxPtr faidx,int64_t index)
{
long pos= faidx->offset +
index / faidx->line_blen * faidx->line_len +
index % faidx->line_blen
;
/* here is the magic: no need to fseek/fread/ftell the file */
return toupper(mapptr[pos]);
}

Mapping the short reads

I've hacked a simple Boyer-Moore-Horspool algorithm from ttp://en.wikipedia.org/wiki/Boyer-Moore-Horspool_algorithm. Of course, you wouldn't use this algorithm to map your short reads for real :-) .

Disposing the mmap

/* close memory mapped map */
if(mapptr!=NULL) munmap(mapptr,buf.st_size);
/* dispose fasta file descriptor */
if(fd!=-1) close(fd);

Compile & run


$ g++ -Wall testmmap.cpp -lz

$ ./a.out -g /path/tp/hg19.fa /path/to/my.fastq.gz

ATCATTTTCCTCCTAACAGATTAAAAATCAAGAAATATAAACCAGATGTAGCAG chr11 93065362
(...)

Source code





That's it,

Pierre

31 July 2011

Storing some SNPs using the leveldb library. My notebook.

In this post I'll describe how I've used the leveldb API, a C++ key-value database, to store some SNPs (key=rs,value=sequence). "LevelDB is a fast key-value storage library written at Google that provides an ordered mapping from string keys to string values.". A benchmark for this engine is available here: http://code.google.com/p/leveldb/source/browse/trunk/doc/benchmark.html.

Download & install

$ svn checkout http://leveldb.googlecode.com/svn/trunk/ leveldb-read-only
$ cd leveldb-read-only/
$ make

Open & close a leveldb database

#include "leveldb/db.h"
(...)
leveldb::DB* db=NULL;
leveldb::Options options;
options.comparator=&rs_comparator;/* custom comparator for ordering the keys, see later */
options.create_if_missing = true;
clog << "[LOG]opening database" << endl;
leveldb::Status status = leveldb::DB::Open(options,db_home, &db);
//check status
if (!status.ok())
{
cerr << "cannot open " << db_home << " : " << status.ToString() << endl;
return EXIT_FAILURE;
}
(...)
/* use the database */
(...)
delete db;

A custom comparator for ordering the keys

Here, the keys are the rs-ids ordered on their numerical value. Both keys and values have a type "leveldb::Slice".
class RsComparator : public leveldb::Comparator
{
private:
/* parses the rsId: rs[0-9]+ */
int rs(const leveldb::Slice& s) const
{
int n=0;
for(size_t i=2;i< s.size();++i)
{
n=n*10+s[i]-'0';
}
return n;
}
public:
virtual int Compare(const leveldb::Slice& a, const leveldb::Slice& b) const
{
return rs(a)-rs(b);
}
(...)
}
(...)
RsComparator rs_comparator;
options.comparator=&rs_comparator;

Inserting a FASTA sequence


(...)
std::string name;
std::string sequence;
(...)
leveldb::Status status = db->Put(leveldb::WriteOptions(), name, sequence);
if(!status.ok())
{
cerr << "Cannot insert "<< name << " "
<< status.ToString() << endl;
return EXIT_FAILURE;
}

Searching for a rs###

(...)
std::string name(argv[optind]);
std::string sequence;
(...)
leveldb::Status status = db->Get(leveldb::ReadOptions(),seqname, &sequence);
if(!status.ok())
{
cerr << "Cannot find " << seqname<< " in " << db_home << endl;
continue;
}
printFasta(seqname,sequence);
(...)

Dumping all the SNPs using an iterator

(...)
leveldb::Iterator* it = db->NewIterator(leveldb::ReadOptions());
for (it->SeekToFirst(); it->Valid(); it->Next())
{
printFasta(it->key(),it->value());
}
delete it;
(...)

Examples

Reading the SNPs:
$ curl -s "ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/rs_fasta/rs_chMT.fas.gz" |\
./rsput -D snp.db

[LOG]opening database
[LOG] added rs8896
[LOG] added rs8936
[LOG] added rs9743
(...)
[LOG]closing database

$ du -h snp.db
336K snp.db
Dump the snps

$ ./rsget  -D snp.db

[LOG]opening database
>rs8896
GGTGTTGGTTCTCTTAATCTTTAACTTAAAAGGTTAATGCTAAGTTAGCTTTACAGTGGG
CTCTAGAGGGGGTAGAGGGGGTGYTATAGGGTAAATACGGGCCCTATTTCAAAGATTTTT
AGGGGAATTAATTCTAGGACGATGGGCATGAAACTGTGGTTTGCTCCACAGATTTCAGAG
CATT
>rs8936
ACTACGGCGGACTAATCTTCAACTCCTACATACTTCCCCCATTATTCCTAGAACCAGGCG
ACCTGCGACTCCTTGACGTTGACAATCGAGTAGTACTCCCGATTGAAGCCCCCATTCGTA
TAATAATTACATCACAAGACGTCTTGCACTCATGAGCTGTCCCCACATTAGGCTTAAAAA
CAGATGCAATTCCCGGACGTHTAAACCAAACCACTTTCACCGCTACACGACCGGGGGTAT
ACTACGGTCAATGCTCTGAAATCTGTGGAGCAAACCACAGTTTCATGCCCATCGTCCTAG
AATTAATTCCCCTAAAAATCTTTGAAATAGGGCCCGTATTTACCCTATAGCACCCCCTCT
ACCCCCTCTAGAGCCCACTGTAAAGCTAACTTAGCATTAAC
>rs9743
CCATGTGATTTCACTTCCACTCCATAACGCTCCTCATACTAGGCCTACTAACCAACACAC
TAACCATATACCAATGATGNCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACA
CACCACCTGTCCAAAAAGGCCTTCGATACGGGATAATCCTATTTATTACCTCAGAANTTT
TTTTCTTCGCAGGATTTTTCTGAGCCTTTTACCACTCCAGCCTAGCCCCTACCCCCCAAN
(...)
[LOG]closing database
Search for some SNPs:
$ ./rsget -D snp.db rs78894381 rs72619361 rs25
[LOG]opening database
[LOG]searching for rs78894381
>rs78894381
CTACTAATCTCATCAACACAACCCCCGCCCATCCTACCCAGCACACACACACCGCTGCTA
ACCCCATACCCCGAACCAACCAAACCCCAAAGACACCCCCNCACAGTTTATGTAGCTTAC
CTCCTCAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATA
GGTTTGGTCCTAGCCTTTCTA
[LOG]searching for rs72619361
>rs72619361
ATGCATTTGGTATTTTAATCTGGGGGGTGTGCACGCGATAGCATTGTGAAACGCTGGCCC
CAGAGCACCCTATGTCGCAGTGTCTGTCTTTGATTCCTGCCYCATCCCATTATTGATCAC
ACCTACATTCAATATCCCAGGCGAGCATACCTATCACAAGGTGTTAATTAATTAATGCTT
GTAGGACATAACAATCAGTAAAC
[LOG]searching for rs25
Cannot find rs25 in snp.db
[LOG]closing database

Source code






That's it,

Pierre

27 July 2011

Controlling IGV through a Port: my notebook.

A protocol for interacting with the The Integrative Genomics Viewer (IGV) is described here: http://www.broadinstitute.org/igv/PortCommands. So, you can save some pictures of a genomic region without interacting with IGV.

Via a shell script


$ cat script.sh

mkdir -p /tmp/igv
sleep 2
echo "snapshotDirectory /tmp/igv"
sleep 2
echo "goto chr3:100"
sleep 2
echo "snapshot"
sleep 2
echo "goto chr3:200-300"
sleep 2
echo "snapshot"
sleep 2

invoke the script
$ sh script.sh |  telnet 127.0.0.1 60151
Trying 127.0.0.1...
Connected to srv-clc-02.irt.univ-nantes.prive (127.0.0.1).
Escape character is '^]'.
INFO [2011-07-27 17:58:16,680] [CommandExecutor.java:124] [Thread-6] OK
OK
INFO [2011-07-27 17:58:18,685] [CommandExecutor.java:124] [Thread-6] OK
OK
INFO [2011-07-27 17:58:18,687] [MacroSnapshotAction.java:85] [Thread-6] Snapshot: chr3_80_119.png
INFO [2011-07-27 17:58:20,787] [CommandExecutor.java:124] [Thread-6] OK
OK
INFO [2011-07-27 17:58:22,792] [CommandExecutor.java:124] [Thread-6] OK
OK
INFO [2011-07-27 17:58:22,794] [MacroSnapshotAction.java:85] [Thread-6] Snapshot: chr3_200_300.png
Connection closed by foreign host.


Via java


The code
import java.io.*;
import java.net.*;
class Test
{
public static void main(String[] args) throws Exception
{
Socket socket = new Socket("127.0.0.1", 60151);
PrintWriter out = new PrintWriter(socket.getOutputStream(), true);
BufferedReader in = new BufferedReader(new InputStreamReader(socket.getInputStream()));


out.println("snapshotDirectory /tmp/igv");
System.out.println(in.readLine());
out.println("goto chr3:100");
System.out.println(in.readLine());

out.println("snapshot");
System.out.println(in.readLine());

in.close();
out.close();
socket.close();
}
}

Compile and execute:
$ javac Test
$ java Test
INFO [2011-07-27 18:01:15,706] [CommandExecutor.java:124] [Thread-6] OK
OK

INFO [2011-07-27 18:01:17,711] [CommandExecutor.java:124] [Thread-6] OK
OK

INFO [2011-07-27 18:01:17,729] [MacroSnapshotAction.java:85] [Thread-6] Snapshot: chr3_80_119.png
INFO [2011-07-27 18:01:20,533] [CommandExecutor.java:124] [Thread-6] OK
OK

26 July 2011

A mysql full-text parser searching for some SNPs

A mysql full-text parser server plugin can be used to replace or modify the built-in full-text parser.
This post a full-text parser plugin named bioparser that extracts the rs#### id from a text.

The source

The source of this plugin is available on github at:
https://github.com/lindenb/bioudf/blob/master/src/bioparser.c.

The work horse of the plugin is a simple function bioparser_parse scanning the SQL query or the TEXT. Each time a word starting with "rs" and followed by one or more number is found, it is added to the list of words to be find:
static int bioparser_parse(MYSQL_FTPARSER_PARAM *param)
{
char *curr=param->doc;
const char *begin=param->doc;
const char *end= begin + param->length;

param->flags = MYSQL_FTFLAGS_NEED_COPY;
while(curr+2<end)
{
if(tolower(*curr)=='r' &&
tolower(*(curr+1))=='s' &&
isdigit(*(curr+2)) &&
(curr==begin || IS_DELIM(*(curr-1) ) )
)
{
char* p=curr+2;
while(p!=end && isdigit(*p))
{
++p;
}
if(p==end || IS_DELIM(*p))
{
my_add_word(param,curr,p-curr);
}
curr=p;
}
else
{
curr++;
}
}

return 0;
}

Install the Plugin


mysql> INSTALL PLUGIN bioparser SONAME 'bioparser.so';
Query OK, 0 rows affected (0.00 sec)

mysql> show plugins;
+-----------------------+----------+--------------------+--------------+---------+
| Name | Status | Type | Library | License |
+-----------------------+----------+--------------------+--------------+---------+
(...)
| partition | ACTIVE | STORAGE ENGINE | NULL | GPL |
| bioparser | ACTIVE | FTPARSER | bioparser.so | GPL |
+-----------------------+----------+--------------------+--------------+---------+
21 rows in set (0.00 sec)

Invoke Plugin


create a table that will use the plugin:
mysql> create table pubmed(
abstract TEXT,
FULLTEXT (abstract) WITH PARSER bioparser
) ENGINE=MyISAM;

Insert some abstracts.
mysql> insert into pubmed(abstract) values("A predictive role in radiation pneumonitis (RP) development was observed for the LIG4 SNP rs1805388 (adjusted hazard ratio, 2.08; 95% confidence interval, 1.04-4.12; P = .037 for the CT/TT genotype vs the CC genotype). In addition, men with the TT genotype of the XRCC4 rs6869366 SNP and women with AG + AA genotypes of the XRCC5 rs3835 SNP also were at increased risk of developing severe RP.");
Query OK, 1 row affected (0.00 sec)

(...)

mysql> select abstract from pubmed;
+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| abstract |
+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| A predictive role in radiation pneumonitis (RP) development was observed for the LIG4 SNP rs1805388 (adjusted hazard ratio, 2.08; 95% confidence interval, 1.04-4.12; P = .037 for the CT/TT genotype vs the CC genotype). In addition, men with the TT genotype of the XRCC4 rs6869366 SNP and women with AG + AA genotypes of the XRCC5 rs3835 SNP also were at increased risk of developing severe RP. |
| Nonhomologous end joining (NHEJ) is a pathway that repairs DNA double-strand breaks (DSBs) to maintain genomic stability in response to irradiation. The authors hypothesized that single nucleotide polymorphisms (SNPs) in NHEJ repair genes may affect clinical outcomes in patients with nonsmall cell lung cancer (NSCLC) who receive definitive radio(chemo)therapy. |
| The authors genotyped 5 potentially functional SNPs-x-ray repair complementing defective repair in Chinese hamster cells 4 (XRCC4) reference SNP (rs) number rs6869366 (-1394 guanine to thymine [-1394G?T] change) and rs28360071 (intron 3, deletion/insertion), XRCC5 rs3835 (guanine to adenine [G?A] change at nucleotide 2408), XRCC6 rs2267437 (-1310 cytosine to guanine [C?G) change], and DNA ligase IV (LIG4) rs1805388 (threonine-to-isoleucine change at codon 9 [T9I])-and estimated their associations with severe radiation pneumonitis (RP) (grade ?3) in 195 patients with NSCLC. |
| The current results indicated that NHEJ genetic polymorphisms, particularly LIG4 rs1805388, may modulate the risk of RP in patients with NSCLC who receive definitive radio(chemo)therapy. Large studies will be needed to confirm these findings. |
| The repair of DNA double-strand breaks (DSBs) is the major mechanism to maintain genomic stability in response to irradiation. We hypothesized that genetic polymorphisms in DSB repair genes may affect clinical outcomes among non-small cell lung cancer (NSCLC) patients treated with definitive radio(chemo)therapy. |
| We also found that RAD51 -135G>C and XRCC2 R188H SNPs were independent prognostic factors for overall survival (adjusted HR?=?1.70, 95% CI, 1.14-2.62, P?=?0.009 for CG/CC vs. GG; and adjusted HR?=?1.70; 95% CI, 1.02-2.85, P?=?0.043 for AG vs. GG, respectively) and that the SNP-survival association was most pronounced in the presence of RP. |
| A total of 291 patients (145 male/146 female, mean age (± S.D.) 52.2 (± 13.1) years) with PsA were examined clinically, by standard laboratory tests and their DNA was genotyped for the SNP rs2476601 (PTPN22 +1858 C/T). Allelic frequencies were determined and compared with 725 controls. |
| this is a test rs2476601, rs1805388, rs3835 and rs25 |
+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
8 rows in set (0.00 sec)

Now, test the plugin:

mysql> select
concat(left(abstract,40),"...") as ABSTRACT,
match(abstract) against("DNA double-strand) as SCORE
from pubmed group by 2 HAVING SCORE!=0;
Empty set (0.00 sec)

mysql> select
concat(left(abstract,40),"...") as ABSTRACT,
match(abstract) against("rs25") as SCORE
from pubmed group by 2 HAVING SCORE!=0;
+---------------------------------------------+--------------------+
| ABSTRACT | SCORE |
+---------------------------------------------+--------------------+
| this is a test rs2476601, rs1805388, rs3... | 1.8603347539901733 |
+---------------------------------------------+--------------------+
1 row in set (0.01 sec)

mysql> select
concat(left(abstract,40),"...") as ABSTRACT,
match(abstract) against("rs2476601 rs1805388 rs6869366") as SCORE
from pubmed group by 2 HAVING SCORE!=0 order by 2 desc;
+---------------------------------------------+--------------------+
| ABSTRACT | SCORE |
+---------------------------------------------+--------------------+
| A total of 291 patients (145 male/146 fe... | 1.086121916770935 |
| A predictive role in radiation pneumonit... | 1.0619741678237915 |
| this is a test rs2476601, rs1805388, rs3... | 1.0502985715866089 |
| The authors genotyped 5 potentially func... | 1.0388768911361694 |
+---------------------------------------------+--------------------+
4 rows in set (0.00 sec)

mysql> select
concat(left(abstract,40),"...") as ABSTRACT,
match(abstract) against("rs25,rs2476601,rs1805388,rs6869366") as SCORE
from pubmed group by 2 HAVING SCORE!=0 order by 2 desc;
+---------------------------------------------+--------------------+
| ABSTRACT | SCORE |
+---------------------------------------------+--------------------+
| this is a test rs2476601, rs1805388, rs3... | 2.9106333255767822 |
| A total of 291 patients (145 male/146 fe... | 1.086121916770935 |
| A predictive role in radiation pneumonit... | 1.0619741678237915 |
| The authors genotyped 5 potentially func... | 1.0388768911361694 |
+---------------------------------------------+--------------------+
4 rows in set (0.00 sec)

uninstall the plugin


mysql> UNINSTALL PLUGIN bioparser;
Query OK, 0 rows affected (0.00 sec)


That's it,

Pierre

25 July 2011

The samtools C API: creating a WIG file for the genome coverage

I wrote a tool which uses the samtools API to generate a WIG file for the genome coverage. This tool, bam2wig, as well as some other tools using this library are available in the following git repository:


My code is based on this original C source hosted on: http://samtools.sourceforge.net/sam-exam.shtml. (I hope I understood this API as many functions/structures remain undocumented).

Example

#compile
$ export SAMDIR=/path/to/samtools-0.1.17
$ cd amtools-utilities
$ make
#invoke bam2wig
$ .bin/bam2wig -t /path/to/my/sample.bam "chr1:120408225-120558920" | head
track name="__TRACK_NAME__" description="__TRACK_DESC__" type="wiggle_0"
fixedStep chrom=chr1 start=120408225 step=1 span=1
1
1
1
2
2
3
4
5

Result

The wig file uploaded as a custom track in the UCSC genome browser:


That's it

Pierre