opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Anton Nekrutenko
2008-03-20 19:38:31 +00:00
parent 211673cf3b
commit a262ac1754
14 changed files with 2038 additions and 148 deletions
Download Patch File Download Diff File Expand all files Collapse all files
scripts/taxonomy/processTaxonomy.sh
+4 -7
View File
@@ -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
scripts/taxonomy/process_NCBI_taxonomy.py
-50
View File
@@ -1,50 +0,0 @@
"""
process_NCBI_taxonomy.py <gi2tax.txt file> <name.txt> <database_name>
"""
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 <gi2tax.txt file> <name.txt> <database_name>\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)
test-data/taxonomyGI.dat
-5
View File
@@ -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
tools/taxonomy/T2PS/t2ps.readme
+19 -2
View File
@@ -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
tools/taxonomy/T2T/t2t.readme
+18 -2
View File
@@ -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
tools/taxonomy/TB/avl.c
+890
View File
@@ -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 <blp@gnu.org> on the Internet, or
write to Ben Pfaff, Stanford University, Computer Science Dept., 353
Serra Mall, Stanford CA 94305, USA.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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 <assert.h>
/* 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;
}
tools/taxonomy/TB/avl.h
+115
View File
@@ -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 <blp@gnu.org> 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 <stddef.h>
/* 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 */
tools/taxonomy/TB/main.c
+788
View File
@@ -0,0 +1,788 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <ctype.h>
#include <string.h>
#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 <pl; addThis++)
d->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 <pl; addThis++)
d->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; k<level; k++)
appendCharBufferToString(resStr,noValue);
}
//appendCharacterToString(resStr,')');
//sprintf (buffer, "%d", currentTag->taxonomy_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; i<upTo; i++)
{
int res = (s1->sData[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->sLength<s2->sLength);
}
/*---------------------------------------------------------------------------------------------------- */
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; i<upTo; i++)
{
int res = (int)buffa[s1->startIndex+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->length<s2->length);
}
/*---------------------------------------------------------------------------------------------------- */
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;
}
tools/taxonomy/TB/tb.readme
+21
View File
@@ -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
tools/taxonomy/find_diag_hits.py
+4 -2
View File
@@ -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()
tools/taxonomy/find_diag_hits.xml
+1 -1
View File
@@ -24,7 +24,7 @@
<option value="subtribe">Subtribe</option>
<option value="genus">Genus</option>
<option value="subgenus">Subgenus</option>
<option value="species">Species</option>
<option selected="true" value="species">Species</option>
<option value="subspecies">Subspecies</option>
</param>
<param name="out_format" type="select" label="Select output format">
tools/taxonomy/gi2taxonomy.py
+172
View File
@@ -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
tools/taxonomy/gi2taxonomy.xml
+6 -6
View File
@@ -1,9 +1,9 @@
<tool id="Fetch Taxonomic Ranks" name="Fetch Taxonomic Ranks" version="1.0.0">
<tool id="Fetch Taxonomic Ranks" name="Fetch Taxonomic Ranks" version="1.1.0">
<description>for a list of GIs</description>
<command interpreter="python">tax.py $input $out_file1 $idField $giField</command>
<command interpreter="python">gi2taxonomy.py $input $giField $idField $out_file1</command>
<inputs>
<param format="tabular" name="input" type="data" label="Fetch sequences corresponding to Query"></param>
<param name="giField" label="GIs column" type="data_column" data_ref="input" help="select column containing GI numbers"/>
<param name="giField" label="GIs column" type="data_column" data_ref="input" numerical="True" help="select column containing GI numbers"/>
<param name="idField" label="Name column" type="data_column" data_ref="input" help="select column containing identifiers you want to include into output"/>
</inputs>
<outputs>
@@ -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
------- ------------------------------------------
tools/taxonomy/tax.py
-73
View File
@@ -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()