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