forester/ruby/evoruby/lib/evo/msa/msa.rb

   1 #
   2 # = lib/evo/msa/msa.rb - Msa class
   3 #
   4 # Copyright::  Copyright (C) 2006-2007 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: msa.rb,v 1.11 2009/01/03 00:42:08 cmzmasek Exp $
   8 #
   9
  10
  11 require 'lib/evo/util/constants'
  12 require 'lib/evo/util/util'
  13 require 'lib/evo/sequence/sequence'
  14
  15 module Evoruby
  16
  17   class Msa
  18
  19     def initialize()
  20       @sequences = Array.new
  21       @identical_seqs_detected = Array.new
  22     end
  23
  24
  25     def add_sequence( sequence )
  26       @sequences.push( sequence )
  27     end
  28
  29     def add( name, molecular_sequence_str )
  30       add_sequence( Sequence.new( name, molecular_sequence_str ) )
  31     end
  32
  33     def get_sequence( index )
  34       if ( index < 0 || index > get_number_of_seqs() - 1 )
  35         error_msg = "attempt to get sequence " <<
  36          index.to_s << " in alignment of " << get_number_of_seqs().to_s <<
  37          " sequences"
  38         raise ArgumentError, error_msg
  39       end
  40       return @sequences[ index ]
  41     end
  42
  43     def remove_sequence!( index )
  44       if ( index < 0 || index > get_number_of_seqs() - 1 )
  45         error_msg = "attempt to remove sequence " <<
  46          index.to_s << " in alignment of " << get_number_of_seqs().to_s <<
  47          " sequences"
  48         raise ArgumentError, error_msg
  49       end
  50       @sequences.delete_at( index )
  51     end
  52
  53     def get_identical_seqs_detected
  54       @identical_seqs_detected
  55     end
  56
  57
  58     def is_aligned()
  59       if ( get_number_of_seqs < 1 )
  60         return false
  61       else
  62         l = @sequences[ 0 ].get_length()
  63         for i in 0 ... get_number_of_seqs()
  64           if ( get_sequence( i ).get_length() != l )
  65             return false
  66           end
  67         end
  68       end
  69       return true
  70     end
  71
  72     def find_by_name( name, case_sensitive, partial_match )
  73       indices = Array.new()
  74       for i in 0 ... get_number_of_seqs()
  75         current_name = get_sequence( i ).get_name()
  76         if !case_sensitive
  77           current_name = current_name.downcase
  78           name = name.downcase
  79         end
  80         if current_name == name ||
  81            ( partial_match && current_name.include?( name ) )
  82           indices.push( i )
  83         end
  84       end
  85       indices
  86     end
  87
  88     def find_by_name_start( name, case_sensitive )
  89       indices = Array.new()
  90       for i in 0 ... get_number_of_seqs()
  91         get_sequence( i ).get_name() =~ /^\s*(\S+)/
  92         current_name = $1
  93         if !case_sensitive
  94           current_name = current_name.downcase
  95           name = name.downcase
  96         end
  97         if  ( current_name == name )
  98           indices.push( i )
  99         end
 100       end
 101       indices
 102     end
 103
 104     def has?( name, case_sensitive = true, partial_match = false )
 105       for i in 0 ... get_number_of_seqs()
 106         current_name = get_sequence( i ).get_name()
 107         if !case_sensitive
 108           current_name = current_name.downcase
 109           name = name.downcase
 110         end
 111         if current_name == name ||
 112            ( partial_match && current_name.include?( name ) )
 113           return true
 114         end
 115       end
 116       false
 117     end
 118
 119     # throws ArgumentError
 120     def get_by_name( name, case_sensitive = true, partial_match = false )
 121       indices = find_by_name( name, case_sensitive, partial_match )
 122       if ( indices.length > 1 )
 123         error_msg = "\"" + name + "\" not unique"
 124         raise ArgumentError, error_msg
 125       elsif ( indices.length < 1 )
 126         error_msg = "\"" + name + "\" not found"
 127         raise ArgumentError, error_msg
 128       end
 129       get_sequence( indices[ 0 ] )
 130     end
 131
 132     # throws ArgumentError
 133     def get_by_name_start( name, case_sensitive = true )
 134       indices = find_by_name_start( name, case_sensitive )
 135       if ( indices.length > 1 )
 136         error_msg = "\"" + name + "\" not unique"
 137         raise ArgumentError, error_msg
 138       elsif ( indices.length < 1 )
 139         error_msg = "\"" + name + "\" not found"
 140         raise ArgumentError, error_msg
 141       end
 142       get_sequence( indices[ 0 ] )
 143     end
 144
 145
 146     def get_sub_alignment( seq_numbers )
 147       msa = Msa.new()
 148       for i in 0 ... seq_numbers.length()
 149         msa.add_sequence( get_sequence( seq_numbers[ i ] ).copy() )
 150       end
 151       return msa
 152     end
 153
 154     def get_number_of_seqs()
 155       @sequences.length
 156     end
 157
 158     def get_length
 159       if !is_aligned()
 160         error_msg = "attempt to get length of unaligned msa"
 161         raise StandardError, error_msg, caller
 162       end
 163       if get_number_of_seqs() < 1
 164         -1
 165       else
 166         @sequences[ 0 ].get_length()
 167       end
 168     end
 169
 170     def to_str
 171       to_fasta
 172     end
 173
 174     def to_fasta
 175       s = String.new
 176       for i in 0...get_number_of_seqs
 177         s += @sequences[ i ].to_fasta + Constants::LINE_DELIMITER
 178       end
 179       s
 180     end
 181
 182
 183     def print_overlap_diagram( min_overlap = 1, io = STDOUT, max_name_length = 10 )
 184       if ( !is_aligned() )
 185         error_msg = "attempt to get overlap diagram of unaligned msa"
 186         raise StandardError, error_msg, caller
 187       end
 188       for i in 0 ... get_number_of_seqs()
 189         io.print( Util.normalize_seq_name( get_sequence( i ).get_name(), max_name_length ) )
 190         for j in 0 ... get_number_of_seqs()
 191           if i == j
 192             io.print( " " )
 193           else
 194             if overlap?( i, j, min_overlap )
 195               io.print( "+" )
 196             else
 197               io.print( "-" )
 198             end
 199           end
 200         end
 201         io.print( Evoruby::Constants::LINE_DELIMITER )
 202       end
 203     end
 204
 205     #returns array of Msa with an overlap of min_overlap
 206     def split_into_overlapping_msa( min_overlap = 1 )
 207       if ( !is_aligned() )
 208         error_msg = "attempt to split into overlapping msas of unaligned msa"
 209         raise StandardError, error_msg, caller
 210       end
 211       msas = Array.new()
 212       bins = get_overlaps( min_overlap )
 213       for i in 0 ... bins.length
 214         msas.push( get_sub_alignment( bins[ i ] ) )
 215       end
 216       msas
 217     end
 218
 219     def overlap?( index_1, index_2, min_overlap = 1 )
 220       seq_1 = get_sequence( index_1 )
 221       seq_2 = get_sequence( index_2 )
 222       overlap_count = 0
 223       for i in 0...seq_1.get_length()
 224         if !Util.is_aa_gap_character?( seq_1.get_character_code( i ) ) &&
 225            !Util.is_aa_gap_character?( seq_2.get_character_code( i ) )
 226           overlap_count += 1
 227           if overlap_count >= min_overlap
 228             return true
 229           end
 230         end
 231       end
 232       return false
 233     end
 234
 235     def calculate_overlap( index_1, index_2 )
 236       seq_1 = get_sequence( index_1 )
 237       seq_2 = get_sequence( index_2 )
 238       overlap_count = 0
 239       for i in 0...seq_1.get_length
 240         if !Util.is_aa_gap_character?( seq_1.get_character_code( i ) ) &&
 241            !Util.is_aa_gap_character?( seq_2.get_character_code( i ) )
 242           overlap_count += 1
 243         end
 244       end
 245       overlap_count
 246     end
 247
 248     def calculate_identities( index_1, index_2 )
 249       seq_1 = get_sequence( index_1 )
 250       seq_2 = get_sequence( index_2 )
 251       identities_count = 0
 252       for i in 0...seq_1.get_length
 253         if !Util.is_aa_gap_character?( seq_1.get_character_code( i ) ) &&
 254            !Util.is_aa_gap_character?( seq_2.get_character_code( i ) ) &&
 255            seq_1.get_character_code( i ) != 63 &&
 256            ( seq_1.get_residue( i ).downcase() ==
 257              seq_2.get_residue( i ).downcase() )
 258           identities_count += 1
 259         end
 260       end
 261       identities_count
 262     end
 263
 264     def remove_gap_only_columns!()
 265       remove_columns!( get_gap_only_columns() )
 266     end
 267
 268     def remove_gap_columns!()
 269       remove_columns!( get_gap_columns() )
 270     end
 271
 272     # removes columns for which seqs with gap / number of sequences > gap_ratio
 273     def remove_gap_columns_w_gap_ratio!( gap_ratio )
 274       remove_columns!( get_gap_columns_w_gap_ratio( gap_ratio ) )
 275     end
 276
 277
 278     def remove_sequences_by_gap_ratio!( gap_ratio )
 279       if ( !is_aligned() )
 280         error_msg = "attempt to remove sequences by gap ratio on unaligned msa"
 281         raise StandardError, error_msg, caller
 282       end
 283       n = get_number_of_seqs
 284       removed = Array.new
 285       for s in 0 ... n
 286         if ( get_sequence( ( n - 1 ) - s  ).get_gap_ratio() > gap_ratio )
 287           if ( Evoruby::Constants::VERBOSE )
 288             puts( "removed: " + get_sequence( ( n - 1 ) - s  ).get_name )
 289           end
 290           removed << get_sequence( ( n - 1 ) - s  ).get_name
 291           remove_sequence!( ( n - 1 ) - s  )
 292         end
 293       end
 294       removed
 295     end
 296
 297
 298     def remove_redundant_sequences!( consider_taxonomy = false, verbose = false )
 299       n = get_number_of_seqs
 300       removed = Array.new
 301       to_be_removed = Set.new
 302       @identical_seqs_detected = Array.new
 303       for i in 0 ... ( n - 1 )
 304         for j in ( i + 1 ) ... n
 305           if !to_be_removed.include?( i ) && !to_be_removed.include?( j )
 306             if  !consider_taxonomy ||
 307                ( ( get_sequence( i ).get_taxonomy == nil && get_sequence( j ).get_taxonomy == nil ) ||
 308                  ( get_sequence( i ).get_taxonomy == get_sequence( j ).get_taxonomy ) )
 309               if Util.clean_seq_str( get_sequence( i ).get_sequence_as_string ) ==
 310                  Util.clean_seq_str( get_sequence( j ).get_sequence_as_string )
 311                 to_be_removed.add( j )
 312
 313                 tax_i = ""
 314                 tax_j = ""
 315                 if get_sequence( i ).get_taxonomy != nil
 316                   tax_i = get_sequence( i ).get_taxonomy.get_name
 317                 end
 318                 if get_sequence( j ).get_taxonomy != nil
 319                   tax_j = get_sequence( j ).get_taxonomy.get_name
 320                 end
 321                 identical_pair = get_sequence( i ).get_name + " [#{tax_i}] == " + get_sequence( j ).get_name + " [#{tax_j}]"
 322                 @identical_seqs_detected.push( identical_pair )
 323                 if verbose
 324                   puts identical_pair
 325                 end
 326               end
 327             end
 328           end
 329         end
 330       end
 331
 332       to_be_removed_ary = to_be_removed.to_a.sort.reverse
 333
 334       to_be_removed_ary.each { | index |
 335         removed.push( get_sequence( index ).get_name )
 336         remove_sequence!( index )
 337       }
 338       removed
 339     end
 340
 341
 342     def remove_sequences_by_non_gap_length!( min_non_gap_length )
 343       if ( !is_aligned() )
 344         error_msg = "attempt to remove sequences by non gap length on unaligned msa"
 345         raise StandardError, error_msg, caller
 346       end
 347       n = get_number_of_seqs
 348       l = get_length
 349       removed = Array.new
 350       for s in 0 ... n
 351         if ( ( l - get_sequence( ( n - 1 ) - s ).get_gap_length ) < min_non_gap_length )
 352           if ( Evoruby::Constants::VERBOSE )
 353             puts( "removed: " + get_sequence( ( n - 1 ) - s  ).get_name )
 354           end
 355           removed << get_sequence( ( n - 1 ) - s  ).get_name
 356           remove_sequence!( ( n - 1 ) - s )
 357         end
 358       end
 359       removed
 360     end
 361
 362     def trim!( first, last )
 363       cols = Array.new()
 364       for i in 0 ... get_length()
 365         if ( i < first || i > last )
 366           cols.push( i )
 367         end
 368       end
 369       remove_columns!( cols )
 370     end
 371
 372     def get_gap_only_columns()
 373       if ( !is_aligned() )
 374         error_msg = "attempt to get gap only columns of unaligned msa"
 375         raise StandardError, error_msg, caller
 376       end
 377       cols = Array.new()
 378       for c in 0 ... get_length
 379         nogap_char_found = false
 380         for s in 0 ... get_number_of_seqs
 381           unless Util.is_aa_gap_character?( get_sequence( s ).get_character_code( c ) )
 382             nogap_char_found = true
 383             break
 384           end
 385         end
 386         unless nogap_char_found
 387           cols.push( c )
 388         end
 389       end
 390       return cols
 391     end
 392
 393     def calculate_gap_proportion()
 394       if ( !is_aligned() )
 395         error_msg = "attempt to get gap only columns of unaligned msa"
 396         raise StandardError, error_msg, caller
 397       end
 398       total_sum = 0.0
 399       gap_sum = 0.0
 400       for c in 0 ... get_length
 401         for s in 0 ... get_number_of_seqs
 402           total_sum = total_sum + 1
 403           if Util.is_aa_gap_character?( get_sequence( s ).get_character_code( c ) )
 404             gap_sum = gap_sum  + 1
 405           end
 406         end
 407
 408       end
 409       return gap_sum / total_sum
 410     end
 411
 412     def get_gap_columns()
 413       if ( !is_aligned() )
 414         error_msg = "attempt to get gap columns of unaligned msa"
 415         raise StandardError, error_msg, caller
 416       end
 417       cols = Array.new()
 418       for c in 0 ... get_length
 419         gap_char_found = false
 420         for s in 0 ... get_number_of_seqs
 421           if Util.is_aa_gap_character?( get_sequence( s ).get_character_code( c ) )
 422             gap_char_found = true
 423             break
 424           end
 425         end
 426         if gap_char_found
 427           cols.push( c )
 428         end
 429       end
 430       return cols
 431     end
 432
 433     # gap_ratio = seqs with gap / number of sequences
 434     # returns column indices for which seqs with gap / number of sequences > gap_ratio
 435     def get_gap_columns_w_gap_ratio( gap_ratio )
 436       if ( !is_aligned() )
 437         error_msg = "attempt to get gap columns with gap_ratio of unaligned msa"
 438         raise StandardError, error_msg, caller
 439       end
 440       if ( gap_ratio < 0 || gap_ratio > 1 )
 441         error_msg = "gap ratio must be between 0 and 1 inclusive"
 442         raise ArgumentError, error_msg, caller
 443       end
 444       cols = Array.new()
 445       for c in 0 ... get_length
 446         gap_chars_found = 0
 447         for s in 0 ... get_number_of_seqs
 448           if Util.is_aa_gap_character?( get_sequence( s ).get_character_code( c ) )
 449             gap_chars_found += 1
 450           end
 451         end
 452         if ( ( gap_chars_found.to_f / get_number_of_seqs ) > gap_ratio )
 453           cols.push( c )
 454         end
 455       end
 456       return cols
 457     end
 458
 459
 460     # Split an alignment into n alignemnts of equal size, except last one
 461     def split( n, verbose = false )
 462       if ( n < 2 || n > get_number_of_seqs )
 463         error_msg = "attempt to split into less than two or more than the number of sequences"
 464         raise StandardError, error_msg, caller
 465       end
 466       msas = Array.new()
 467       r = get_number_of_seqs % n
 468       x = get_number_of_seqs / n
 469       for i in 0 ... n
 470         msa = Msa.new()
 471         s = 0
 472
 473         if ( ( r > 0 ) && ( i == ( n - 1 ) ) )
 474           y = x + r
 475           if ( verbose )
 476             puts( i.to_s + ": " + y.to_s )
 477           end
 478           for j in 0 ... y
 479             msa.add_sequence( get_sequence( ( i * x ) + j ) )
 480           end
 481         else
 482           if ( verbose )
 483             puts( i.to_s + ": " + x.to_s )
 484           end
 485           for j in 0 ... x
 486             msa.add_sequence( get_sequence( ( i * x ) + j ) )
 487           end
 488         end
 489         msas.push( msa )
 490       end
 491       msas
 492     end
 493
 494
 495     private
 496
 497     def get_overlaps( min_overlap = 1 )
 498       if ( !is_aligned() )
 499         error_msg = "attempt to get overlaps of unaligned msa"
 500         raise StandardError, error_msg, caller
 501       end
 502       bins = Array.new()
 503       for i in 0 ... get_number_of_seqs()
 504         found_bin = false
 505         for j in 0 ... bins.length
 506           current_bin = bins[ j ]
 507           # does seq i overlap with all seqs in current_bin?
 508           all_overlap = true
 509           for z in 0 ... current_bin.length
 510             unless overlap?( i, current_bin[ z ], min_overlap )
 511               all_overlap = false
 512               break
 513             end
 514           end
 515           if all_overlap
 516             current_bin.push( i )
 517             found_bin = true
 518           end
 519         end
 520         unless found_bin
 521           new_bin = Array.new()
 522           new_bin.push( i )
 523           bins.push( new_bin )
 524         end
 525       end
 526       return bins
 527     end
 528
 529     def remove_columns!( cols )
 530       if ( !is_aligned() )
 531         error_msg = "attempt to remove columns of unaligned msa"
 532         raise StandardError, error_msg, caller
 533       end
 534       cols.reverse!()
 535       for c in 0 ... cols.length()
 536         col = cols[ c ]
 537         for s in 0 ... get_number_of_seqs()
 538           get_sequence( s ).delete_residue!( col )
 539         end
 540       end
 541       return self
 542     end
 543
 544
 545   end # class Msa
 546
 547 end # module Evoruby
 548