#NEWLINE# def main( argv = None ) : #NEWLINE# #INDENT# 'script main.\n\n    parses command line options in sys.argv, unless *argv* is given.\n    ' #NEWLINE# if ( not argv ) : #NEWLINE# #INDENT# argv = sys . argv #NEWLINE# #UNINDENT# parser = optparse . OptionParser ( version = '%prog version: $Id: script_template.py 2871 2010-03-03 10:20:44Z andreas $' , usage = globals ( ) [ '__doc__' ] ) #NEWLINE# parser . add_option ( '-g' , '--gi-accessions' , dest = 'gi_accessions' , type = 'string' , help = 'list of gi accession numbers' ) #NEWLINE# parser . add_option ( '-m' , '--ncbi-map' , dest = 'ncbi_map' , type = 'string' , help = 'ncbi.map file downloaded from the MEGAN website' ) #NEWLINE# parser . add_option ( '-n' , '--nucl-map' , dest = 'nucl_map' , type = 'string' , help = 'gi mapping to tax id downloaded from ncbi website' ) #NEWLINE# parser . add_option ( '-c' , '--taxa-code' , dest = 'taxa_code' , type = 'string' , help = 'code for different levels of the taxonomy downloaded from the MEGAN website' ) #NEWLINE# parser . add_option ( '-t' , '--tree' , dest = 'tree' , type = 'string' , help = 'description of parents in the taxonomy' ) #NEWLINE# ( options , args ) = E . Start ( parser , argv = argv ) #NEWLINE# E . info ( 'reading gi accession numbers' ) #NEWLINE# gi_accessions = set ( ) #NEWLINE# for line in open ( options . gi_accessions ) . readlines ( ) : #NEWLINE# #INDENT# gi_accessions . add ( line [ : ( - 1 ) ] ) #NEWLINE# #UNINDENT# E . info ( 'read gi accession numbers' ) #NEWLINE# E . info ( 'building gi2taxid map' ) #NEWLINE# gi2taxid = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# c_gi = 0 #NEWLINE# for line in IOTools . openFile ( options . nucl_map ) . readlines ( ) : #NEWLINE# #INDENT# data = line [ : ( - 1 ) ] . split ( '\t' ) #NEWLINE# if ( data [ 0 ] not in gi_accessions ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# c_gi += 1 #NEWLINE# gi2taxid [ data [ 0 ] ] = data [ 1 ] #NEWLINE# #UNINDENT# #UNINDENT# E . info ( ( 'built gi2taxid for %i gi accession numbers' % c_gi ) ) #NEWLINE# E . info ( 'building code map' ) #NEWLINE# code2taxa = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for line in open ( options . taxa_code ) . readlines ( ) : #NEWLINE# #INDENT# data = line [ : ( - 1 ) ] . split ( '\t' ) #NEWLINE# code2taxa [ data [ 0 ] ] = data [ 1 ] #NEWLINE# #UNINDENT# E . info ( 'built taxa code map' ) #NEWLINE# E . info ( 'building taxa2name map' ) #NEWLINE# taxid2name = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for line in open ( options . ncbi_map ) . readlines ( ) : #NEWLINE# #INDENT# data = line [ : ( - 1 ) ] . split ( '\t' ) #NEWLINE# taxid2name [ data [ 0 ] ] = ( data [ 1 ] , data [ 3 ] ) #NEWLINE# #UNINDENT# E . info ( 'built taxa2name map' ) #NEWLINE# E . info ( 'build taxid2parentmap' ) #NEWLINE# taxid2parents = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for line in open ( options . tree ) . readlines ( ) : #NEWLINE# #INDENT# data = line [ : ( - 1 ) ] . split ( '\t' ) #NEWLINE# data = [ x for x in data if ( x != '|' ) ] #NEWLINE# taxid2parents [ data [ 0 ] ] = data [ 1 ] #NEWLINE# #UNINDENT# E . info ( 'built taxid2parentmap' ) #NEWLINE# E . info ( 'retrieving parents for each gi accession number' ) #NEWLINE# options . stdout . write ( 'gi\tsub_species\tspecies\tgenus\tfamily\torder\tclass\tphylum\n' ) #NEWLINE# for ( gi , taxid ) in gi2taxid . iteritems ( ) : #NEWLINE# #INDENT# parents = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# sub_species = taxid2name [ taxid ] [ 0 ] #NEWLINE# for i in range ( len ( code2taxa . keys ( ) ) ) : #NEWLINE# #INDENT# parent_taxid = taxid2parents [ taxid ] #NEWLINE# parent_name = taxid2name [ parent_taxid ] [ 0 ] #NEWLINE# parent_code = taxid2name [ parent_taxid ] [ 1 ] #NEWLINE# if ( parent_code not in code2taxa . keys ( ) ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# parent_taxa = code2taxa [ parent_code ] #NEWLINE# parents [ parent_taxa ] = genus #NEWLINE# taxid = parent_taxid #NEWLINE# #UNINDENT# if ( 'genus' not in parents ) : #NEWLINE# #INDENT# genus = 'NA' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# genus = parents [ 'genus' ] #NEWLINE# #UNINDENT# if ( 'family' not in parents ) : #NEWLINE# #INDENT# family = 'NA' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# family = parents [ 'family' ] #NEWLINE# #UNINDENT# if ( 'order' not in parents ) : #NEWLINE# #INDENT# order = 'NA' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# order = parents [ 'order' ] #NEWLINE# #UNINDENT# if ( 'class' not in parents ) : #NEWLINE# #INDENT# _class = 'NA' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# _class = parents [ 'class' ] #NEWLINE# #UNINDENT# if ( 'phylum' not in parents ) : #NEWLINE# #INDENT# phylum = 'NA' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# phylum = parents [ 'phylum' ] #NEWLINE# if ( phylum . find ( '<phylum>' ) != ( - 1 ) ) : #NEWLINE# #INDENT# phylum = phylum . replace ( ' <phylum>' , '' ) #NEWLINE# #UNINDENT# #UNINDENT# if ( 'species' not in parents ) : #NEWLINE# #INDENT# species = 'NA' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# species = parents [ 'species' ] #NEWLINE# #UNINDENT# options . stdout . write ( ( '\t' . join ( [ gi , sub_species . replace ( ' ' , '_' ) , species . replace ( ' ' , '_' ) , genus , family , order , _class , phylum ] ) + '\n' ) ) #NEWLINE# #UNINDENT# E . Stop ( )