From a262ac1754f53447daf71a4281a73b1e7f29adfb Mon Sep 17 00:00:00 2001 From: Anton Nekrutenko Date: Thu, 20 Mar 2008 19:38:31 +0000 Subject: [PATCH] Finalizing taxonomy. Removed taxonomy processing from scripts. This does not require sqlite anymore. Script for fetching taxonomy now uses flat files directly and is much faster --- scripts/taxonomy/processTaxonomy.sh | 11 +- scripts/taxonomy/process_NCBI_taxonomy.py | 50 -- test-data/taxonomyGI.dat | 5 - tools/taxonomy/T2PS/t2ps.readme | 21 +- tools/taxonomy/T2T/t2t.readme | 20 +- tools/taxonomy/TB/avl.c | 890 ++++++++++++++++++++++ tools/taxonomy/TB/avl.h | 115 +++ tools/taxonomy/TB/main.c | 788 +++++++++++++++++++ tools/taxonomy/TB/tb.readme | 21 + tools/taxonomy/find_diag_hits.py | 6 +- tools/taxonomy/find_diag_hits.xml | 2 +- tools/taxonomy/gi2taxonomy.py | 172 +++++ tools/taxonomy/gi2taxonomy.xml | 12 +- tools/taxonomy/tax.py | 73 -- 14 files changed, 2038 insertions(+), 148 deletions(-) delete mode 100755 scripts/taxonomy/process_NCBI_taxonomy.py delete mode 100644 test-data/taxonomyGI.dat create mode 100644 tools/taxonomy/TB/avl.c create mode 100644 tools/taxonomy/TB/avl.h create mode 100644 tools/taxonomy/TB/main.c create mode 100644 tools/taxonomy/TB/tb.readme create mode 100644 tools/taxonomy/gi2taxonomy.py delete mode 100644 tools/taxonomy/tax.py diff --git a/scripts/taxonomy/processTaxonomy.sh b/scripts/taxonomy/processTaxonomy.sh index c6cf7a76b04..95c6c973a69 100755 --- a/scripts/taxonomy/processTaxonomy.sh +++ b/scripts/taxonomy/processTaxonomy.sh @@ -1,5 +1,3 @@ -PYTHONPATH="../../lib:../../eggs:../../eggs/`../check_python_ucs.py`" -export PYTHONPATH echo "Getting files from NCBI..." wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid_nucl.dmp.gz @@ -9,9 +7,8 @@ gunzip -c taxdump.tar.gz | tar xvf - gunzip gi_taxid_nucl.dmp.gz gunzip gi_taxid_prot.dmp.gz cat gi_taxid_nucl.dmp gi_taxid_prot.dmp > gi_taxid_all.dmp -rm gi_taxid_nucl.dmp gi_taxid_prot.dmp -echo "Parsing names.dmg" -cat names.dmp | cut -f 1,2,4 -d "|" | tr -s "\t" "|" | tr "|" "\t" | sed s/\"//g > names.txt -python process_NCBI_taxonomy.py gi_taxid_all.dmp names.txt taxonomy.db -echo "Done!.." +echo "Sorting gi2tax files..." +sort -n -k 1 gi_taxid_all.dmp > gi_taxid_sorted.txt +rm gi_taxid_nucl.dmp gi_taxid_prot.dmp gi_taxid_all.dmp + diff --git a/scripts/taxonomy/process_NCBI_taxonomy.py b/scripts/taxonomy/process_NCBI_taxonomy.py deleted file mode 100755 index 1f73449b2e0..00000000000 --- a/scripts/taxonomy/process_NCBI_taxonomy.py +++ /dev/null @@ -1,50 +0,0 @@ -""" -process_NCBI_taxonomy.py -""" - -import pkg_resources -pkg_resources.require( 'pysqlite' ) -from pysqlite2 import dbapi2 as sqlite -import string, sys, tempfile - -def stop_err(msg): - sys.stderr.write(msg) - sys.exit() - - -try: - gi2tax = open(sys.argv[1], 'r') - names = open(sys.argv[2], 'r') - db_name = sys.argv[3] -except: - stop_err('Check arguments: process_NCBI_taxonomy.py \n') - - -try: - con = sqlite.connect(db_name) - cur = con.cursor() - cur.execute('create table gi2tax(gi int unsigned not null, taxId int unsigned not null)') - cur.execute('create table t_names(taxId int unsigned not null, name text not null)') - cur.execute('create table names(taxId int unsigned not null, name text not null)') - - for line in gi2tax: - fields = string.split(line.rstrip(), '\t') - cur.execute('insert into gi2tax values(%s, %s)' % ( fields[0], fields[1] ) ) - - gi2tax.close() - - for line in names: - fields = string.split(line.rstrip(), '\t') - cur.execute('insert into t_names values(%s, "%s")' % ( fields[0], fields[1] ) ) - - names.close() - - cur.execute('create index gi_i on gi2tax(gi)') - cur.execute('insert into names select * from t_names group by name') - cur.execute('drop table t_names') - cur.execute('create index name_i on names(name)') - cur.execute('vacuum') - con.commit() - con.close() -except Exception, e: - stop_err("%s\n" % e) diff --git a/test-data/taxonomyGI.dat b/test-data/taxonomyGI.dat deleted file mode 100644 index d246b34882b..00000000000 --- a/test-data/taxonomyGI.dat +++ /dev/null @@ -1,5 +0,0 @@ -33001686 9443 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n Euarchontoglires Primates n n n n n n n n n n -23236241 9604 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n Euarchontoglires Primates Haplorrhini Hominoidea Hominidae n n n n n n n -12583 9606 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n Euarchontoglires Primates Haplorrhini Hominoidea Hominidae n n n Homo n Homo sapiens n -410771 40674 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n n n n n n n n n n n n n -2286205 63221 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n Euarchontoglires Primates Haplorrhini Hominoidea Hominidae n n n Homo n Homo sapiens Homo sapiens neanderthalensis diff --git a/tools/taxonomy/T2PS/t2ps.readme b/tools/taxonomy/T2PS/t2ps.readme index d1f65cf32ba..aece1e1716a 100644 --- a/tools/taxonomy/T2PS/t2ps.readme +++ b/tools/taxonomy/T2PS/t2ps.readme @@ -1,5 +1,22 @@ -To compile run: +COMPILE +======= +On Max OS X: $gcc -o tree2PS-fast -fast avl.c tree2ps.c -Then place the binary into a directory that is listed in setup_paths.sh script located in the root of Galaxy installation +On Linux: + +Insert the following near the top of tree2ps.c + + #define isnumber (c) ( (c>='0') && (c<='9')) + +Replace -fast with -O3 (the letter 'O'): + + $gcc -o tree2PS-fast -O3 avl.c tree2ps.c + +INSTALL +======= + +Place the binary into a directory that is listed in setup_paths.sh script located in the root of Galaxy installation + + diff --git a/tools/taxonomy/T2T/t2t.readme b/tools/taxonomy/T2T/t2t.readme index fbbb82b8ed1..057408024db 100644 --- a/tools/taxonomy/T2T/t2t.readme +++ b/tools/taxonomy/T2T/t2t.readme @@ -1,6 +1,22 @@ -To compile run: +COMPILE +======= +On Max OS X: $gcc -o taxonomy2tree -fast avl.c taxonomy2tree.c -Then place the binary into a directory that is listed in setup_paths.sh script located in the root of Galaxy installation +On Linux: + +Insert the following near the top of taxonomy2tree.c + + #define isnumber (c) ( (c>='0') && (c<='9')) + +Replace -fast with -O3 (the letter 'O'): + + $gcc -o taxonomy2tree -O3 avl.c taxonomy2tree.c + +INSTALL +======= + +Place the binary into a directory that is listed in setup_paths.sh script located in the root of Galaxy installation + diff --git a/tools/taxonomy/TB/avl.c b/tools/taxonomy/TB/avl.c new file mode 100644 index 00000000000..8af208a526e --- /dev/null +++ b/tools/taxonomy/TB/avl.c @@ -0,0 +1,890 @@ +/* Produced by texiweb from libavl.w on 2002/08/24 at 13:21. */ + +/* libavl - library for manipulation of binary trees. + Copyright (C) 1998-2002 Free Software Foundation, Inc. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + See the GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA + 02111-1307, USA. + + The author may be contacted at on the Internet, or + write to Ben Pfaff, Stanford University, Computer Science Dept., 353 + Serra Mall, Stanford CA 94305, USA. +*/ + +#include +#include +#include +#include +#include "avl.h" + +/* Creates and returns a new table + with comparison function |compare| using parameter |param| + and memory allocator |allocator|. + Returns |NULL| if memory allocation failed. */ +struct avl_table * +avl_create (avl_comparison_func *compare, void *param, + struct libavl_allocator *allocator) +{ + struct avl_table *tree; + + assert (compare != NULL); + + if (allocator == NULL) + allocator = &avl_allocator_default; + + tree = allocator->libavl_malloc (allocator, sizeof *tree); + if (tree == NULL) + return NULL; + + tree->avl_root = NULL; + tree->avl_compare = compare; + tree->avl_param = param; + tree->avl_alloc = allocator; + tree->avl_count = 0; + tree->avl_generation = 0; + + return tree; +} + +/* Search |tree| for an item matching |item|, and return it if found. + Otherwise return |NULL|. */ +void * +avl_find (const struct avl_table *tree, const void *item) +{ + const struct avl_node *p; + + assert (tree != NULL && item != NULL); + for (p = tree->avl_root; p != NULL; ) + { + int cmp = tree->avl_compare (item, p->avl_data, tree->avl_param); + + if (cmp < 0) + p = p->avl_link[0]; + else if (cmp > 0) + p = p->avl_link[1]; + else /* |cmp == 0| */ + return p->avl_data; + } + + return NULL; +} + +/* Inserts |item| into |tree| and returns a pointer to |item|'s address. + If a duplicate item is found in the tree, + returns a pointer to the duplicate without inserting |item|. + Returns |NULL| in case of memory allocation failure. */ +void ** +avl_probe (struct avl_table *tree, void *item) +{ + struct avl_node *y, *z; /* Top node to update balance factor, and parent. */ + struct avl_node *p, *q; /* Iterator, and parent. */ + struct avl_node *n; /* Newly inserted node. */ + struct avl_node *w; /* New root of rebalanced subtree. */ + int dir; /* Direction to descend. */ + + unsigned char da[AVL_MAX_HEIGHT]; /* Cached comparison results. */ + int k = 0; /* Number of cached results. */ + + assert (tree != NULL && item != NULL); + + z = (struct avl_node *) &tree->avl_root; + y = tree->avl_root; + dir = 0; + for (q = z, p = y; p != NULL; q = p, p = p->avl_link[dir]) + { + int cmp = tree->avl_compare (item, p->avl_data, tree->avl_param); + if (cmp == 0) + return &p->avl_data; + + if (p->avl_balance != 0) + z = q, y = p, k = 0; + da[k++] = dir = cmp > 0; + } + + n = q->avl_link[dir] = + tree->avl_alloc->libavl_malloc (tree->avl_alloc, sizeof *n); + if (n == NULL) + return NULL; + + tree->avl_count++; + n->avl_data = item; + n->avl_link[0] = n->avl_link[1] = NULL; + n->avl_balance = 0; + if (y == NULL) + return &n->avl_data; + + for (p = y, k = 0; p != n; p = p->avl_link[da[k]], k++) + if (da[k] == 0) + p->avl_balance--; + else + p->avl_balance++; + + if (y->avl_balance == -2) + { + struct avl_node *x = y->avl_link[0]; + if (x->avl_balance == -1) + { + w = x; + y->avl_link[0] = x->avl_link[1]; + x->avl_link[1] = y; + x->avl_balance = y->avl_balance = 0; + } + else + { + assert (x->avl_balance == +1); + w = x->avl_link[1]; + x->avl_link[1] = w->avl_link[0]; + w->avl_link[0] = x; + y->avl_link[0] = w->avl_link[1]; + w->avl_link[1] = y; + if (w->avl_balance == -1) + x->avl_balance = 0, y->avl_balance = +1; + else if (w->avl_balance == 0) + x->avl_balance = y->avl_balance = 0; + else /* |w->avl_balance == +1| */ + x->avl_balance = -1, y->avl_balance = 0; + w->avl_balance = 0; + } + } + else if (y->avl_balance == +2) + { + struct avl_node *x = y->avl_link[1]; + if (x->avl_balance == +1) + { + w = x; + y->avl_link[1] = x->avl_link[0]; + x->avl_link[0] = y; + x->avl_balance = y->avl_balance = 0; + } + else + { + assert (x->avl_balance == -1); + w = x->avl_link[0]; + x->avl_link[0] = w->avl_link[1]; + w->avl_link[1] = x; + y->avl_link[1] = w->avl_link[0]; + w->avl_link[0] = y; + if (w->avl_balance == +1) + x->avl_balance = 0, y->avl_balance = -1; + else if (w->avl_balance == 0) + x->avl_balance = y->avl_balance = 0; + else /* |w->avl_balance == -1| */ + x->avl_balance = +1, y->avl_balance = 0; + w->avl_balance = 0; + } + } + else + return &n->avl_data; + z->avl_link[y != z->avl_link[0]] = w; + + tree->avl_generation++; + return &n->avl_data; +} + +/* Inserts |item| into |table|. + Returns |NULL| if |item| was successfully inserted + or if a memory allocation error occurred. + Otherwise, returns the duplicate item. */ +void * +avl_insert (struct avl_table *table, void *item) +{ + void **p = avl_probe (table, item); + return p == NULL || *p == item ? NULL : *p; +} + +/* Inserts |item| into |table|, replacing any duplicate item. + Returns |NULL| if |item| was inserted without replacing a duplicate, + or if a memory allocation error occurred. + Otherwise, returns the item that was replaced. */ +void * +avl_replace (struct avl_table *table, void *item) +{ + void **p = avl_probe (table, item); + if (p == NULL || *p == item) + return NULL; + else + { + void *r = *p; + *p = item; + return r; + } +} + +/* Deletes from |tree| and returns an item matching |item|. + Returns a null pointer if no matching item found. */ +void * +avl_delete (struct avl_table *tree, const void *item) +{ + /* Stack of nodes. */ + struct avl_node *pa[AVL_MAX_HEIGHT]; /* Nodes. */ + unsigned char da[AVL_MAX_HEIGHT]; /* |avl_link[]| indexes. */ + int k; /* Stack pointer. */ + + struct avl_node *p; /* Traverses tree to find node to delete. */ + int cmp; /* Result of comparison between |item| and |p|. */ + + assert (tree != NULL && item != NULL); + + k = 0; + p = (struct avl_node *) &tree->avl_root; + for (cmp = -1; cmp != 0; + cmp = tree->avl_compare (item, p->avl_data, tree->avl_param)) + { + int dir = cmp > 0; + + pa[k] = p; + da[k++] = dir; + + p = p->avl_link[dir]; + if (p == NULL) + return NULL; + } + item = p->avl_data; + + if (p->avl_link[1] == NULL) + pa[k - 1]->avl_link[da[k - 1]] = p->avl_link[0]; + else + { + struct avl_node *r = p->avl_link[1]; + if (r->avl_link[0] == NULL) + { + r->avl_link[0] = p->avl_link[0]; + r->avl_balance = p->avl_balance; + pa[k - 1]->avl_link[da[k - 1]] = r; + da[k] = 1; + pa[k++] = r; + } + else + { + struct avl_node *s; + int j = k++; + + for (;;) + { + da[k] = 0; + pa[k++] = r; + s = r->avl_link[0]; + if (s->avl_link[0] == NULL) + break; + + r = s; + } + + s->avl_link[0] = p->avl_link[0]; + r->avl_link[0] = s->avl_link[1]; + s->avl_link[1] = p->avl_link[1]; + s->avl_balance = p->avl_balance; + + pa[j - 1]->avl_link[da[j - 1]] = s; + da[j] = 1; + pa[j] = s; + } + } + + tree->avl_alloc->libavl_free (tree->avl_alloc, p); + + assert (k > 0); + while (--k > 0) + { + struct avl_node *y = pa[k]; + + if (da[k] == 0) + { + y->avl_balance++; + if (y->avl_balance == +1) + break; + else if (y->avl_balance == +2) + { + struct avl_node *x = y->avl_link[1]; + if (x->avl_balance == -1) + { + struct avl_node *w; + assert (x->avl_balance == -1); + w = x->avl_link[0]; + x->avl_link[0] = w->avl_link[1]; + w->avl_link[1] = x; + y->avl_link[1] = w->avl_link[0]; + w->avl_link[0] = y; + if (w->avl_balance == +1) + x->avl_balance = 0, y->avl_balance = -1; + else if (w->avl_balance == 0) + x->avl_balance = y->avl_balance = 0; + else /* |w->avl_balance == -1| */ + x->avl_balance = +1, y->avl_balance = 0; + w->avl_balance = 0; + pa[k - 1]->avl_link[da[k - 1]] = w; + } + else + { + y->avl_link[1] = x->avl_link[0]; + x->avl_link[0] = y; + pa[k - 1]->avl_link[da[k - 1]] = x; + if (x->avl_balance == 0) + { + x->avl_balance = -1; + y->avl_balance = +1; + break; + } + else + x->avl_balance = y->avl_balance = 0; + } + } + } + else + { + y->avl_balance--; + if (y->avl_balance == -1) + break; + else if (y->avl_balance == -2) + { + struct avl_node *x = y->avl_link[0]; + if (x->avl_balance == +1) + { + struct avl_node *w; + assert (x->avl_balance == +1); + w = x->avl_link[1]; + x->avl_link[1] = w->avl_link[0]; + w->avl_link[0] = x; + y->avl_link[0] = w->avl_link[1]; + w->avl_link[1] = y; + if (w->avl_balance == -1) + x->avl_balance = 0, y->avl_balance = +1; + else if (w->avl_balance == 0) + x->avl_balance = y->avl_balance = 0; + else /* |w->avl_balance == +1| */ + x->avl_balance = -1, y->avl_balance = 0; + w->avl_balance = 0; + pa[k - 1]->avl_link[da[k - 1]] = w; + } + else + { + y->avl_link[0] = x->avl_link[1]; + x->avl_link[1] = y; + pa[k - 1]->avl_link[da[k - 1]] = x; + if (x->avl_balance == 0) + { + x->avl_balance = +1; + y->avl_balance = -1; + break; + } + else + x->avl_balance = y->avl_balance = 0; + } + } + } + } + + tree->avl_count--; + tree->avl_generation++; + return (void *) item; +} + +/* Refreshes the stack of parent pointers in |trav| + and updates its generation number. */ +static void +trav_refresh (struct avl_traverser *trav) +{ + assert (trav != NULL); + + trav->avl_generation = trav->avl_table->avl_generation; + + if (trav->avl_node != NULL) + { + avl_comparison_func *cmp = trav->avl_table->avl_compare; + void *param = trav->avl_table->avl_param; + struct avl_node *node = trav->avl_node; + struct avl_node *i; + + trav->avl_height = 0; + for (i = trav->avl_table->avl_root; i != node; ) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + assert (i != NULL); + + trav->avl_stack[trav->avl_height++] = i; + i = i->avl_link[cmp (node->avl_data, i->avl_data, param) > 0]; + } + } +} + +/* Initializes |trav| for use with |tree| + and selects the null node. */ +void +avl_t_init (struct avl_traverser *trav, struct avl_table *tree) +{ + trav->avl_table = tree; + trav->avl_node = NULL; + trav->avl_height = 0; + trav->avl_generation = tree->avl_generation; +} + +/* Initializes |trav| for |tree| + and selects and returns a pointer to its least-valued item. + Returns |NULL| if |tree| contains no nodes. */ +void * +avl_t_first (struct avl_traverser *trav, struct avl_table *tree) +{ + struct avl_node *x; + + assert (tree != NULL && trav != NULL); + + trav->avl_table = tree; + trav->avl_height = 0; + trav->avl_generation = tree->avl_generation; + + x = tree->avl_root; + if (x != NULL) + while (x->avl_link[0] != NULL) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = x; + x = x->avl_link[0]; + } + trav->avl_node = x; + + return x != NULL ? x->avl_data : NULL; +} + +/* Initializes |trav| for |tree| + and selects and returns a pointer to its greatest-valued item. + Returns |NULL| if |tree| contains no nodes. */ +void * +avl_t_last (struct avl_traverser *trav, struct avl_table *tree) +{ + struct avl_node *x; + + assert (tree != NULL && trav != NULL); + + trav->avl_table = tree; + trav->avl_height = 0; + trav->avl_generation = tree->avl_generation; + + x = tree->avl_root; + if (x != NULL) + while (x->avl_link[1] != NULL) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = x; + x = x->avl_link[1]; + } + trav->avl_node = x; + + return x != NULL ? x->avl_data : NULL; +} + +/* Searches for |item| in |tree|. + If found, initializes |trav| to the item found and returns the item + as well. + If there is no matching item, initializes |trav| to the null item + and returns |NULL|. */ +void * +avl_t_find (struct avl_traverser *trav, struct avl_table *tree, void *item) +{ + struct avl_node *p, *q; + + assert (trav != NULL && tree != NULL && item != NULL); + trav->avl_table = tree; + trav->avl_height = 0; + trav->avl_generation = tree->avl_generation; + for (p = tree->avl_root; p != NULL; p = q) + { + int cmp = tree->avl_compare (item, p->avl_data, tree->avl_param); + + if (cmp < 0) + q = p->avl_link[0]; + else if (cmp > 0) + q = p->avl_link[1]; + else /* |cmp == 0| */ + { + trav->avl_node = p; + return p->avl_data; + } + + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = p; + } + + trav->avl_height = 0; + trav->avl_node = NULL; + return NULL; +} + +/* Attempts to insert |item| into |tree|. + If |item| is inserted successfully, it is returned and |trav| is + initialized to its location. + If a duplicate is found, it is returned and |trav| is initialized to + its location. No replacement of the item occurs. + If a memory allocation failure occurs, |NULL| is returned and |trav| + is initialized to the null item. */ +void * +avl_t_insert (struct avl_traverser *trav, struct avl_table *tree, void *item) +{ + void **p; + + assert (trav != NULL && tree != NULL && item != NULL); + + p = avl_probe (tree, item); + if (p != NULL) + { + trav->avl_table = tree; + trav->avl_node = + ((struct avl_node *) + ((char *) p - offsetof (struct avl_node, avl_data))); + trav->avl_generation = tree->avl_generation - 1; + return *p; + } + else + { + avl_t_init (trav, tree); + return NULL; + } +} + +/* Initializes |trav| to have the same current node as |src|. */ +void * +avl_t_copy (struct avl_traverser *trav, const struct avl_traverser *src) +{ + assert (trav != NULL && src != NULL); + + if (trav != src) + { + trav->avl_table = src->avl_table; + trav->avl_node = src->avl_node; + trav->avl_generation = src->avl_generation; + if (trav->avl_generation == trav->avl_table->avl_generation) + { + trav->avl_height = src->avl_height; + memcpy (trav->avl_stack, (const void *) src->avl_stack, + sizeof *trav->avl_stack * trav->avl_height); + } + } + + return trav->avl_node != NULL ? trav->avl_node->avl_data : NULL; +} + +/* Returns the next data item in inorder + within the tree being traversed with |trav|, + or if there are no more data items returns |NULL|. */ +void * +avl_t_next (struct avl_traverser *trav) +{ + struct avl_node *x; + + assert (trav != NULL); + + if (trav->avl_generation != trav->avl_table->avl_generation) + trav_refresh (trav); + + x = trav->avl_node; + if (x == NULL) + { + return avl_t_first (trav, trav->avl_table); + } + else if (x->avl_link[1] != NULL) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = x; + x = x->avl_link[1]; + + while (x->avl_link[0] != NULL) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = x; + x = x->avl_link[0]; + } + } + else + { + struct avl_node *y; + + do + { + if (trav->avl_height == 0) + { + trav->avl_node = NULL; + return NULL; + } + + y = x; + x = trav->avl_stack[--trav->avl_height]; + } + while (y == x->avl_link[1]); + } + trav->avl_node = x; + + return x->avl_data; +} + +/* Returns the previous data item in inorder + within the tree being traversed with |trav|, + or if there are no more data items returns |NULL|. */ +void * +avl_t_prev (struct avl_traverser *trav) +{ + struct avl_node *x; + + assert (trav != NULL); + + if (trav->avl_generation != trav->avl_table->avl_generation) + trav_refresh (trav); + + x = trav->avl_node; + if (x == NULL) + { + return avl_t_last (trav, trav->avl_table); + } + else if (x->avl_link[0] != NULL) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = x; + x = x->avl_link[0]; + + while (x->avl_link[1] != NULL) + { + assert (trav->avl_height < AVL_MAX_HEIGHT); + trav->avl_stack[trav->avl_height++] = x; + x = x->avl_link[1]; + } + } + else + { + struct avl_node *y; + + do + { + if (trav->avl_height == 0) + { + trav->avl_node = NULL; + return NULL; + } + + y = x; + x = trav->avl_stack[--trav->avl_height]; + } + while (y == x->avl_link[0]); + } + trav->avl_node = x; + + return x->avl_data; +} + +/* Returns |trav|'s current item. */ +void * +avl_t_cur (struct avl_traverser *trav) +{ + assert (trav != NULL); + + return trav->avl_node != NULL ? trav->avl_node->avl_data : NULL; +} + +/* Replaces the current item in |trav| by |new| and returns the item replaced. + |trav| must not have the null item selected. + The new item must not upset the ordering of the tree. */ +void * +avl_t_replace (struct avl_traverser *trav, void *new) +{ + void *old; + + assert (trav != NULL && trav->avl_node != NULL && new != NULL); + old = trav->avl_node->avl_data; + trav->avl_node->avl_data = new; + return old; +} + +static void +copy_error_recovery (struct avl_node **stack, int height, + struct avl_table *new, avl_item_func *destroy) +{ + assert (stack != NULL && height >= 0 && new != NULL); + + for (; height > 2; height -= 2) + stack[height - 1]->avl_link[1] = NULL; + avl_destroy (new, destroy); +} + +/* Copies |org| to a newly created tree, which is returned. + If |copy != NULL|, each data item in |org| is first passed to |copy|, + and the return values are inserted into the tree, + with |NULL| return values taken as indications of failure. + On failure, destroys the partially created new tree, + applying |destroy|, if non-null, to each item in the new tree so far, + and returns |NULL|. + If |allocator != NULL|, it is used for allocation in the new tree. + Otherwise, the same allocator used for |org| is used. */ +struct avl_table * +avl_copy (const struct avl_table *org, avl_copy_func *copy, + avl_item_func *destroy, struct libavl_allocator *allocator) +{ + struct avl_node *stack[2 * (AVL_MAX_HEIGHT + 1)]; + int height = 0; + + struct avl_table *new; + const struct avl_node *x; + struct avl_node *y; + + assert (org != NULL); + new = avl_create (org->avl_compare, org->avl_param, + allocator != NULL ? allocator : org->avl_alloc); + if (new == NULL) + return NULL; + new->avl_count = org->avl_count; + if (new->avl_count == 0) + return new; + + x = (const struct avl_node *) &org->avl_root; + y = (struct avl_node *) &new->avl_root; + for (;;) + { + while (x->avl_link[0] != NULL) + { + assert (height < 2 * (AVL_MAX_HEIGHT + 1)); + + y->avl_link[0] = + new->avl_alloc->libavl_malloc (new->avl_alloc, + sizeof *y->avl_link[0]); + if (y->avl_link[0] == NULL) + { + if (y != (struct avl_node *) &new->avl_root) + { + y->avl_data = NULL; + y->avl_link[1] = NULL; + } + + copy_error_recovery (stack, height, new, destroy); + return NULL; + } + + stack[height++] = (struct avl_node *) x; + stack[height++] = y; + x = x->avl_link[0]; + y = y->avl_link[0]; + } + y->avl_link[0] = NULL; + + for (;;) + { + y->avl_balance = x->avl_balance; + if (copy == NULL) + y->avl_data = x->avl_data; + else + { + y->avl_data = copy (x->avl_data, org->avl_param); + if (y->avl_data == NULL) + { + y->avl_link[1] = NULL; + copy_error_recovery (stack, height, new, destroy); + return NULL; + } + } + + if (x->avl_link[1] != NULL) + { + y->avl_link[1] = + new->avl_alloc->libavl_malloc (new->avl_alloc, + sizeof *y->avl_link[1]); + if (y->avl_link[1] == NULL) + { + copy_error_recovery (stack, height, new, destroy); + return NULL; + } + + x = x->avl_link[1]; + y = y->avl_link[1]; + break; + } + else + y->avl_link[1] = NULL; + + if (height <= 2) + return new; + + y = stack[--height]; + x = stack[--height]; + } + } +} + +/* Frees storage allocated for |tree|. + If |destroy != NULL|, applies it to each data item in inorder. */ +void +avl_destroy (struct avl_table *tree, avl_item_func *destroy) +{ + struct avl_node *p, *q; + + assert (tree != NULL); + + for (p = tree->avl_root; p != NULL; p = q) + if (p->avl_link[0] == NULL) + { + q = p->avl_link[1]; + if (destroy != NULL && p->avl_data != NULL) + destroy (p->avl_data, tree->avl_param); + tree->avl_alloc->libavl_free (tree->avl_alloc, p); + } + else + { + q = p->avl_link[0]; + p->avl_link[0] = q->avl_link[1]; + q->avl_link[1] = p; + } + + tree->avl_alloc->libavl_free (tree->avl_alloc, tree); +} + +/* Allocates |size| bytes of space using |malloc()|. + Returns a null pointer if allocation fails. */ +void * +avl_malloc (struct libavl_allocator *allocator, size_t size) +{ + assert (allocator != NULL && size > 0); + return malloc (size); +} + +/* Frees |block|. */ +void +avl_free (struct libavl_allocator *allocator, void *block) +{ + assert (allocator != NULL && block != NULL); + free (block); +} + +/* Default memory allocator that uses |malloc()| and |free()|. */ +struct libavl_allocator avl_allocator_default = + { + avl_malloc, + avl_free + }; + +#undef NDEBUG +#include + +/* Asserts that |avl_insert()| succeeds at inserting |item| into |table|. */ +void +(avl_assert_insert) (struct avl_table *table, void *item) +{ + void **p = avl_probe (table, item); + assert (p != NULL && *p == item); +} + +/* Asserts that |avl_delete()| really removes |item| from |table|, + and returns the removed item. */ +void * +(avl_assert_delete) (struct avl_table *table, void *item) +{ + void *p = avl_delete (table, item); + assert (p != NULL); + return p; +} + diff --git a/tools/taxonomy/TB/avl.h b/tools/taxonomy/TB/avl.h new file mode 100644 index 00000000000..a9cf3c1b6d0 --- /dev/null +++ b/tools/taxonomy/TB/avl.h @@ -0,0 +1,115 @@ +/* Produced by texiweb from libavl.w on 2002/08/24 at 13:21. */ + +/* libavl - library for manipulation of binary trees. + Copyright (C) 1998-2002 Free Software Foundation, Inc. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + See the GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA + 02111-1307, USA. + + The author may be contacted at on the Internet, or + write to Ben Pfaff, Stanford University, Computer Science Dept., 353 + Serra Mall, Stanford CA 94305, USA. +*/ + +#ifndef AVL_H +#define AVL_H 1 + +#include + +/* Function types. */ +typedef int avl_comparison_func (const void *avl_a, const void *avl_b, + void *avl_param); +typedef void avl_item_func (void *avl_item, void *avl_param); +typedef void *avl_copy_func (void *avl_item, void *avl_param); + +#ifndef LIBAVL_ALLOCATOR +#define LIBAVL_ALLOCATOR +/* Memory allocator. */ +struct libavl_allocator + { + void *(*libavl_malloc) (struct libavl_allocator *, size_t libavl_size); + void (*libavl_free) (struct libavl_allocator *, void *libavl_block); + }; +#endif + +/* Default memory allocator. */ +extern struct libavl_allocator avl_allocator_default; +void *avl_malloc (struct libavl_allocator *, size_t); +void avl_free (struct libavl_allocator *, void *); + +/* Maximum AVL height. */ +#ifndef AVL_MAX_HEIGHT +#define AVL_MAX_HEIGHT 32 +#endif + +/* Tree data structure. */ +struct avl_table + { + struct avl_node *avl_root; /* Tree's root. */ + avl_comparison_func *avl_compare; /* Comparison function. */ + void *avl_param; /* Extra argument to |avl_compare|. */ + struct libavl_allocator *avl_alloc; /* Memory allocator. */ + size_t avl_count; /* Number of items in tree. */ + unsigned long avl_generation; /* Generation number. */ + }; + +/* An AVL tree node. */ +struct avl_node + { + struct avl_node *avl_link[2]; /* Subtrees. */ + void *avl_data; /* Pointer to data. */ + signed char avl_balance; /* Balance factor. */ + }; + +/* AVL traverser structure. */ +struct avl_traverser + { + struct avl_table *avl_table; /* Tree being traversed. */ + struct avl_node *avl_node; /* Current node in tree. */ + struct avl_node *avl_stack[AVL_MAX_HEIGHT]; + /* All the nodes above |avl_node|. */ + size_t avl_height; /* Number of nodes in |avl_parent|. */ + unsigned long avl_generation; /* Generation number. */ + }; + +/* Table functions. */ +struct avl_table *avl_create (avl_comparison_func *, void *, + struct libavl_allocator *); +struct avl_table *avl_copy (const struct avl_table *, avl_copy_func *, + avl_item_func *, struct libavl_allocator *); +void avl_destroy (struct avl_table *, avl_item_func *); +void **avl_probe (struct avl_table *, void *); +void *avl_insert (struct avl_table *, void *); +void *avl_replace (struct avl_table *, void *); +void *avl_delete (struct avl_table *, const void *); +void *avl_find (const struct avl_table *, const void *); +void avl_assert_insert (struct avl_table *, void *); +void *avl_assert_delete (struct avl_table *, void *); + +#define avl_count(table) ((size_t) (table)->avl_count) + +/* Table traverser functions. */ +void avl_t_init (struct avl_traverser *, struct avl_table *); +void *avl_t_first (struct avl_traverser *, struct avl_table *); +void *avl_t_last (struct avl_traverser *, struct avl_table *); +void *avl_t_find (struct avl_traverser *, struct avl_table *, void *); +void *avl_t_insert (struct avl_traverser *, struct avl_table *, void *); +void *avl_t_copy (struct avl_traverser *, const struct avl_traverser *); +void *avl_t_next (struct avl_traverser *); +void *avl_t_prev (struct avl_traverser *); +void *avl_t_cur (struct avl_traverser *); +void *avl_t_replace (struct avl_traverser *, void *); + +#endif /* avl.h */ diff --git a/tools/taxonomy/TB/main.c b/tools/taxonomy/TB/main.c new file mode 100644 index 00000000000..8a62e3938c3 --- /dev/null +++ b/tools/taxonomy/TB/main.c @@ -0,0 +1,788 @@ +#include +#include +#include +#include +#include +#include + +#include "avl.h" + +#define NUMBER_OF_TAX_FIELDS 22 + +char *rankLabels [NUMBER_OF_TAX_FIELDS] = + {"root" , + "superkingdom", + "kingdom" , + "subkingdom" , + "superphylum" , + "phylum" , + "subphylum" , + "superclass" , + "class" , + "subclass" , + "superorder" , + "order" , + "suborder" , + "superfamily" , + "family" , + "subfamily" , + "tribe" , + "subtribe" , + "genus" , + "subgenus" , + "species" , + "subspecies" +}; + +char noValue[] = "n\t", + rootValue[] = "toor\t", + usage[] = "Incorrect number of arguments.\nExpected arguments: tax_id->label file, tax_id->hierarchy file, input file (gid taxid pairs), output file (24 columns)"; + +/*#define DEBUG_ME*/ + +/*---------------------------------------------------------------------------------------------------- */ + +void check_pointer (void*); + +/*---------------------------------------------------------------------------------------------------- */ + +struct avl_table * nameTagAVL = NULL, + * rankLabelAVL = NULL; + +#define DEFAULT_STRING_ALLOC 16L + +/*---------------------------------------------------------------------------------------------------- */ + +struct bufferedString +{ + char *sData; + long sLength, + saLength; +} +*globalNameBuffer; + + +/*---------------------------------------------------------------------------------------------------- */ + +struct bufferedString *allocateNewString (void) +{ + struct bufferedString *newS = (struct bufferedString*)malloc (sizeof (struct bufferedString)); + check_pointer (newS); + check_pointer (newS->sData = (char*)malloc (DEFAULT_STRING_ALLOC+1)); + newS->sLength = 0; + newS->saLength = DEFAULT_STRING_ALLOC; + newS->sData[0] = 0; + return newS; +} + +/*---------------------------------------------------------------------------------------------------- */ + +void clear_buffered_string (struct bufferedString* theString) +{ + theString->sLength = 0; +} + +/*---------------------------------------------------------------------------------------------------- */ + +void appendCharacterToString (struct bufferedString * s, const char c) +{ + long addThis; + if (s->saLength == s->sLength) + { + addThis = s->saLength / 8; + if (DEFAULT_STRING_ALLOC > addThis) + addThis = DEFAULT_STRING_ALLOC; + s->saLength += addThis; + check_pointer (s->sData = realloc (s->sData,s->saLength+1)); + } + s->sData[s->sLength] = c; + s->sData[++s->sLength] = 0; +} + +/*---------------------------------------------------------------------------------------------------- */ + +long appendRangeToString (struct bufferedString * d, struct bufferedString *s, long from, long to) +{ + long addThis, + pl = to-from+1; + + if (pl<=0) + return -1; + + if (d->saLength < d->sLength + pl) + { + addThis = d->saLength / 8; + + if (DEFAULT_STRING_ALLOC > addThis) + addThis = DEFAULT_STRING_ALLOC; + if (addThis < pl) + addThis = pl; + + d->saLength += addThis; + check_pointer (d->sData = realloc (d->sData,d->saLength+1)); + } + for (addThis = from; addThis <=to; addThis++) + d->sData[d->sLength++] = s->sData[addThis]; + + d->sData[d->sLength] = 0; + + return pl; +} + +/*---------------------------------------------------------------------------------------------------- */ + +void appendCharRangeToString (struct bufferedString * d, char * buffer) +{ + long addThis, + pl = strlen(buffer); + + if (pl<=0) + return; + + if (d->saLength < d->sLength + pl) + { + addThis = d->saLength / 8; + + if (DEFAULT_STRING_ALLOC > addThis) + addThis = DEFAULT_STRING_ALLOC; + if (addThis < pl) + addThis = pl; + + d->saLength += addThis; + check_pointer (d->sData = realloc (d->sData,d->saLength+1)); + } + for (addThis = 0; addThis sData[d->sLength++] = buffer[addThis]; + + d->sData[d->sLength] = 0; + +} + +/*---------------------------------------------------------------------------------------------------- */ + +long appendCharBufferToString (struct bufferedString * d, const char * b) +{ + long addThis, + pl = strlen (b); + + if (pl<=0) + return -1; + + if (d->saLength < d->sLength + pl) + { + addThis = d->saLength / 8; + + if (DEFAULT_STRING_ALLOC > addThis) + addThis = DEFAULT_STRING_ALLOC; + if (addThis < pl) + addThis = pl; + + d->saLength += addThis; + check_pointer (d->sData = realloc (d->sData,d->saLength+1)); + } + for (addThis = 0; addThis sData[d->sLength++] = b[addThis]; + + d->sData[d->sLength] = 0; + return pl; +} + + +/*---------------------------------------------------------------------------------------------------- */ +/*---------------------------------------------------------------------------------------------------- */ + +struct storedNameTag +{ + long taxID, + startIndex, + length; + + char taxonomy_level; + + struct storedNameTag * parent; +}; + +/*---------------------------------------------------------------------------------------------------- */ + +struct storedNameTag * allocateNameTag (void) +{ + struct storedNameTag *newS = (struct storedNameTag*)malloc (sizeof (struct storedNameTag)); + check_pointer (newS); + newS->taxID = 0; + newS->startIndex = 0; + newS->length = 0; + newS->parent = NULL; + newS->taxonomy_level = -1; + return newS; +} + +/*---------------------------------------------------------------------------------------------------- */ + +int compare_tags (const void *avl_a, const void *avl_b, void * xtra) +{ + long n1 = ((struct storedNameTag*)avl_a)->taxID; + long n2 = ((struct storedNameTag*)avl_b)->taxID; + + if (n1 > n2) return 1; + if (n1 < n2) return -1; + return 0; +} + +/*---------------------------------------------------------------------------------------------------- */ + +struct bufferedString * nameByID (long taxID) +{ + static struct storedNameTag queryTag; + struct storedNameTag * res; + struct bufferedString * resStr = NULL; + + queryTag.taxID = taxID; + res = (struct storedNameTag *)avl_find(nameTagAVL, &queryTag); + if (res) + { + resStr = allocateNewString(); + appendRangeToString(resStr,globalNameBuffer,res->startIndex,res->startIndex+res->length-1); + } + return resStr; +} + +/*---------------------------------------------------------------------------------------------------- */ + +struct storedNameTag * tagByID (long taxID) +{ + static struct storedNameTag queryTag; + queryTag.taxID = taxID; + return (struct storedNameTag *)avl_find(nameTagAVL, &queryTag); +} + +/*---------------------------------------------------------------------------------------------------- */ + +struct bufferedString * walkPath (long taxID) +{ + struct storedNameTag * currentTag = tagByID(taxID); + if (!currentTag) + return NULL; + + struct bufferedString * resStr = allocateNewString(); + long k, + level = NUMBER_OF_TAX_FIELDS-1; + + char // buffer [32], + first_record = 1; + + while (currentTag) { + + if (currentTag->taxonomy_level >= 0) + { + if (first_record) + { + for (k = currentTag->taxonomy_level+1; k < NUMBER_OF_TAX_FIELDS; k++) + appendCharBufferToString(resStr,noValue); + first_record = 0; + } + else + { + for (k=currentTag->taxonomy_level+1; ktaxonomy_level); + //for (k=strlen(buffer)-1; k>=0; k--) + // appendCharacterToString(resStr,buffer[k]); + //appendCharacterToString(resStr,'('); + for (k=currentTag->startIndex+currentTag->length-1; k>=currentTag->startIndex; k--) + appendCharacterToString(resStr,globalNameBuffer->sData[k]); + + appendCharacterToString(resStr,'\t'); + level = currentTag->taxonomy_level; + } + currentTag = currentTag->parent; + } + + for (k = level-1; k ; k--) + appendCharBufferToString(resStr,noValue); + + appendCharBufferToString(resStr,rootValue); + + return resStr; +} + +/*---------------------------------------------------------------------------------------------------- */ +/*---------------------------------------------------------------------------------------------------- */ + +void check_pointer (void * p) +{ + if (p == NULL) + { + fprintf (stderr,"Memory allocation error\n"); + exit (1); + } +} + +/*---------------------------------------------------------------------------------------------------- */ + +int compare_strings (struct bufferedString * s1, struct bufferedString * s2) +{ + long upTo, + i; + + if (s1->sLength>s2->sLength) + upTo = s2->sLength; + else + upTo = s1->sLength; + + for (i=0; isData[i]-s2->sData[i]); + if (res < 0) + return -1; + else + if (res>0) + return 1; + } + + if (s1->sLength == s2->sLength) + return 0; + + return 1-2*(s1->sLengthsLength); +} + +/*---------------------------------------------------------------------------------------------------- */ + +int compare_tag_strings (const void *avl_a, const void *avl_b, void * xtra) +{ + struct storedNameTag* s1 = (struct storedNameTag*)avl_a; + struct storedNameTag* s2 = (struct storedNameTag*)avl_b; + char *buffa = ((struct bufferedString*)xtra)->sData; + + long upTo, + i; + + if (s1->length>s2->length) + upTo = s2->length; + else + upTo = s1->length; + + for (i=0; istartIndex+i]-buffa[s2->startIndex+i]; + if (res < 0) + return -1; + else + if (res>0) + return 1; + } + + if (s1->length == s2->length) + return 0; + + return 1-2*(s1->lengthlength); +} + + + +/*---------------------------------------------------------------------------------------------------- */ + +void destroy_string (struct bufferedString* aStr) +{ + free (aStr->sData); + free (aStr); +} + + +/*---------------------------------------------------------------------------------------------------- */ + +void reportError (char * theMessage) +{ + fprintf (stderr, "\nERROR: %s\n", theMessage); + exit (1); +} + +/*---------------------------------------------------------------------------------------------------- */ + +void reportErrorLine (char * theMessage, long lineID) +{ + fprintf (stderr, "\nERROR in line %d: %s\n", lineID, theMessage); + exit (1); +} + + +/*---------------------------------------------------------------------------------------------------- */ + +int main (int argc, const char * argv[]) +{ + FILE *inFile, + *structFile, + *queryFile, + *outFile; + + struct bufferedString *scientificName = allocateNewString(), + **currentBuffers, + *qry; + + struct storedNameTag *aTag, + *aTag2, + *aTag3 = allocateNameTag(); + + char automatonState = 0, + currentField = 0, + currentChar = 0, + taxonBuffer [1024], + firstLine; + + long currentLineID = 1, + expectedFields = 4, + indexer, + indexer2; + + + + if (argc != 5) + { + fprintf (stderr,"%s\n", usage); + return 1; + } + + appendCharBufferToString(scientificName, "scientific name"); + globalNameBuffer = allocateNewString(); + currentBuffers = (struct bufferedString**)malloc (expectedFields*sizeof (struct bufferedString*)); + nameTagAVL = avl_create (compare_tags, NULL, NULL); + rankLabelAVL = avl_create (compare_tag_strings, globalNameBuffer, NULL); + + for (indexer = 0; indexer < expectedFields; indexer++) + currentBuffers[indexer] = allocateNewString(); + + for (indexer = 0; indexer < NUMBER_OF_TAX_FIELDS; indexer++) + { + aTag = allocateNameTag (); + aTag->startIndex = globalNameBuffer->sLength; + appendCharBufferToString (globalNameBuffer,rankLabels[indexer]); + aTag->taxonomy_level = indexer; + aTag->length = strlen(rankLabels[indexer]); + if ((aTag2 = *avl_probe(rankLabelAVL, aTag)) != aTag) + reportErrorLine ("Duplicate taxonomic rank name", aTag2->taxonomy_level); + } + + + + inFile = fopen (argv[1], "rb"); + + if (!inFile) + { + fprintf (stderr,"Failed to open input file: %s\n", argv[1]); + return 1; + } + + structFile = fopen (argv[2], "rb"); + + if (!structFile) + { + fprintf (stderr,"Failed to open input file: %s\n", argv[2]); + return 1; + } + + queryFile = fopen (argv[3], "rb"); + + if (!queryFile) + { + fprintf (stderr,"Failed to open input file: %s\n", argv[3]); + return 1; + } + + outFile = fopen (argv[4], "w"); + + if (!outFile) + { + fprintf (stderr,"Failed to open output file: %s\n", argv[3]); + return 1; + } + + currentChar = fgetc(inFile); + while (!feof(inFile)) + { + switch (automatonState) + { + case 0: /* start of the line; expecting numbers */ + if (currentChar >= '0' && currentChar <='9') + { + automatonState = 1; /* reading sequence ID */ + appendCharacterToString(currentBuffers[currentField],currentChar); + } + else + if (!(currentChar == '\n' || currentChar == '\r')) + reportErrorLine ("Could not find a valid sequence ID to start the line",currentLineID); + break; + + case 1: /* reading sequence ID */ + if (currentChar >= '0' && currentChar <='9') + appendCharacterToString(currentBuffers[currentField],currentChar); + else + if (currentChar == '\t') + automatonState = 2; + else + reportErrorLine ("Expected a tab following the tax ID",currentLineID); + break; + + case 2: /* looking for a | */ + if (currentChar != '|') + reportErrorLine ("Expected a '|' following the tab",currentLineID); + else + automatonState = 3; + break; + + case 3: /* looking for a \t or a \n|\r*/ + if (currentChar == '\t') + { + automatonState = 4; + currentField ++; + if (currentField == expectedFields) + reportErrorLine ("Too many fields",currentLineID); + } + else + if (currentChar == '\n' || currentChar == '\r') + { + if (currentField < expectedFields-1) + reportErrorLine ("Too few fields",currentLineID); + automatonState = 0; + currentLineID ++; + currentField = 0; + if (compare_strings(currentBuffers[3],scientificName)==0) + { + aTag = allocateNameTag(); + aTag->taxID = atoi(currentBuffers[0]->sData); + aTag->startIndex = globalNameBuffer->sLength; + aTag->length = appendRangeToString(globalNameBuffer,currentBuffers[1],0,currentBuffers[1]->sLength-1); + if (aTag->length <= 0) + reportErrorLine ("Empty name tag",currentLineID); + if (*avl_probe(nameTagAVL,aTag) != aTag) + reportErrorLine ("Duplicate name tag",currentLineID); + + } + for (indexer = 0; indexer < expectedFields; indexer++) + clear_buffered_string(currentBuffers[indexer]); + } + else + reportErrorLine ("Expected a tab following the '|'",currentLineID); + break; + + case 4: /* read a field */ + if (currentChar == '\t') + automatonState = 2; + else + if (currentChar == '\n' || currentChar == '\r') + reportErrorLine ("Unexpected end-of-line",currentLineID); + else + appendCharacterToString(currentBuffers[currentField],currentChar); + break; + + } + currentChar = fgetc(inFile); + } + + fclose (inFile); + + for (indexer = 0; indexer < expectedFields; indexer++) + destroy_string(currentBuffers[indexer]); + + free(currentBuffers); + expectedFields = 13; + currentBuffers = (struct bufferedString**)malloc (expectedFields*sizeof (struct bufferedString*)); + + for (indexer = 0; indexer < expectedFields; indexer++) + currentBuffers[indexer] = allocateNewString(); + + automatonState = 0; + currentLineID = 1; + currentField = 0; + + currentChar = fgetc(structFile); + + while (!feof(structFile)) + { + switch (automatonState) + { + case 0: /* start of the line; expecting numbers */ + if (currentChar >= '0' && currentChar <='9') + { + automatonState = 1; /* reading sequence ID */ + appendCharacterToString(currentBuffers[currentField],currentChar); + } + else + if (!(currentChar == '\n' || currentChar == '\r')) + reportErrorLine ("Could not find a valid sequence ID to start the line",currentLineID); + break; + + case 1: /* reading sequence ID */ + if (currentChar >= '0' && currentChar <='9') + appendCharacterToString(currentBuffers[currentField],currentChar); + else + if (currentChar == '\t') + automatonState = 2; + else + reportErrorLine ("Expected a tab following the tax ID",currentLineID); + break; + + case 2: /* looking for a | */ + if (currentChar != '|') + reportErrorLine ("Expected a '|' following the tab",currentLineID); + else + automatonState = 3; + break; + + case 3: /* looking for a \t or a \n|\r*/ + if (currentChar == '\t') + { + automatonState = 4; + currentField ++; + if (currentField == expectedFields) + reportErrorLine ("Too many fields",currentLineID); + } + else + if (currentChar == '\n' || currentChar == '\r') + { + if (currentField < expectedFields-1) + reportErrorLine ("Too few fields",currentLineID); + + aTag = tagByID(atoi(currentBuffers[0]->sData)); + aTag2 = tagByID(atoi(currentBuffers[1]->sData)); + + if (! (aTag && aTag2)) + reportErrorLine ("Invalid ID tag",currentLineID); + + if (aTag2 != aTag) + { + aTag->parent = aTag2; + aTag3->startIndex = globalNameBuffer->sLength; + aTag3->length = currentBuffers[2]->sLength; + appendRangeToString (globalNameBuffer, currentBuffers[2], 0, currentBuffers[2]->sLength-1); + aTag2 = avl_find (rankLabelAVL, aTag3); + if (aTag2) + aTag->taxonomy_level = aTag2->taxonomy_level; + globalNameBuffer->sLength = aTag3->startIndex; + } + + currentField = 0; + automatonState = 0; + currentLineID ++; + + for (indexer = 0; indexer < expectedFields; indexer++) + clear_buffered_string(currentBuffers[indexer]); + } + else + reportErrorLine ("Expected a tab following the '|'",currentLineID); + break; + + case 4: /* read a field */ + if (currentChar == '\t') + automatonState = 2; + else + if (currentChar == '\n' || currentChar == '\r') + reportErrorLine ("Unexpected end-of-line",currentLineID); + else + appendCharacterToString(currentBuffers[currentField],currentChar); + break; + + } + currentChar = fgetc(structFile); + } + + fclose (structFile); + + for (indexer = 0; indexer < expectedFields; indexer++) + destroy_string(currentBuffers[indexer]); + free(currentBuffers); + expectedFields = 2; + currentBuffers = (struct bufferedString**)malloc (expectedFields*sizeof (struct bufferedString*)); + + for (indexer = 0; indexer <= expectedFields; indexer++) + currentBuffers[indexer] = allocateNewString(); + + automatonState = 0; + currentLineID = 1; + currentField = 0; + firstLine = 1; + + currentChar = fgetc(queryFile); + + while (!feof(queryFile)) + { + switch (automatonState) + { + case 0: /* start of the line; expecting numbers */ + if (!isspace(currentChar)) + { + automatonState = 1; /* reading sequence ID */ + appendCharacterToString(currentBuffers[currentField],currentChar); + } + else + if (!(currentChar == '\n' || currentChar == '\r')) + reportErrorLine ("Could not find a valid sequence ID to start the line",currentLineID); + break; + + case 1: /* reading sequence ID */ + if ((currentField == 1 && currentChar >= '0' && currentChar <='9')||(currentField == 0 && currentChar != '\n' && currentChar != '\r' && currentChar != '\t')) + appendCharacterToString(currentBuffers[currentField],currentChar); + else + if (currentChar == '\t') + { + automatonState = 1; + currentField ++; + if (currentField == expectedFields) + automatonState = 2; + } + else + if (currentChar == '\n' || currentChar == '\r') + { + currentField++; + automatonState = 2; + ungetc (currentChar,queryFile); + } + else + reportErrorLine ("Expected a tab following a Tax/GID",currentLineID); + break; + + case 2: + if (currentChar == '\n' || currentChar == '\r') + { + if (currentField == expectedFields) + { + + indexer2 = atoi (currentBuffers[1]->sData); // TaxID + qry = walkPath (indexer2); + if (firstLine) + firstLine = 0; + else + fputc ('\n',outFile); + + fprintf (outFile,"%s\t%d", currentBuffers[0]->sData, indexer2); + if (qry) + { + for (indexer = qry->sLength-1; indexer >= 0; indexer--) + fputc (qry->sData[indexer], outFile); + destroy_string(qry); + } + else + fprintf (outFile, "\tNULL"); + if (currentBuffers[expectedFields]->sLength) + fprintf (outFile, "\t%s", currentBuffers[expectedFields]->sData); + + } + currentLineID ++; + currentField = 0; + for (indexer = 0; indexer <= expectedFields; indexer++) + clear_buffered_string(currentBuffers[indexer]); + automatonState = 0; + } + else + appendCharacterToString(currentBuffers[currentField],currentChar); + + + } + currentChar = fgetc(queryFile); + } + + fclose (structFile); + fclose (queryFile); + fclose (outFile); + return 0; +} diff --git a/tools/taxonomy/TB/tb.readme b/tools/taxonomy/TB/tb.readme new file mode 100644 index 00000000000..6848a96f1f1 --- /dev/null +++ b/tools/taxonomy/TB/tb.readme @@ -0,0 +1,21 @@ +COMPILE +======= + +On Max OS X: +$gcc -o taxBuilder -fast avl.c main.c + +On Linux: + +Insert the following near the top of main.c + + #define isnumber (c) ( (c>='0') && (c<='9')) + +Replace -fast with -O3 (the letter 'O'): + + $gcc -o taxBuilder -O3 avl.c main.c + +INSTALL +======= + +Place the binary into a directory that is listed in setup_paths.sh script located in the root of Galaxy installation + diff --git a/tools/taxonomy/find_diag_hits.py b/tools/taxonomy/find_diag_hits.py index dee1da580aa..dc31297d1c9 100644 --- a/tools/taxonomy/find_diag_hits.py +++ b/tools/taxonomy/find_diag_hits.py @@ -97,12 +97,14 @@ try: elif sys.argv[4] == 'counts': out_format = False else: - stop_err('Plase specify "reads" or "counts" for output format\n') + stop_err('Please specify "reads" or "counts" for output format\n') out_file = open(sys.argv[5], 'w') except: stop_err('Check arguments\n') +if taxa[0] == 'None': stop_err('Please, use checkboxes to specify taxonomic ranks.\n') + sql = "" for i in range(len(taxa)): if taxa[i] == 'order': taxa[i] = 'ord' # SQL does not like fields to be named 'order' @@ -133,7 +135,7 @@ try: val_string = "insert into tax values(" + val_string + ")" cur.execute(val_string) except Exception, e: - stop_err(e) + stop_err('%s\n' % e) tax_file.close() diff --git a/tools/taxonomy/find_diag_hits.xml b/tools/taxonomy/find_diag_hits.xml index 6de4309cf5c..5b528d41719 100644 --- a/tools/taxonomy/find_diag_hits.xml +++ b/tools/taxonomy/find_diag_hits.xml @@ -24,7 +24,7 @@ - + diff --git a/tools/taxonomy/gi2taxonomy.py b/tools/taxonomy/gi2taxonomy.py new file mode 100644 index 00000000000..a3eea5af29a --- /dev/null +++ b/tools/taxonomy/gi2taxonomy.py @@ -0,0 +1,172 @@ +import sys +import string +import tempfile +import subprocess + +# ----------------------------------------------------------------------------------- + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit() + +# ----------------------------------------------------------------------------------- +def gi_name_to_sorted_list(file_name, gi_col, name_col): + """ Suppose input file looks like this: + a 2 + b 4 + c 5 + d 5 + where column 1 is gi_col and column 0 is name_col + output of this function will look like this: + [[2, 'a'], [4, 'b'], [5, 'c'], [5, 'd']] + """ + + result = [] + try: + F = open( file_name, 'r' ) + try: + for line in F: + file_cols = string.split(line.rstrip(), '\t') + file_cols[gi_col] = int( file_cols[gi_col] ) + result.append( [ file_cols[gi_col], file_cols[name_col] ] ) + except: + print >>sys.stderr, 'Non numeric GI field...skipping' + + except Exception, e: + stop_err('%s\n' % e) + F.close() + result.sort() + return result + +# ----------------------------------------------------------------------------------- + +def collapse_repeating_gis( L ): + """ Accepts 2-d array of gi-key pairs such as this + L = [ + [gi1, 'key1'], + [gi1, 'key2'], + [gi2','key3'] + ] + + Returns this: + [ [gi1, 'key1', 'key2'], + [gi2, 'key3' ] + ] + + The first value in each sublist MUST be int + """ + gi = [] + i = 0 + result = [] + + try: + for item in L: + if i == 0: + prev = item[0] + + if prev != item[0]: + prev_L = [] + prev_L.append( prev ) + result.append( prev_L + gi ) + prev = item[0] + gi =[] + + gi.append( item[1] ) + i += 1 + + except Exception, e: + stop_err('%s\n' % e) + + prev_L = [] + prev_L.append( prev ) + result.append( prev_L + gi ) + del(L) + return result + +# ----------------------------------------------------------------------------------- + +def get_taxId( gi2tax_file, gi_name_list, out_file ): + """ Maps GI numbers from gi_name_list to TaxId identifiers from gi2tax_file and + prints result to out_file + + gi2tax_file MUST be sorted on GI column + + gi_name_list is a list that look slike this: + [[1,'a'], [2,'b','x'], [7,'c'], [10,'d'], [90,'f']] + where the first element of each sublist is a GI number + this list MUST also be sorted on GI + + This function searches through 117,000,000 rows of gi2taxId file from NCBI + in approximately 4 minutes. This time is not dependent on the length of + gi_name_list + """ + + L = gi_name_list.pop(0) + my_gi = L[0] + F = open( out_file, 'w' ) + gi = 0 + for line in file( gi2tax_file ): + line = line.rstrip() + gi, taxId = string.split( line, '\t' ) + gi = int( gi ) + + if gi > my_gi: + try: + while ( my_gi < gi ): + L = gi_name_list.pop(0) + my_gi = L[0] + except: + break + + if gi == my_gi: + for i in range( 1,len( L ) ): + print >>F, '%s\t%s\t%d' % (L[i], taxId, gi) + try: + L = gi_name_list.pop(0) + my_gi = L[0] + except: + break + +# ----------------------------------------------------------------------------------- + + +try: + in_f = sys.argv[1] # input file with GIs + gi_col = int( sys.argv[2] ) - 1 # column in input containing GIs + name_col = int( sys.argv[3] ) - 1 # column containing sequence names + out_f = sys.argv[4] # output file +except: + stor_err('Check arguments\n') + +# GI2TAX point to a file produced by concatenation of: +# ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid_nucl.zip +# and +# ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid_prot.zip +# a sorting using this command: +# sort -n -k 1 + +GI2TAX = '/Volumes/METAG/processTaxTest/gi_taxId_sorted.txt' + +# NAME_FILE and NODE_FILE point to names.dmg and nodes.dmg +# files contained within: +# ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz + +NAME_FILE = '/Volumes/METAG/processTaxTest/names.dmp' +NODE_FILE = '/Volumes/METAG/processTaxTest/nodes.dmp' + +g2n = gi_name_to_sorted_list(in_f, gi_col, name_col) + +if len(g2n) == 0: + stop_err('No valid GI-containing fields. Please, check your column assignments.\n') + +tb_F = tempfile.NamedTemporaryFile('w') + +get_taxId( GI2TAX, collapse_repeating_gis( g2n ), tb_F.name ) + +try: + tb_cmd = 'taxBuilder %s %s %s %s' % ( NAME_FILE, NODE_FILE, tb_F.name, out_f ) + retcode = subprocess.call( tb_cmd, shell=True ) + if retcode < 0: + print >>sys.stderr, "Execution of taxBuilder terminated by signal", -retcode +except OSError, e: + print >>sys.stderr, "Execution of taxBuilder2tree failed:", e diff --git a/tools/taxonomy/gi2taxonomy.xml b/tools/taxonomy/gi2taxonomy.xml index ca26eaa437f..0f14b3feff7 100644 --- a/tools/taxonomy/gi2taxonomy.xml +++ b/tools/taxonomy/gi2taxonomy.xml @@ -1,9 +1,9 @@ - + for a list of GIs - tax.py $input $out_file1 $idField $giField + gi2taxonomy.py $input $giField $idField $out_file1 - + @@ -49,10 +49,10 @@ and you want to obtain full taxonomic representation for GIs listed in *targetGI the tool will generate the following output (you may need to scroll sideways to see the entire line):: - 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 - 1L_EYKX4VC01BXWX1_265 1430919 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n Euarchontoglires Primates Haplorrhini Hominoidea Hominidae n n n Homo n Homo sapiens n + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 + 1L_EYKX4VC01BXWX1_265 9606 root Eukaryota Metazoa n n Chordata Craniata Gnathostomata Mammalia n Euarchontoglires Primates Haplorrhini Hominoidea Hominidae n n n Homo n Homo sapiens n 1430919 -In other words the tool printed *Name column*, *GI column*, and appended 22 columns containing taxonomic ranks from Superkingdom to Subspecies. Below is a formal definition of the output columns:: +In other words the tool printed *Name column*, *taxonomy Id*, appended 22 columns containing taxonomic ranks from Superkingdom to Subspecies and added *GI* as the last column. Below is a formal definition of the output columns:: Column Definition ------- ------------------------------------------ diff --git a/tools/taxonomy/tax.py b/tools/taxonomy/tax.py deleted file mode 100644 index bcd6295a90f..00000000000 --- a/tools/taxonomy/tax.py +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env python - -""" -Identify full taxonomic standing for sequences identified by gi number - -usage: tax.py gi_list_file out_file columnNumber - - gi_list_file - input file containing GI identifiers - out_file - output file - idColumnNumber - integer corresponding to column containing identifier that will be included into output - giColumnNumber - integer corresponding to column in gi_list_file containing GIs (column numbers start with 1) - -""" - -import pkg_resources -pkg_resources.require( 'bx-python' ) -pkg_resources.require( 'pysqlite' ) -import traceback -import fileinput -from pysqlite2 import dbapi2 as sqlite -from warnings import warn -import string, sys - - - -TAXONOMY = '/Users/anton/galaxy/static/taxonomy/taxonomy.db' - -# database containing collapsed NCBI taxonomy generated by prepareTaxonomy.sh script -# distributed with Galaxy. See prepareTaxonomy.readme (in scripts/taxonomy ditrectory) for information on how to generate -# necessary files - -def main(): - - try: - gi_fname = sys.argv[1] - out_fname = sys.argv[2] - idCol = int( sys.argv[3] ) - 1 - giCol = int( sys.argv[4] ) - 1 - except: - sys.stderr.write('Not enough arguments\n') - sys.exit(0) - - try: - con = sqlite.connect(TAXONOMY) - except: - sys.stderr.write('Cannot connect to database\n') - sys.exit(0) - - cur = con.cursor() - fg = open(gi_fname, 'r') - of = open( out_fname, "w" ) - - try: - for line in fg: - try: - field = string.split(line.rstrip(), '\t') - sqlTemplate = string.Template('select gi2tax.gi, tax.* from gi2tax left join tax on gi2tax.taxId = tax.taxId where gi2tax.gi = $gi') - sql = sqlTemplate.substitute(gi = int(field[giCol])) - cur.execute(sql) - - for item in cur.fetchall(): - ranks = string.split(item[2], ",") - print >> of, field[idCol] + "\t" + str(item[0]) + "\t" + "\t".join(ranks) - - except: - pass - - finally: - fg.close() - of.close() - -if __name__ == "__main__": - main()