diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 47cf0c1e734..e05bdfdfe6a 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -155,6 +155,7 @@ +
diff --git a/tools/taxonomy/lca.py b/tools/taxonomy/lca.py index 0b9785d9dd5..d70425db366 100644 --- a/tools/taxonomy/lca.py +++ b/tools/taxonomy/lca.py @@ -1,10 +1,9 @@ #!/usr/bin/env python #Guruprasad Ananda """ -This tool provides the SQL "group by" functionality. +Least Common Ancestor tool. """ import sys, string, re, commands, tempfile, random -#from rpy import * def stop_err(msg): sys.stderr.write(msg) @@ -42,8 +41,8 @@ def main(): """ except: stop_err("Syntax error: Use correct syntax: program infile outfile") + group_col = 0 - tmpfile = tempfile.NamedTemporaryFile() try: @@ -68,10 +67,6 @@ def main(): prev_vals = [] remaining_vals = [] skipped_lines = 0 - first_invalid_line = 0 - invalid_line = '' - invalid_value = '' - invalid_column = 0 fout = open(outfile, "w") cols = range(1,25) block_valid = False @@ -104,9 +99,10 @@ def main(): out_list[0] = str(prev_item) out_list[1] = str(prev_vals[0][0]) out_list[2] = str(prev_vals[1][0]) - out_list[24] = str(prev_vals[23][0]) - #print >> fout, prev_vals - #sys.exit() + try: + out_list[24] = str(prev_vals[23][0]) + except: + pass for k, col in enumerate(cols): if col >= 3 and col < 24: if len(set(prev_vals[k])) == 1: @@ -116,17 +112,14 @@ def main(): while k < 23: out_list[k+1] = 'n' k += 1 - - # print >>fout, '\t'.join(out_list) if rank_bound == 0: - print >>fout, ''.join(out_list) - print 'n'*( 24 - rank_bound ) + print >>fout, '\t'.join(out_list).strip() + #print 'n'*( 24 - rank_bound ) else: - print '\t'.join(out_list[rank_bound:24]) + #print '\t'.join(out_list[rank_bound:24]) if ''.join(out_list[rank_bound:24]) != 'n'*( 24 - rank_bound ): - print >>fout, '\t'.join(out_list) - + print >>fout, '\t'.join(out_list).strip() block_valid = True prev_item = item @@ -145,21 +138,21 @@ def main(): val_list.append(fields[col].strip()) prev_vals.append(val_list) - except Exception, exc: + except: skipped_lines += 1 - if not first_invalid_line: - first_invalid_line = ii+1 else: skipped_lines += 1 - if not first_invalid_line: - first_invalid_line = ii+1 - + # Handle the last grouped value out_list = ['']*25 out_list[0] = str(prev_item) out_list[1] = str(prev_vals[0][0]) out_list[2] = str(prev_vals[1][0]) - out_list[24] = str(prev_vals[23][0]) + try: + out_list[24] = str(prev_vals[23][0]) + except: + pass + for k, col in enumerate(cols): if col >= 3 and col < 24: if len(set(prev_vals[k])) == 1: @@ -171,16 +164,15 @@ def main(): k += 1 if rank_bound == 0: - print >>fout, '\t'.join(out_list) + print >>fout, '\t'.join(out_list).strip() else: - print ''.join(out_list[rank_bound:24]) - print 'n'*( 24 - rank_bound ) + #print ''.join(out_list[rank_bound:24]) + #print 'n'*( 24 - rank_bound ) if ''.join(out_list[rank_bound:24]) != 'n'*( 24 - rank_bound ): - print >>fout, '\t'.join(out_list) + print >>fout, '\t'.join(out_list).strip() if skipped_lines > 0: - msg= "Skipped %d invalid lines starting with line %d. Value '%s' in column %d is not numeric." % ( skipped_lines, first_invalid_line, invalid_value, invalid_column ) - print msg + print "Skipped %d invalid lines." % ( skipped_lines ) if __name__ == "__main__": main() \ No newline at end of file diff --git a/tools/taxonomy/lca.xml b/tools/taxonomy/lca.xml index f608a16cfc3..ca5e246d41a 100644 --- a/tools/taxonomy/lca.xml +++ b/tools/taxonomy/lca.xml @@ -1,12 +1,12 @@ - + lca.py $input1 $out_file1 $rank_bound - - - + + + @@ -32,13 +32,54 @@ - + + + + + + + + - + **What it does** -When performing metagenomic analyses it is often necessary to identify sequence reads corresponding to a particular taxonomic group, or, in other words, diagnostic of a particular taxonomic rank. This utility performs this analysis. It takes data generated by *Taxonomy manipulation->Fetch Taxonomic Ranks* as input and outputs either a list of sequence reads unique to a particular taxonomic rank, or a list of taxonomic ranks and the count of unique reads corresponding to each rank. +This tool identifies the lowest taxonomic rank for which a mategenomic sequencing read is diagnostic. It takes datasets produced by *Fetch Taxonomic Ranks* tool (aka Taxonomy format) as the input. + +------- + +**Example** + +Suppose you have two reads, **read_1** and **read_2**, with the following taxonomic profiles (scroll sideways to see the entire dataset):: + + read_1 1 root superkingdom1 kingdom1 subkingdom1 superphylum1 phylum1 subphylum1 superclass1 class1 subclass1 superorder1 order1 suborder1 superfamily1 family1 subfamily1 tribe1 subtribe1 genus1 subgenus1 species1 subspecies1 + read_1 2 root superkingdom1 kingdom1 subkingdom1 superphylum1 phylum1 subphylum1 superclass1 class1 subclass1 superorder1 order1 suborder1 superfamily1 family1 subfamily1 tribe1 subtribe1 genus2 subgenus2 species2 subspecies2 + read_2 3 root superkingdom1 kingdom1 subkingdom1 superphylum1 phylum3 subphylum3 superclass3 class3 subclass3 superorder3 order3 suborder3 superfamily3 family3 subfamily3 tribe3 subtribe3 genus3 subgenus3 species3 subspecies3 + read_2 4 root superkingdom1 kingdom1 subkingdom1 superphylum1 phylum4 subphylum4 superclass4 class4 subclass4 superorder4 order4 suborder4 superfamily4 family4 subfamily4 tribe4 subtribe4 genus4 subgenus4 species4 subspecies4 + +For **read_1** taxonomic labels are consistent until the genus level, where the taxonomy splits into two branches, one ending with *subspecies1* and the other with *subspecies2*. This implies **that the lowest taxomomic rank read_1 can identify is SUBTRIBE**. Similarly, read_2 is diagnostic up until the **superphylum** level. As a results the output of this tool will be:: + + read_1 2 root superkingdom1 kingdom1 subkingdom1 superphylum1 phylum1 subphylum1 superclass1 class1 subclass1 superorder1 order1 suborder1 superfamily1 family1 subfamily1 tribe1 subtribe1 n n n n + read_2 3 root superkingdom1 kingdom1 subkingdom1 superphylum1 n n n n n n n n n n n n n n n n n + +where, **n** means *EMPTY*. + +-------- + +**What's up with the drop down?** + +Why do we need the *require the lowest rank to be at least* dropdown? Let's look at the above example again. Suppose you need to find only those reads that are diagnostic on at least phylum level. To do this you need to set the *require the lowest rank to be at least* to **phylum**. As a result your output will look like this:: + + read_1 2 root superkingdom1 kingdom1 subkingdom1 superphylum1 phylum1 subphylum1 superclass1 class1 subclass1 superorder1 order1 suborder1 superfamily1 family1 subfamily1 tribe1 subtribe1 n n n n + +.. class:: infomark + +Note, that **read_2** is now omitted as it matches two phyla (**phylum3** and **phylum4**) and therefore is not diagnostic (but rather cosmopolitan) on *phylum* level. + + + + \ No newline at end of file