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

   1 #
   2 # = lib/evo/apps/msa_processor.rb - MsaProcessor class
   3 #
   4 # Copyright::  Copyright (C) 2006-2007 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: msa_processor.rb,v 1.33 2010/12/13 19:00:10 cmzmasek Exp $
   8 #
   9
  10 require 'date'
  11 require 'set'
  12 require 'lib/evo/util/constants'
  13 require 'lib/evo/util/util'
  14 require 'lib/evo/util/command_line_arguments'
  15 require 'lib/evo/msa/msa_factory'
  16 require 'lib/evo/io/msa_io'
  17 require 'lib/evo/io/writer/phylip_sequential_writer'
  18 require 'lib/evo/io/writer/nexus_writer'
  19 require 'lib/evo/io/writer/fasta_writer'
  20 require 'lib/evo/io/parser/fasta_parser'
  21 require 'lib/evo/io/parser/general_msa_parser'
  22 require 'lib/evo/io/writer/msa_writer'
  23
  24 module Evoruby
  25   class MsaProcessor
  26
  27     PRG_NAME       = "msa_pro"
  28     PRG_DATE       = "170609"
  29     PRG_DESC       = "processing of multiple sequence alignments"
  30     PRG_VERSION    = "1.09"
  31     WWW            = "https://sites.google.com/site/cmzmasek/home/software/forester"
  32
  33     NAME_LENGTH_DEFAULT                = 10
  34     WIDTH_DEFAULT_FASTA                = 60
  35     INPUT_TYPE_OPTION                  = "i"
  36     OUTPUT_TYPE_OPTION                 = "o"
  37     MAXIMAL_NAME_LENGTH_OPTION         = "n"
  38     DIE_IF_NAME_TOO_LONG               = "d"
  39     WIDTH_OPTION                       = "w"
  40     CLEAN_UP_SEQ_OPTION                = "c"
  41     REM_RED_OPTION                     = "rem_red"
  42     REMOVE_GAP_COLUMNS_OPTION          = "rgc"
  43     REMOVE_GAP_ONLY_COLUMNS            = "rgoc"
  44     REMOVE_COLUMNS_GAP_RATIO_OPTION    = "rr"
  45     REMOVE_ALL_GAP_CHARACTERS_OPTION   = "rg"
  46     REMOVE_ALL_SEQUENCES_LISTED_OPTION = "r"
  47     KEEP_ONLY_SEQUENCES_LISTED_OPTION  = "k"
  48
  49     KEEP_MATCHING_SEQUENCES_OPTION     = "mk"
  50     REMOVE_MATCHING_SEQUENCES_OPTION   = "mr"
  51
  52     TRIM_OPTION                        = "t"
  53     SLIDING_EXTRACTION_OPTION          = "se"
  54     REMOVE_SEQS_GAP_RATIO_OPTION       = "rsgr"
  55     REMOVE_SEQS_NON_GAP_LENGTH_OPTION  = "rsl"
  56     SPLIT                              = "split"
  57     SPLIT_BY_OS                        = "split_by_os"
  58     LOG_SUFFIX                         = "_msa_pro.log"
  59     HELP_OPTION_1                      = "help"
  60     HELP_OPTION_2                      = "h"
  61     def initialize()
  62       @input_format_set = false
  63       @output_format_set = false
  64       @fasta_input      = false
  65       @phylip_input     = true
  66       @name_length      = NAME_LENGTH_DEFAULT
  67       @name_length_set  = false
  68       @width            = WIDTH_DEFAULT_FASTA     # fasta only
  69       @pi_output        = true
  70       @fasta_output     = false
  71       @nexus_output     = false
  72       @clean            = false  # phylip only
  73       @rgc              = false
  74       @rgoc             = false
  75       @rg               = false  # fasta only
  76       @rem_red          = false
  77       @die_if_name_too_long  = false
  78       @rgr              = -1
  79       @rsgr             = -1
  80       @rsl              = -1
  81       @remove_matching  = nil
  82       @keep_matching    = nil
  83
  84       @seqs_name_file   = nil
  85       @remove_seqs      = false
  86       @keep_seqs        = false
  87       @trim             = false
  88       @split            = -1
  89       @split_by_os      = false
  90       @first            = -1
  91       @last             = -1
  92       @window = false
  93       @step = -1
  94       @size =  -1
  95     end
  96
  97     def run()
  98
  99       Util.print_program_information( PRG_NAME,
 100       PRG_VERSION,
 101       PRG_DESC,
 102       PRG_DATE,
 103       WWW,
 104       STDOUT )
 105
 106       if ( ARGV == nil || ARGV.length < 1 )
 107         Util.print_message( PRG_NAME, "Illegal number of arguments" )
 108         print_help
 109         exit( -1 )
 110       end
 111
 112       begin
 113         cla = CommandLineArguments.new( ARGV )
 114       rescue ArgumentError => e
 115         Util.fatal_error( PRG_NAME, "Error: " + e.to_s, STDOUT )
 116       end
 117
 118       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
 119       cla.is_option_set?( HELP_OPTION_2 ) )
 120         print_help
 121         exit( 0 )
 122       end
 123
 124       if ( cla.get_number_of_files != 2 || ARGV.length < 2 )
 125         Util.print_message( PRG_NAME, "Illegal number of arguments" )
 126         print_help
 127         exit( -1 )
 128       end
 129
 130       allowed_opts = Array.new
 131       allowed_opts.push( INPUT_TYPE_OPTION )
 132       allowed_opts.push( OUTPUT_TYPE_OPTION )
 133       allowed_opts.push( MAXIMAL_NAME_LENGTH_OPTION )
 134       allowed_opts.push( WIDTH_OPTION )
 135       allowed_opts.push( CLEAN_UP_SEQ_OPTION )
 136       allowed_opts.push( REMOVE_GAP_COLUMNS_OPTION )
 137       allowed_opts.push( REMOVE_GAP_ONLY_COLUMNS )
 138       allowed_opts.push( REMOVE_COLUMNS_GAP_RATIO_OPTION )
 139       allowed_opts.push( REMOVE_ALL_GAP_CHARACTERS_OPTION )
 140       allowed_opts.push( REMOVE_ALL_SEQUENCES_LISTED_OPTION )
 141       allowed_opts.push( KEEP_ONLY_SEQUENCES_LISTED_OPTION )
 142       allowed_opts.push( TRIM_OPTION )
 143       allowed_opts.push( REMOVE_SEQS_GAP_RATIO_OPTION )
 144       allowed_opts.push( REMOVE_SEQS_NON_GAP_LENGTH_OPTION )
 145       allowed_opts.push( SPLIT )
 146       allowed_opts.push( SPLIT_BY_OS )
 147       allowed_opts.push( REM_RED_OPTION )
 148       allowed_opts.push( KEEP_MATCHING_SEQUENCES_OPTION )
 149       allowed_opts.push( REMOVE_MATCHING_SEQUENCES_OPTION )
 150       allowed_opts.push( DIE_IF_NAME_TOO_LONG )
 151       allowed_opts.push( SLIDING_EXTRACTION_OPTION )
 152
 153       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
 154       if ( disallowed.length > 0 )
 155         Util.fatal_error( PRG_NAME,
 156         "unknown option(s): " + disallowed )
 157       end
 158
 159       input = cla.get_file_name( 0 )
 160       output = cla.get_file_name( 1 )
 161
 162       analyze_command_line( cla )
 163
 164       begin
 165         Util.check_file_for_readability( input )
 166       rescue IOError => e
 167         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 168       end
 169
 170       begin
 171         Util.check_file_for_writability( output )
 172       rescue IOError => e
 173         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 174       end
 175
 176       if ( @rg )
 177         set_pi_output( false )
 178         set_fasta_output( true )
 179         set_nexus_output( false )
 180       end
 181
 182       if ( !@input_format_set )
 183         fasta_like = false
 184         begin
 185           fasta_like = Util.looks_like_fasta?( input )
 186         rescue ArgumentError => e
 187           Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 188         end
 189         @fasta_input = fasta_like
 190         @phylip_input = !fasta_like
 191         if ( !@output_format_set )
 192           @fasta_output = fasta_like
 193           @pi_output = !fasta_like
 194           @nexus_output = false
 195         end
 196       end
 197
 198       ld = Constants::LINE_DELIMITER
 199       log = PRG_NAME + " " + PRG_VERSION + " [" + PRG_DATE + "]" + " LOG" + ld
 200       now = DateTime.now
 201       log << "Date/time: " + now.to_s + ld
 202
 203       puts()
 204       puts( "Input alignment  : " + input )
 205       log << "Input alignment  : " + input + ld
 206       puts( "Output alignment : " + output )
 207       log << "Output alignment : " + output + ld
 208       if ( @phylip_input )
 209         puts( "Input is         : Phylip, or something like it" )
 210         log << "Input is         : Phylip, or something like it" + ld
 211       elsif ( @fasta_input )
 212         puts( "Input is         : Fasta" )
 213         log << "Input is         : Fasta" + ld
 214       end
 215       if( @rgr >= 0 )
 216         puts( "Max col gap ratio: " + @rgr.to_s )
 217         log << "Max col gap ratio: " + @rgr.to_s + ld
 218       elsif ( @rgc )
 219         puts( "Remove gap colums" )
 220         log << "Remove gap colums" + ld
 221       elsif( @rgoc )
 222         puts( "Remove gap only colums" )
 223         log << "Remove gap only colums" + ld
 224       end
 225       if ( @clean )
 226         puts( "Clean up         : true" )
 227         log << "Clean up         : true" + ld
 228       end
 229
 230       if ( @pi_output )
 231         puts( "Output is        : Phylip interleaved" )
 232         log << "Output is        : Phylip interleaved" + ld
 233       elsif ( @fasta_output )
 234         puts( "Output is        : Fasta" )
 235         log << "Output is        : Fasta" + ld
 236         if ( @width )
 237           puts( "Width            : " + @width.to_s )
 238           log << "Width            : " + @width.to_s + ld
 239         end
 240         if ( @rg )
 241           puts( "Remove all gap characters (alignment is destroyed)" )
 242           log << "Remove all gap characters (alignment is destroyed)" + ld
 243         end
 244       elsif ( @nexus_output )
 245         puts( "Output is        : Nexus" )
 246         log << "Output is        : Nexus" + ld
 247       end
 248       if ( @name_length_set || !@fasta_output )
 249         puts( "Max name length  : " + @name_length.to_s )
 250         log << "Max name length  : " + @name_length.to_s + ld
 251       end
 252       if( @rsgr >= 0 )
 253         puts( "Remove sequences for which the gap ratio > " + @rsgr.to_s )
 254         log << "Remove sequences for which the gap ratio > " + @rsgr.to_s + ld
 255       end
 256       if( @rsl >= 0 )
 257         puts( "Remove sequences with less than "  + @rsl.to_s + " non-gap characters" )
 258         log << "Remove sequences with less than "  + @rsl.to_s + " non-gap characters" + ld
 259       end
 260       if ( @remove_seqs )
 261         puts( "Remove sequences listed in: " + @seqs_name_file )
 262         log << "Remove sequences listed in: " + @seqs_name_file + ld
 263       elsif ( @keep_seqs )
 264         puts( "Keep only sequences listed in: " + @seqs_name_file )
 265         log << "Keep only sequences listed in: " + @seqs_name_file + ld
 266       end
 267       if ( @trim )
 268         puts( "Keep only columns from: "+ @first.to_s + " to " + @last.to_s )
 269         log << "Keep only columns from: "+ @first.to_s + " to " + @last.to_s + ld
 270       end
 271       if ( @rem_red )
 272         puts( "Remove redundant sequences: true" )
 273         log << "Remove redundant sequences: true" + ld
 274       end
 275       if ( @split_by_os )
 276         puts( "Split by OS      : true"  )
 277         log << "Split            : true" + ld
 278       end
 279       if ( @split > 0 )
 280         puts( "Split            : " + @split.to_s )
 281         log << "Split            : " + @split.to_s + ld
 282       end
 283       if @window
 284         puts( "Sliding window extraction: true"  )
 285         log << "Sliding window extraction: true" + ld
 286         puts( "Sliding window step      : " + @step.to_s )
 287         log << "Sliding window step      : " + @step.to_s + ld
 288         puts( "Sliding window size      : " +  @size.to_s )
 289         log << "Sliding window size      : " +  @size.to_s + ld
 290       end
 291
 292       puts()
 293
 294       f = MsaFactory.new()
 295
 296       msa = nil
 297
 298       begin
 299         if ( @phylip_input )
 300           msa = f.create_msa_from_file( input, GeneralMsaParser.new() )
 301         elsif ( @fasta_input )
 302           msa = f.create_msa_from_file( input, FastaParser.new() )
 303         end
 304       rescue Exception => e
 305         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 306       end
 307
 308       if ( msa.is_aligned() )
 309         Util.print_message( PRG_NAME, "Length of original alignment         : " + msa.get_length.to_s )
 310         log << "Length of original alignment         : " + msa.get_length.to_s + ld
 311         gp = msa.calculate_gap_proportion
 312         Util.print_message( PRG_NAME, "Gap-proportion of original alignment : " + gp.to_s )
 313         log << "Gap-proportion of original alignment : " +  gp.to_s + ld
 314       else
 315         Util.print_message( PRG_NAME, "Input is not aligned" )
 316         log << "Input is not aligned" + ld
 317       end
 318
 319       all_names = Set.new()
 320       for i in 0 ... msa.get_number_of_seqs()
 321         current_name = msa.get_sequence( i ).get_name
 322         if all_names.include?( current_name )
 323           Util.print_warning_message( PRG_NAME, "sequence name [" + current_name + "] is not unique" )
 324         else
 325           all_names.add( current_name )
 326         end
 327       end
 328
 329       begin
 330
 331         if ( @remove_seqs || @keep_seqs )
 332           names = Util.file2array( @seqs_name_file, true )
 333           if ( names == nil ||  names.length() < 1 )
 334             error_msg = "file \"" + @seqs_name_file.to_s + "\" appears empty"
 335             Util.fatal_error( PRG_NAME, error_msg )
 336           end
 337
 338           if ( @remove_seqs )
 339             c = 0
 340             for i in 0 ... names.length()
 341               to_delete = msa.find_by_name( names[ i ], true, false )
 342               if ( to_delete.length() < 1 )
 343                 error_msg = "sequence name \"" + names[ i ] + "\" not found"
 344                 Util.fatal_error( PRG_NAME, error_msg )
 345               elsif ( to_delete.length() > 1 )
 346                 error_msg = "sequence name \"" + names[ i ] + "\" is not unique"
 347                 Util.fatal_error( PRG_NAME, error_msg )
 348               else
 349                 msa.remove_sequence!( to_delete[ 0 ] )
 350                 c += 1
 351               end
 352             end
 353             Util.print_message( PRG_NAME, "Removed " + c.to_s + " sequences" )
 354             log <<  "Removed " + c.to_s + " sequences" + ld
 355           elsif ( @keep_seqs )
 356             msa_new = Msa.new()
 357             r = 0
 358             k = 0
 359             for j in 0 ... msa.get_number_of_seqs()
 360               if ( names.include?( msa.get_sequence( j ).get_name() ) )
 361                 msa_new.add_sequence( msa.get_sequence( j ) )
 362                 k += 1
 363               else
 364                 r += 1
 365               end
 366             end
 367             msa = msa_new
 368             Util.print_message( PRG_NAME, "Kept    " + k.to_s + " sequences" )
 369             log << "Kept    " + k.to_s + " sequences" + ld
 370             Util.print_message( PRG_NAME, "Removed " + r.to_s + " sequences" )
 371             log << "removed " + r.to_s + " sequences" + ld
 372           end
 373         end
 374
 375         if ( @trim )
 376           msa.trim!( @first, @last, '_S' )
 377         end
 378
 379         if @window
 380           msas = msa.sliding_extraction( @step, @size )
 381           begin
 382             io = MsaIO.new()
 383             w = MsaWriter
 384             if ( @pi_output )
 385               w = PhylipSequentialWriter.new()
 386               w.clean( @clean )
 387               w.set_max_name_length( @name_length )
 388             elsif( @fasta_output )
 389               w = FastaWriter.new()
 390               w.set_line_width( @width )
 391               w.clean( @clean )
 392               if ( @name_length_set )
 393                 w.set_max_name_length( @name_length )
 394               end
 395             elsif( @nexus_output )
 396               w = NexusWriter.new()
 397               w.clean( @clean )
 398               w.set_max_name_length( @name_length )
 399             end
 400             i = 0
 401             for m in msas
 402               i = i + 1
 403               name = output + "_" + m.get_name
 404               if @fasta_output
 405                 name += ".fasta"
 406               elsif @nexus_output
 407                 name += ".nex"
 408               end
 409               io.write_to_file( m, name, w )
 410             end
 411             Util.print_message( PRG_NAME, "wrote " + msas.length.to_s + " files"  )
 412             log << "wrote " + msas.length.to_s + " files" + ld
 413           rescue Exception => e
 414             Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 415           end
 416         end
 417
 418         if( @rgr >= 0 )
 419           msa.remove_gap_columns_w_gap_ratio!( @rgr )
 420         elsif ( @rgc )
 421           msa.remove_gap_columns!()
 422         elsif( @rgoc )
 423           msa.remove_gap_only_columns!()
 424         end
 425         if( @rsgr >= 0 )
 426           n = msa.get_number_of_seqs()
 427           removed = msa.remove_sequences_by_gap_ratio!( @rsgr )
 428           k = msa.get_number_of_seqs()
 429           r = n - k
 430           Util.print_message( PRG_NAME, "Kept    " + k.to_s + " sequences" )
 431           log << "Kept    " + k.to_s + " sequences" + ld
 432           Util.print_message( PRG_NAME, "Removed " + r.to_s + " sequences"  )
 433           log << "Removed " + r.to_s + " sequences:" + ld
 434           removed.each { | seq_name |
 435             log << "         " + seq_name  + ld
 436           }
 437         end
 438         if( @rsl >= 0 )
 439           n = msa.get_number_of_seqs()
 440           removed = msa.remove_sequences_by_non_gap_length!( @rsl )
 441           k = msa.get_number_of_seqs()
 442           r = n - k
 443           Util.print_message( PRG_NAME, "Kept    " + k.to_s + " sequences" )
 444           log << "Kept    " + k.to_s + " sequences" + ld
 445           Util.print_message( PRG_NAME, "Removed " + r.to_s + " sequences" )
 446           log << "Removed " + r.to_s + " sequences:" + ld
 447           removed.each { | seq_name |
 448             log << "         " + seq_name  + ld
 449           }
 450         end
 451         if ( @keep_matching )
 452           n = msa.get_number_of_seqs
 453           to_be_removed = Set.new
 454           for ii in 0 ...  n
 455             seq = msa.get_sequence( ii )
 456             if !seq.get_name.downcase.index( @keep_matching.downcase )
 457               to_be_removed.add( ii )
 458             end
 459           end
 460           to_be_removed_ary = to_be_removed.to_a.sort.reverse
 461           to_be_removed_ary.each { | index |
 462             msa.remove_sequence!( index )
 463           }
 464           # msa = sort( msa )
 465         end
 466         if ( @remove_matching )
 467           n = msa.get_number_of_seqs
 468           to_be_removed = Set.new
 469           for iii in 0 ... n
 470
 471             seq = msa.get_sequence( iii )
 472
 473             if seq.get_name.downcase.index( @remove_matching.downcase )
 474               to_be_removed.add( iii )
 475             end
 476           end
 477           to_be_removed_ary = to_be_removed.to_a.sort.reverse
 478           to_be_removed_ary.each { | index |
 479             msa.remove_sequence!( index )
 480           }
 481           msa = sort( msa )
 482         end
 483
 484         if ( @split_by_os )
 485           begin
 486             msa_hash = msa.split_by_os(true)
 487             io = MsaIO.new()
 488             w = MsaWriter
 489             if ( @pi_output )
 490               w = PhylipSequentialWriter.new()
 491               w.clean( @clean )
 492               w.set_max_name_length( @name_length )
 493             elsif( @fasta_output )
 494               w = FastaWriter.new()
 495               w.set_line_width( @width )
 496               if ( @rg )
 497                 w.remove_gap_chars( true )
 498                 Util.print_warning_message( PRG_NAME, "removing gap character, the output is likely to become unaligned" )
 499                 log << "removing gap character, the output is likely to become unaligned" + ld
 500               end
 501               w.clean( @clean )
 502               if ( @name_length_set )
 503                 w.set_max_name_length( @name_length )
 504               end
 505             elsif( @nexus_output )
 506               w = NexusWriter.new()
 507               w.clean( @clean )
 508               w.set_max_name_length( @name_length )
 509             end
 510             msa_hash.each do |os, m|
 511               my_os = os.gsub(' ', '_').gsub('/', '_').gsub('(', '_').gsub(')', '_')
 512               io.write_to_file( m, output + '_' + my_os, w )
 513             end
 514
 515             Util.print_message( PRG_NAME, "wrote " + msa_hash.length.to_s + " files"  )
 516             log << "wrote " + msa_hash.length.to_s + " files" + ld
 517           rescue Exception => e
 518             Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 519           end
 520
 521         elsif ( @split > 0 )
 522           begin
 523             msas = msa.split( @split, true )
 524             io = MsaIO.new()
 525             w = MsaWriter
 526             if ( @pi_output )
 527               w = PhylipSequentialWriter.new()
 528               w.clean( @clean )
 529               w.set_max_name_length( @name_length )
 530             elsif( @fasta_output )
 531               w = FastaWriter.new()
 532               w.set_line_width( @width )
 533               if ( @rg )
 534                 w.remove_gap_chars( true )
 535                 Util.print_warning_message( PRG_NAME, "removing gap character, the output is likely to become unaligned" )
 536                 log << "removing gap character, the output is likely to become unaligned" + ld
 537               end
 538               w.clean( @clean )
 539               if ( @name_length_set )
 540                 w.set_max_name_length( @name_length )
 541               end
 542             elsif( @nexus_output )
 543               w = NexusWriter.new()
 544               w.clean( @clean )
 545               w.set_max_name_length( @name_length )
 546             end
 547             i = 0
 548             for m in msas
 549               i = i + 1
 550               io.write_to_file( m, output + "_" + i.to_s, w )
 551             end
 552             Util.print_message( PRG_NAME, "wrote " + msas.length.to_s + " files"  )
 553             log << "wrote " + msas.length.to_s + " files" + ld
 554           rescue Exception => e
 555             Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 556           end
 557         end
 558       rescue Exception => e
 559         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 560       end
 561
 562       if (@split <= 0) && (!@split_by_os) && (!@window)
 563
 564         unless ( @rg )
 565           if ( msa.is_aligned() )
 566             Util.print_message( PRG_NAME, "Length of processed alignment        : " + msa.get_length.to_s )
 567             log <<  "Length of processed alignment        : " + msa.get_length.to_s + ld
 568             gp = msa.calculate_gap_proportion
 569             Util.print_message( PRG_NAME, "Gap-proportion of processed alignment: " + gp.to_s )
 570             log << "Gap-proportion of processed alignment: " +  gp.to_s + ld
 571           else
 572             min = 0
 573             max = 0
 574             sum = 0
 575             first = true
 576             for s in 0 ... msa.get_number_of_seqs
 577               seq = msa.get_sequence( s )
 578               l = seq.get_length
 579               sum += l
 580               if l > max
 581                 max = l
 582               end
 583               if first || l < min
 584                 min = l
 585               end
 586               first = false
 587             end
 588             avg = sum / msa.get_number_of_seqs
 589             Util.print_message( PRG_NAME, "Output is not aligned" )
 590             log << "Output is not aligned" + ld
 591             Util.print_message( PRG_NAME, "Shortest sequence                    : " + min.to_s )
 592             log <<  "Shortest sequence                    : " + min.to_s + ld
 593             Util.print_message( PRG_NAME, "Longest sequence                     : " + max.to_s )
 594             log <<  "Longest sequence                     : " + max.to_s + ld
 595             Util.print_message( PRG_NAME, "Average length                       : " + avg.to_s )
 596             log <<  "Average length                       : " + avg.to_s + ld
 597
 598           end
 599         end
 600
 601         if @rem_red
 602           removed = msa.remove_redundant_sequences!( true, false )
 603           if removed.size > 0
 604             identicals = msa.get_identical_seqs_detected
 605             log << "the following " + identicals.size.to_s + " sequences are identical:" + ld
 606             identicals.each { | identical |
 607               log << identical + ld
 608             }
 609             log << "ignoring the following " + removed.size.to_s + " redundant sequences:" + ld
 610             removed.each { | seq_name |
 611               log << seq_name + ld
 612             }
 613             Util.print_message( PRG_NAME, "will store " + msa.get_number_of_seqs.to_s + " non-redundant sequences" )
 614             log << "will store " + msa.get_number_of_seqs.to_s + " non-redundant sequences" + ld
 615           end
 616         end
 617
 618         io = MsaIO.new()
 619
 620         w = MsaWriter
 621
 622         if ( @pi_output )
 623           w = PhylipSequentialWriter.new()
 624           w.clean( @clean )
 625           w.set_max_name_length( @name_length )
 626           w.set_exception_if_name_too_long( @die_if_name_too_long )
 627         elsif( @fasta_output )
 628           w = FastaWriter.new()
 629           w.set_line_width( @width )
 630           if ( @rg )
 631             w.remove_gap_chars( true )
 632             Util.print_warning_message( PRG_NAME, "removing gap characters, the output is likely to become unaligned"  )
 633             log << "removing gap character, the output is likely to become unaligned" + ld
 634           end
 635           w.clean( @clean )
 636           if ( @name_length_set )
 637             w.set_max_name_length( @name_length )
 638             w.set_exception_if_name_too_long( @die_if_name_too_long )
 639           end
 640         elsif( @nexus_output )
 641           w = NexusWriter.new()
 642           w.clean( @clean )
 643           w.set_max_name_length( @name_length )
 644           w.set_exception_if_name_too_long( @die_if_name_too_long )
 645         end
 646
 647         begin
 648           io.write_to_file( msa, output, w )
 649         rescue Exception => e
 650           Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 651         end
 652
 653         Util.print_message( PRG_NAME, "Number of sequences in output        : " + msa.get_number_of_seqs.to_s )
 654         log << "Number of sequences in output        : " + msa.get_number_of_seqs.to_s + ld
 655
 656         begin
 657           f = File.open( output + LOG_SUFFIX, 'a' )
 658           f.print( log )
 659           f.close
 660         rescue Exception => e
 661           Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 662         end
 663
 664       end
 665       Util.print_message( PRG_NAME, "OK" )
 666       puts
 667     end
 668
 669     private
 670
 671     def sort( msa )
 672       names = Set.new
 673       for i in 0 ... msa.get_number_of_seqs
 674         name = msa.get_sequence( i ).get_name
 675         names.add( name )
 676       end
 677       sorted_ary = names.to_a.sort
 678       new_msa = Msa.new
 679       sorted_ary.each { | seq_name |
 680         seq = msa.get_sequence( msa.find_by_name( seq_name, true, false )[ 0 ] )
 681         new_msa.add_sequence( seq )
 682       }
 683       new_msa
 684     end
 685
 686     def set_fasta_input( fi = true )
 687       @fasta_input = fi
 688       @input_format_set = true
 689     end
 690
 691     def set_phylip_input( pi = true )
 692       @phylip_input = pi
 693       @input_format_set = true
 694     end
 695
 696     def set_name_length( i )
 697       @name_length = i
 698       @name_length_set = true
 699     end
 700
 701     def set_width( i )
 702       @width = i
 703     end
 704
 705     def set_fasta_output( fo = true )
 706       @fasta_output = fo
 707       @output_format_set = true
 708     end
 709
 710     def set_pi_output( pso = true )
 711       @pi_output = pso
 712       @output_format_set = true
 713     end
 714
 715     def set_nexus_output( nexus = true )
 716       @nexus_output = nexus
 717       @output_format_set = true
 718     end
 719
 720     def set_clean( c = true )
 721       @clean = c
 722     end
 723
 724     def set_remove_gap_columns( rgc = true )
 725       @rgc = rgc
 726     end
 727
 728     def set_remove_gap_only_columns( rgoc = true )
 729       @rgoc = rgoc
 730     end
 731
 732     def set_remove_gaps( rg = true )
 733       @rg = rg
 734     end
 735
 736     def set_remove_gap_ratio( rgr )
 737       @rgr = rgr
 738     end
 739
 740     def set_remove_seqs_gap_ratio( rsgr )
 741       @rsgr = rsgr
 742     end
 743
 744     def set_remove_seqs_min_non_gap_length( rsl )
 745       @rsl = rsl
 746     end
 747
 748     def set_remove_seqs( file )
 749       @seqs_name_file = file
 750       @remove_seqs    = true
 751       @keep_seqs      = false
 752     end
 753
 754     def set_keep_seqs( file )
 755       @seqs_name_file = file
 756       @keep_seqs      = true
 757       @remove_seqs    = false
 758     end
 759
 760     def set_trim( first, last )
 761       @trim            = true
 762       @first           = first
 763       @last            = last
 764     end
 765
 766     def set_remove_matching( remove )
 767       @remove_matching  = remove
 768     end
 769
 770     def set_keep_matching( keep )
 771       @keep_matching = keep
 772     end
 773
 774     def set_rem_red( rr )
 775       @rem_red = rr
 776     end
 777
 778     def set_split( s )
 779       if ( s > 0 )
 780         @split            = s
 781         @split_by_os      = false
 782         @clean            = false  # phylip only
 783         @rgc              = false
 784         @rgoc             = false
 785         @rg               = false  # fasta only
 786         @rgr              = -1
 787         @rsgr             = -1
 788         @rsl              = -1
 789         @seqs_name_file   = nil
 790         @remove_seqs      = false
 791         @keep_seqs        = false
 792         @trim             = false
 793         @first            = -1
 794         @last             = -1
 795       end
 796     end
 797
 798     def set_split_by_os()
 799       @split            = -1
 800       @split_by_os      = true
 801       @clean            = false  # phylip only
 802       @rgc              = false
 803       @rgoc             = false
 804       @rg               = false  # fasta only
 805       @rgr              = -1
 806       @rsgr             = -1
 807       @rsl              = -1
 808       @seqs_name_file   = nil
 809       @remove_seqs      = false
 810       @keep_seqs        = false
 811       @trim             = false
 812       @first            = -1
 813       @last             = -1
 814       @window           = false
 815     end
 816
 817     def set_window()
 818       @split            = -1
 819       @split_by_os      = false
 820       @rgc              = false
 821       @rgoc             = false
 822       @rg               = false  # fasta only
 823       @rgr              = -1
 824       @rsgr             = -1
 825       @rsl              = -1
 826       @seqs_name_file   = nil
 827       @remove_seqs      = false
 828       @keep_seqs        = false
 829       @trim             = false
 830       @first            = -1
 831       @last             = -1
 832       @window           = true
 833     end
 834
 835     def analyze_command_line( cla )
 836       if ( cla.is_option_set?( INPUT_TYPE_OPTION ) )
 837         begin
 838           type = cla.get_option_value( INPUT_TYPE_OPTION )
 839           if ( type == "p" )
 840             set_phylip_input( true )
 841             set_fasta_input( false )
 842           elsif ( type == "f" )
 843             set_fasta_input( true )
 844             set_phylip_input( false )
 845           end
 846         rescue ArgumentError => e
 847           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 848         end
 849       end
 850       if ( cla.is_option_set?( OUTPUT_TYPE_OPTION ) )
 851         begin
 852           type = cla.get_option_value( OUTPUT_TYPE_OPTION )
 853           if ( type == "p" )
 854             set_pi_output( true )
 855             set_fasta_output( false )
 856             set_nexus_output( false )
 857           elsif ( type == "f" )
 858             set_pi_output( false )
 859             set_fasta_output( true )
 860             set_nexus_output( false )
 861           elsif ( type == "n" )
 862             set_pi_output( false )
 863             set_fasta_output( false )
 864             set_nexus_output( true )
 865           end
 866         rescue ArgumentError => e
 867           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 868         end
 869       end
 870       if ( cla.is_option_set?( MAXIMAL_NAME_LENGTH_OPTION ) )
 871         begin
 872           l = cla.get_option_value_as_int( MAXIMAL_NAME_LENGTH_OPTION )
 873           set_name_length( l )
 874         rescue ArgumentError => e
 875           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 876         end
 877       end
 878       if ( cla.is_option_set?( WIDTH_OPTION ) )
 879         begin
 880           w = cla.get_option_value_as_int( WIDTH_OPTION )
 881           set_width( w )
 882         rescue ArgumentError => e
 883           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 884         end
 885       end
 886       if ( cla.is_option_set?( CLEAN_UP_SEQ_OPTION ) )
 887         set_clean( true )
 888       end
 889       if ( cla.is_option_set?( REMOVE_GAP_COLUMNS_OPTION ) )
 890         set_remove_gap_columns( true )
 891       end
 892       if ( cla.is_option_set?( REM_RED_OPTION ) )
 893         set_rem_red( true )
 894       end
 895       if ( cla.is_option_set?( REMOVE_GAP_ONLY_COLUMNS ) )
 896         set_remove_gap_only_columns( true )
 897       end
 898       if ( cla.is_option_set?( REMOVE_ALL_GAP_CHARACTERS_OPTION ) )
 899         set_remove_gaps( true )
 900       end
 901       if ( cla.is_option_set?( REMOVE_COLUMNS_GAP_RATIO_OPTION ) )
 902         begin
 903           f = cla.get_option_value_as_float( REMOVE_COLUMNS_GAP_RATIO_OPTION )
 904           set_remove_gap_ratio( f )
 905         rescue ArgumentError => e
 906           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 907         end
 908       end
 909       if ( cla.is_option_set?( REMOVE_ALL_SEQUENCES_LISTED_OPTION ) )
 910         begin
 911           s = cla.get_option_value( REMOVE_ALL_SEQUENCES_LISTED_OPTION )
 912           set_remove_seqs( s )
 913         rescue ArgumentError => e
 914           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 915         end
 916       end
 917       if ( cla.is_option_set?( KEEP_ONLY_SEQUENCES_LISTED_OPTION ) )
 918         begin
 919           s = cla.get_option_value( KEEP_ONLY_SEQUENCES_LISTED_OPTION )
 920           set_keep_seqs( s )
 921         rescue ArgumentError => e
 922           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 923         end
 924       end
 925       if ( cla.is_option_set?( TRIM_OPTION ) )
 926         begin
 927           s = cla.get_option_value( TRIM_OPTION )
 928           if ( s =~ /(\d+)-(\d+)/ )
 929             set_trim( $1.to_i(), $2.to_i() )
 930           else
 931             puts( "illegal argument" )
 932             print_help
 933             exit( -1 )
 934           end
 935         rescue ArgumentError => e
 936           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 937         end
 938       end
 939       if ( cla.is_option_set?( SLIDING_EXTRACTION_OPTION ) )
 940         begin
 941           s = cla.get_option_value( SLIDING_EXTRACTION_OPTION )
 942           if ( s =~ /(\d+)\/(\d+)/ )
 943             set_window
 944             @window = true
 945             @step = $1.to_i()
 946             @size = $2.to_i()
 947           else
 948             puts( "illegal argument" )
 949             print_help
 950             exit( -1 )
 951           end
 952           if (@step <= 0) || (@size <= 0)
 953             puts( "illegal argument" )
 954             print_help
 955             exit( -1 )
 956           end
 957         rescue ArgumentError => e
 958           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 959         end
 960       end
 961
 962       if ( cla.is_option_set?( REMOVE_SEQS_GAP_RATIO_OPTION ) )
 963         begin
 964           f = cla.get_option_value_as_float( REMOVE_SEQS_GAP_RATIO_OPTION )
 965           set_remove_seqs_gap_ratio( f )
 966         rescue ArgumentError => e
 967           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 968         end
 969       end
 970       if ( cla.is_option_set?( REMOVE_SEQS_NON_GAP_LENGTH_OPTION ) )
 971         begin
 972           f = cla.get_option_value_as_int( REMOVE_SEQS_NON_GAP_LENGTH_OPTION )
 973           set_remove_seqs_min_non_gap_length( f )
 974         rescue ArgumentError => e
 975           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 976         end
 977       end
 978       if cla.is_option_set?( SPLIT_BY_OS )
 979         begin
 980           set_split_by_os()
 981         rescue ArgumentError => e
 982           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 983         end
 984       elsif ( cla.is_option_set?( SPLIT ) )
 985         begin
 986           s = cla.get_option_value_as_int( SPLIT )
 987           set_split( s )
 988         rescue ArgumentError => e
 989           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 990         end
 991       end
 992       if ( cla.is_option_set?( REMOVE_MATCHING_SEQUENCES_OPTION ) )
 993         begin
 994           s = cla.get_option_value( REMOVE_MATCHING_SEQUENCES_OPTION )
 995           set_remove_matching( s )
 996         rescue ArgumentError => e
 997           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 998         end
 999       end
1000       if ( cla.is_option_set?( KEEP_MATCHING_SEQUENCES_OPTION ) )
1001         begin
1002           s = cla.get_option_value( KEEP_MATCHING_SEQUENCES_OPTION )
1003           set_keep_matching( s )
1004         rescue ArgumentError => e
1005           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
1006         end
1007       end
1008       if ( cla.is_option_set?( DIE_IF_NAME_TOO_LONG ) )
1009         @die_if_name_too_long = true
1010       end
1011
1012     end
1013
1014     def print_help()
1015       puts()
1016       puts( "Usage:" )
1017       puts()
1018       puts( "  " + PRG_NAME + ".rb [options] <input alignment> <output>" )
1019       puts()
1020       puts( "  options: -" + INPUT_TYPE_OPTION + "=<input type>: f for fasta (default), p for phylip/selex type" )
1021       puts( "           -" + OUTPUT_TYPE_OPTION + "=<output type>: f for fasta (default), n for nexus, p for phylip sequential" )
1022       puts( "           -" + MAXIMAL_NAME_LENGTH_OPTION + "=<n>: n=maximal name length (default for phylip 10, for fasta: unlimited )" )
1023       puts( "           -" + DIE_IF_NAME_TOO_LONG + ": die if sequence name too long" )
1024       puts( "           -" + WIDTH_OPTION + "=<n>: n=width (fasta output only, default is 60)" )
1025       puts( "           -" + CLEAN_UP_SEQ_OPTION + ": clean up sequences" )
1026       puts( "           -" + REMOVE_GAP_COLUMNS_OPTION + ": remove gap columns" )
1027       puts( "           -" + REMOVE_GAP_ONLY_COLUMNS + ": remove gap-only columns" )
1028       puts( "           -" + REMOVE_COLUMNS_GAP_RATIO_OPTION + "=<n>: remove columns for which ( seqs with gap / number of sequences > n )" )
1029       puts( "           -" + REMOVE_ALL_GAP_CHARACTERS_OPTION + ": remove all gap characters (destroys alignment, fasta output only)" )
1030       puts( "           -" + REMOVE_ALL_SEQUENCES_LISTED_OPTION + "=<file>: remove all sequences listed in file" )
1031       puts( "           -" + KEEP_ONLY_SEQUENCES_LISTED_OPTION + "=<file>: keep only sequences listed in file" )
1032       puts( "           -" + TRIM_OPTION + "=<first>-<last>: remove columns before first and after last" )
1033       puts( "           -" + REMOVE_SEQS_GAP_RATIO_OPTION + "=<n>: remove sequences for which the gap ratio > n (after column operations)" )
1034       puts( "           -" + REMOVE_SEQS_NON_GAP_LENGTH_OPTION + "=<n> remove sequences with less than n non-gap characters (after column operations)" )
1035       puts( "           -" + REMOVE_MATCHING_SEQUENCES_OPTION + "=<s> remove all sequences with names containing s" )
1036       puts( "           -" + KEEP_MATCHING_SEQUENCES_OPTION + "=<s> keep only sequences with names containing s" )
1037       puts( "           -" + SPLIT + "=<n> split a fasta file into n files of equal number of sequences (expect for " )
1038       puts( "            last one), cannot be used with other options" )
1039       puts( "           -" + SLIDING_EXTRACTION_OPTION + "=<step>/<window size>: sliding window extraction, cannot be used with other options" )
1040       puts( "           -" + REM_RED_OPTION + ": remove redundant sequences" )
1041       puts()
1042     end
1043
1044   end # class MsaProcessor
1045
1046 end # module Evoruby