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

   1 #
   2 # = lib/evo/apps/taxonomy_processor - TaxonomyProcessor class
   3 #
   4 # Copyright::    Copyright (C) 2017 Christian M. Zmasek
   5 # License::      GNU Lesser General Public License (LGPL)
   6
   7 require 'lib/evo/util/constants'
   8 require 'lib/evo/util/util'
   9 require 'lib/evo/msa/msa_factory'
  10 require 'lib/evo/msa/msa'
  11 require 'lib/evo/io/msa_io'
  12 require 'lib/evo/io/parser/fasta_parser'
  13 require 'lib/evo/io/parser/general_msa_parser'
  14 require 'lib/evo/io/writer/fasta_writer'
  15 require 'lib/evo/io/writer/phylip_sequential_writer'
  16 require 'lib/evo/util/command_line_arguments'
  17
  18 module Evoruby
  19   class TaxonomyProcessor
  20
  21     PRG_NAME       = "tap"
  22     PRG_DATE       = "170214"
  23     PRG_DESC       = "Replacement of labels in multiple sequence files"
  24     PRG_VERSION    = "2.004"
  25     WWW            = "https://sites.google.com/site/cmzmasek/home/software/forester"
  26
  27     EXTRACT_TAXONOMY_OPTION = "t"
  28     ANNOTATION_OPTION       = "a"
  29     HELP_OPTION_1           = "help"
  30     HELP_OPTION_2           = "h"
  31     def run()
  32
  33       Util.print_program_information( PRG_NAME,
  34       PRG_VERSION,
  35       PRG_DESC,
  36       PRG_DATE,
  37       WWW,
  38       STDOUT )
  39
  40       if ( ARGV == nil || ( ARGV.length < 1 ) )
  41         print_help()
  42         exit( -1 )
  43       end
  44
  45       begin
  46         cla = CommandLineArguments.new( ARGV )
  47       rescue ArgumentError => e
  48         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
  49       end
  50
  51       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
  52       cla.is_option_set?( HELP_OPTION_2 ) )
  53         print_help
  54         exit( 0 )
  55       end
  56
  57       input      = nil
  58       output     = nil
  59       list_file  = nil
  60
  61       if cla.get_number_of_files == 3
  62         input     = cla.get_file_name( 0 )
  63         output    = cla.get_file_name( 1 )
  64         list_file = cla.get_file_name( 2 )
  65       elsif cla.get_number_of_files == 1
  66         input     = cla.get_file_name( 0 )
  67         i = nil
  68         if input.downcase.end_with?( ".fasta" )
  69           i = input[ 0 .. input.length - 7 ]
  70         elsif input.downcase.end_with?( ".fsa" )
  71           i = input[ 0 .. input.length - 5 ]
  72         else
  73           i = input
  74         end
  75         output    = i + Constants::ID_NORMALIZED_FASTA_FILE_SUFFIX
  76         list_file = i + Constants::ID_MAP_FILE_SUFFIX
  77       else
  78         print_help()
  79         exit(-1)
  80       end
  81
  82       allowed_opts = Array.new
  83       allowed_opts.push( EXTRACT_TAXONOMY_OPTION )
  84       allowed_opts.push( ANNOTATION_OPTION )
  85
  86       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
  87       if ( disallowed.length > 0 )
  88         Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
  89       end
  90
  91       extract_taxonomy = false
  92       if ( cla.is_option_set?( EXTRACT_TAXONOMY_OPTION ) )
  93         extract_taxonomy = true
  94       end
  95
  96       annotation = nil
  97       if ( cla.is_option_set?( ANNOTATION_OPTION ) )
  98         annotation = cla.get_option_value( ANNOTATION_OPTION )
  99       end
 100
 101       if ( File.exist?( output ) )
 102         Util.fatal_error( PRG_NAME, "outfile [" + output + "] already exists" )
 103       end
 104       if ( File.exist?( list_file ) )
 105         Util.fatal_error( PRG_NAME, "list file [" + list_file + "] already exists" )
 106       end
 107       if ( !File.exist?( input) )
 108         Util.fatal_error( PRG_NAME, "infile [" + input + "] does not exist" )
 109       end
 110
 111       fasta_like = Util.looks_like_fasta?( input )
 112
 113       puts()
 114       puts( "Input alignment : " + input )
 115       puts( "Output alignment: " + output )
 116       puts( "Name list       : " + list_file )
 117       if ( fasta_like )
 118         puts( "Format          : Fasta"  )
 119       else
 120         puts( "Format          : Phylip like" )
 121       end
 122       if ( extract_taxonomy )
 123         puts( "Extract taxonomy: true"  )
 124       end
 125       if ( annotation != nil )
 126         puts( "Annotation      : " + annotation )
 127       end
 128       puts()
 129
 130       f = MsaFactory.new()
 131       begin
 132         if ( fasta_like )
 133           msa = f.create_msa_from_file( input, FastaParser.new() )
 134         else
 135           msa = f.create_msa_from_file( input, GeneralMsaParser.new() )
 136         end
 137       rescue Exception => e
 138         Util.fatal_error( PRG_NAME, "failed to read file: " + e.to_s )
 139       end
 140
 141       if ( msa == nil || msa.get_number_of_seqs() < 1 )
 142         Util.fatal_error( PRG_NAME, "failed to read MSA" )
 143       end
 144       begin
 145         Util.check_file_for_writability( list_file )
 146       rescue Exception => e
 147         Util.fatal_error( PRG_NAME, "error: " + e.to_, STDOUT )
 148       end
 149
 150       lf = File.open( list_file, "a" )
 151       for i in 0 ... msa.get_number_of_seqs
 152         seq = msa.get_sequence( i )
 153         seq.set_name( modify_name( seq.get_name(), i, lf, extract_taxonomy, annotation ) )
 154       end
 155       io = MsaIO.new()
 156       w = nil
 157       if ( fasta_like )
 158         w = FastaWriter.new()
 159       else
 160         w = PhylipSequentialWriter.new()
 161       end
 162       w.set_max_name_length( 9 )
 163       w.clean( true )
 164       begin
 165         io.write_to_file( msa, output, w )
 166       rescue Exception => e
 167         Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
 168       end
 169       lf.close()
 170       Util.print_message( PRG_NAME, "wrote: " + list_file )
 171       Util.print_message( PRG_NAME, "wrote: " + output )
 172       Util.print_message( PRG_NAME, "next steps in standard analysis pipeline: hmmscan followed by hsp.rb")
 173       Util.print_message( PRG_NAME, "hmmscan example: hmmscan --max --domtblout P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10 -E 10 Pfam-A.hmm P53_ni.fasta")
 174
 175       Util.print_message( PRG_NAME, "OK" )
 176     end
 177
 178     private
 179
 180     def modify_name( desc, counter, file, extract_taxonomy, annotation )
 181       new_desc = nil
 182       desc.gsub!( /\s+/, ' ' )
 183       if extract_taxonomy
 184         if desc =~/\s\[(([A-Z9][A-Z]{2}[A-Z0-9]{2})|RAT|PIG|PEA|CAP)\]/
 185           new_desc = counter.to_s( 16 ) + "_" + $1
 186         else
 187           Util.fatal_error( PRG_NAME, "could not get taxonomy from: " + desc )
 188         end
 189       else
 190         new_desc = counter.to_s( 16 )
 191       end
 192       if (annotation != nil)
 193         new_desc = new_desc + annotation
 194         file.print( new_desc + "\t" + desc + " " + annotation + "\n" )
 195       else
 196         file.print( new_desc + "\t" + desc + "\n" )
 197       end
 198       if ( new_desc.length > 9)
 199         Util.fatal_error( PRG_NAME, "shortened identifier [" +
 200         new_desc + "] is too long (" + new_desc.length.to_s + " characters)" )
 201       end
 202       new_desc
 203     end
 204
 205     def print_help()
 206       puts( "Usage:" )
 207       puts()
 208       puts( "  " + PRG_NAME + ".rb [options] <input sequences> [output sequences] [output id list]" )
 209       puts()
 210       puts( "  options: -" + EXTRACT_TAXONOMY_OPTION + "    : to extract taxonomy information from bracketed expressions" )
 211       puts( "           -" + ANNOTATION_OPTION + "=<s>: to add an annotation to all entries" )
 212       puts()
 213       puts( "  [next steps in standard analysis pipeline: hmmscan followed by hsp.rb]")
 214       puts( "  [hmmscan example: hmmscan --max --domtblout P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10 -E 10 Pfam-A.hmm P53_ni.fasta]")
 215       puts()
 216       puts( "Example:" )
 217       puts()
 218       puts( "  " + PRG_NAME + ".rb P53.fasta" )
 219       puts()
 220     end
 221
 222   end # class TaxonomyProcessor
 223
 224 end # module Evoruby