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