forester/ruby/evoruby/lib/evo/tool/taxonomy_processor.rb

   1 #
   2 # = lib/evo/apps/taxonomy_processor - TaxonomyProcessor class
   3 #
   4 # Copyright::  Copyright (C) 2006-2007 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: taxonomy_processor.rb,v 1.26 2010/12/13 19:00:11 cmzmasek Exp $
   8
   9
  10 require 'lib/evo/util/util'
  11 require 'lib/evo/msa/msa_factory'
  12 require 'lib/evo/msa/msa'
  13 require 'lib/evo/io/msa_io'
  14 require 'lib/evo/io/parser/fasta_parser'
  15 require 'lib/evo/io/parser/general_msa_parser'
  16 require 'lib/evo/io/writer/fasta_writer'
  17 require 'lib/evo/io/writer/phylip_sequential_writer'
  18 require 'lib/evo/util/command_line_arguments'
  19
  20 module Evoruby
  21
  22   class TaxonomyProcessor
  23
  24     PRG_NAME       = "tap"
  25     PRG_DATE       = "2012.09.27"
  26     PRG_DESC       = "replacement of species names in multiple sequence files"
  27     PRG_VERSION    = "1.02"
  28     COPYRIGHT      = "2012 Christian M Zmasek"
  29     CONTACT        = "phylosoft@gmail.com"
  30     WWW            = "www.phylosoft.org"
  31
  32     SIMPLE = false
  33
  34     EXTRACT_TAXONOMY_OPTION = "t"
  35
  36     def initialize()
  37       @taxonomies = Hash.new()
  38     end
  39
  40     def run()
  41
  42       Util.print_program_information( PRG_NAME,
  43         PRG_VERSION,
  44         PRG_DESC,
  45         PRG_DATE,
  46         COPYRIGHT,
  47         CONTACT,
  48         WWW,
  49         STDOUT )
  50
  51       if ( ARGV == nil || ( ARGV.length != 1 && ARGV.length != 2 && ARGV.length != 3 && ARGV.length != 4 && ARGV.length != 5 && ARGV.length != 6 ) )
  52         puts( "Usage: #{PRG_NAME}.rb [options] [input map file] <input sequences> [output sequences] [output id list]" )
  53         puts()
  54         puts( "  options: -" + EXTRACT_TAXONOMY_OPTION + ": to extract taxonomy information from bracketed expression" )
  55         puts()
  56         exit( -1 )
  57       end
  58
  59       begin
  60         cla = CommandLineArguments.new( ARGV )
  61       rescue ArgumentError => e
  62         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
  63       end
  64
  65
  66       mapfile   = nil
  67       input     = nil
  68       output    = nil
  69       list_file = nil
  70
  71
  72
  73       if cla.get_number_of_files == 4
  74         mapfile   = cla.get_file_name( 0 )
  75         input     = cla.get_file_name( 1 )
  76         output    = cla.get_file_name( 2 )
  77         list_file = cla.get_file_name( 3 )
  78       elsif cla.get_number_of_files == 3
  79         input     = cla.get_file_name( 0 )
  80         output    = cla.get_file_name( 1 )
  81         list_file = cla.get_file_name( 2 )
  82       elsif cla.get_number_of_files == 1
  83         input     = cla.get_file_name( 0 )
  84         i = nil
  85         if input.downcase.end_with?( ".fasta" )
  86           i = input[ 0 .. input.length - 7 ]
  87         elsif input.downcase.end_with?( ".fsa" )
  88           i = input[ 0 .. input.length - 5 ]
  89         else
  90           i = input
  91         end
  92         output    = i + "_ni.fasta"
  93         list_file = i + ".nim"
  94       end
  95
  96
  97       allowed_opts = Array.new
  98       allowed_opts.push( EXTRACT_TAXONOMY_OPTION )
  99
 100       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
 101       if ( disallowed.length > 0 )
 102         Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
 103       end
 104
 105       extract_taxonomy = false
 106       if ( cla.is_option_set?( EXTRACT_TAXONOMY_OPTION ) )
 107         extract_taxonomy = true
 108       end
 109
 110       if ( File.exists?( output ) )
 111         Util.fatal_error( PRG_NAME, "outfile [" + output + "] already exists" )
 112       end
 113       if ( File.exists?( list_file ) )
 114         Util.fatal_error( PRG_NAME, "list file [" + list_file + "] already exists" )
 115       end
 116       if ( !File.exists?( input) )
 117         Util.fatal_error( PRG_NAME, "infile [" + input + "] does not exist" )
 118       end
 119       if ( mapfile != nil && !File.exists?( mapfile ) )
 120         Util.fatal_error( PRG_NAME, "mapfile [" + mapfile + "] does not exist" )
 121       end
 122
 123       fasta_like = Util.looks_like_fasta?( input )
 124
 125       puts()
 126       if mapfile != nil
 127         puts( "Map file        : " + mapfile )
 128       end
 129       puts( "Input alignment : " + input )
 130       puts( "Output alignment: " + output )
 131       puts( "Name list       : " + list_file )
 132       if ( fasta_like )
 133         puts( "Format          : Fasta"  )
 134       else
 135         puts( "Format          : Phylip like" )
 136       end
 137       if ( extract_taxonomy )
 138         puts( "Extract taxonomy: true"  )
 139       end
 140       puts()
 141
 142       species_map = Hash.new
 143       if mapfile != nil
 144         File.open( mapfile ) do | file |
 145           while line = file.gets
 146             if ( line =~/(.+)#(.+)/ || line =~/(.+)\s+(.+)/ )
 147               species_map[ $1 ] = $2
 148               Util.print_message( PRG_NAME, "mapping: " + $1 + ' => ' + $2 )
 149             end
 150           end
 151         end
 152       end
 153
 154       f = MsaFactory.new()
 155       begin
 156         if ( fasta_like )
 157           msa = f.create_msa_from_file( input, FastaParser.new() )
 158         else
 159           msa = f.create_msa_from_file( input, GeneralMsaParser.new() )
 160         end
 161       rescue Exception => e
 162         Util.fatal_error( PRG_NAME, "failed to read file: " + e.to_s )
 163       end
 164
 165       if ( msa == nil || msa.get_number_of_seqs() < 1 )
 166         Util.fatal_error( PRG_NAME, "failed to read MSA" )
 167       end
 168       begin
 169         Util.check_file_for_writability( list_file )
 170       rescue Exception => e
 171         Util.fatal_error( PRG_NAME, "error: " + e.to_, STDOUT )
 172       end
 173
 174       #removed = msa.remove_redundant_sequences!( true )
 175       #if removed.size > 0
 176       #  Util.print_message( PRG_NAME, "going to ignore the following " + removed.size.to_s + " redundant sequences:" )
 177       #  removed.each { | seq_name |
 178       #    puts seq_name
 179       #  }
 180       #  Util.print_message( PRG_NAME, "will process " + msa.get_number_of_seqs.to_s + " non redundant sequences" )
 181       #end
 182
 183       lf = File.open( list_file, "a" )
 184       for i in 0 ... msa.get_number_of_seqs
 185         seq  = msa.get_sequence( i )
 186         seq.set_name( Util::normalize_seq_name( modify_name( seq.get_name(), i, lf, species_map, extract_taxonomy ), 10 ) )
 187       end
 188
 189       io = MsaIO.new()
 190       w = nil
 191       if ( fasta_like )
 192         w = FastaWriter.new()
 193       else
 194         w = PhylipSequentialWriter.new()
 195       end
 196       w.set_max_name_length( 10 )
 197       w.clean( true )
 198       begin
 199         io.write_to_file( msa, output, w )
 200       rescue Exception => e
 201         Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
 202       end
 203       lf.close()
 204       if ( @taxonomies.length > 0 )
 205         Util.print_message( PRG_NAME, "number of unique taxonomies: " + @taxonomies.length.to_s )
 206       end
 207       Util.print_message( PRG_NAME, "wrote: " + list_file )
 208       Util.print_message( PRG_NAME, "wrote: " + output )
 209       Util.print_message( PRG_NAME, "OK" )
 210     end
 211
 212     private
 213
 214     def modify_name( desc, counter, file, species_map, extract_taxonomy )
 215       new_desc = nil
 216       my_species = nil
 217       # if desc =~ /^>?\s*\S{1,10}_([0-9A-Z]{3,5})/
 218       desc.gsub!( /:\s+/, ":" ) #new
 219       desc.gsub!( /\s+/, " " ) #new
 220       if desc =~ /^>?\s*\S{1,10}_([A-Z]{3,5})/
 221         new_desc = counter.to_s( 16 ) + "_" + $1
 222       elsif SIMPLE
 223         new_desc = counter.to_s( 16 )
 224       elsif extract_taxonomy
 225
 226         species = nil
 227         species_map.each_key do | key |
 228           if desc =~ /[\b|_]#{key}\b/  # Added boundaries to prevent e.g. RAT matching ARATH.
 229             species = species_map[ key ]
 230             new_desc = counter.to_s( 16 ) + "_" + species
 231             break
 232           end
 233         end
 234         if species == nil
 235           #if desc =~/.*\[(\S{3,}?)\]/
 236           if desc =~/\[([A-Z0-9]{3,6})\]\s*$/ #new
 237             species = $1
 238             species.strip!
 239             species.upcase!
 240             species.gsub!( /\s+/, " " )
 241             species.gsub!( /-/, "" )
 242             species.gsub!( /\)/, "" )
 243             species.gsub!( /\(/, "" )
 244             species.gsub!( /\'/, "" )
 245             if species =~ /\S+\s\S+/ || species =~ /\S{3,5}/
 246               if species =~ /(\S+)\s(\S+)/
 247                 code = $1[ 0..2 ] + $2[ 0..1 ]
 248               elsif  species =~ /\S{3,5}/
 249                 code = species
 250               elsif species.count( " " ) > 2
 251                 species =~ /(\S+)\s+(\S+)\s+(\S+)$/
 252                 third_last = $1
 253                 second_last = $2
 254                 last = $3
 255                 code = code[ 0 ] + third_last[ 0 ] + second_last[ 0 ] + last[ 0 ] + last[ last.size - 1 ]
 256               elsif species.count( " " ) > 1
 257                 species =~ /(\S+)\s+(\S+)$/
 258                 second_last = $1
 259                 last = $2
 260                 code = code[ 0..1 ] + second_last[ 0 ] + last[ 0 ] + last[ last.size - 1 ]
 261               end
 262               new_desc = counter.to_s( 16 ) + "_" + code
 263               if @taxonomies.has_key?( code )
 264                 if ( !@taxonomies.has_value?( species ) )
 265                   Util.fatal_error( PRG_NAME, "code [#{code}] is not unique in [#{desc}]" )
 266                 end
 267               else
 268                 if ( @taxonomies.has_value?( species ) )
 269                   Util.fatal_error( PRG_NAME, "genome [#{species}] is not unique in [#{desc}]" )
 270                 else
 271                   @taxonomies[ code ] = species
 272                 end
 273               end
 274             else
 275               Util.fatal_error( PRG_NAME, "illegal format [#{species}] in: " + desc )
 276             end
 277           else
 278             Util.fatal_error( PRG_NAME, "illegal format in: " + desc )
 279           end
 280         end
 281       else
 282         species = nil
 283         my_species = nil
 284         species_map.each_key do | key |
 285           if desc =~ /#{key}/
 286             species = species_map[ key ]
 287             species = species.gsub( /\s+/, "" )
 288             species = species.gsub( /_/, " " )
 289             my_species = species
 290             if species =~ /(\S+)\s+(\S+)/
 291               species = $1[0..2] + $2[0..1]
 292             end
 293             species = species.gsub( /\s+/, "" )
 294             species = species.slice(0, 5)
 295             species.upcase!
 296             break
 297           end
 298         end
 299         if species == nil
 300           Util.fatal_error( PRG_NAME, "species not found in: " + desc  )
 301         else
 302           new_desc = counter.to_s( 16 ) + "_" + species
 303         end
 304       end
 305       if new_desc == nil
 306         Util.fatal_error( PRG_NAME, "failed to extract species from: " + desc  )
 307       end
 308       if my_species != nil
 309         file.print( new_desc + ": " + desc + " [" + my_species + "]" + "\n" )
 310       else
 311         file.print( new_desc + ": " + desc + "\n" )
 312       end
 313       new_desc
 314     end
 315
 316   end # class TaxonomyProcessor
 317
 318 end # module Evoruby