forester/ruby/evoruby/exe/select_same_gn.rb

   1 #!/usr/local/bin/ruby -w
   2
   3
   4 require 'lib/evo/sequence/sequence'
   5 require 'lib/evo/msa/msa'
   6 require 'lib/evo/msa/msa_factory'
   7 require 'lib/evo/io/writer/fasta_writer'
   8 require 'lib/evo/io/parser/fasta_parser'
   9
  10
  11 module Evoruby
  12
  13   if ARGV.length != 1
  14     puts "usage: select_same_gn.rb <fasta formatted multiple sequence file>"
  15     exit
  16   end
  17
  18   input = ARGV[ 0 ]
  19   f = MsaFactory.new()
  20
  21   IGNORE_SEQS_LACKING_GN = false
  22   IGNORE_FRAGMENTS = false
  23   IGNORE_SPECIES = true
  24
  25   msa = nil
  26
  27   begin
  28     msa = f.create_msa_from_file( input, FastaParser.new() )
  29   rescue Exception => e
  30     puts "error: " + e.to_s
  31     exit
  32   end
  33
  34   outbase = input.sub( /\..+/, "" )
  35
  36   outfile = File.open(outbase + "_log.txt" , "w")
  37
  38   all_names = Set.new
  39   all_seqs_per_species = Hash.new
  40   all_msa_per_species = Hash.new
  41   gn_to_seqs = Hash.new
  42   unique_genes_msa = Msa.new
  43   longest_non_unique_genes_msa = Msa.new
  44   gn_re = /GN=(\S+)/
  45   fragment_re = /fragment/i
  46   species_re = /\[([A-Z0-9]{3,5})\]$/
  47
  48   frag_counter = 0
  49   no_gn_counter = 0
  50   same_seq_counter = 0
  51
  52   for i in 0 ... msa.get_number_of_seqs()
  53     seq = msa.get_sequence( i )
  54     name = seq.get_name
  55     if all_names.include?( name )
  56       puts "error: sequence name \"" + name + "\" is not unique (#" + i.to_s + ")"
  57       exit
  58     else
  59       all_names << name
  60     end
  61
  62     if IGNORE_FRAGMENTS && fragment_re.match( name )
  63       outfile.puts("ignored because fragment: " + name)
  64       frag_counter += 1
  65       next
  66     end
  67
  68     species = nil
  69     if IGNORE_SPECIES || species_re.match( name )
  70       unless IGNORE_SPECIES
  71         species = species_re.match( name )[ 1 ]
  72       else
  73         species = "XXXXX"
  74       end
  75
  76       unless all_seqs_per_species.has_key?( species )
  77         all_seqs_per_species[ species ] = Set.new
  78       end
  79       all_seqs = all_seqs_per_species[ species ]
  80       mol_seq = seq.get_sequence_as_string.upcase
  81       if all_seqs.include?( mol_seq )
  82         outfile.puts("ignored because identical sequence in same species: " + name )
  83
  84         same_seq_counter += 1
  85         next
  86       else
  87         all_seqs << mol_seq
  88       end
  89     else
  90       puts "error: no species for: " + name
  91       exit
  92     end
  93
  94     gn_match = gn_re.match( name )
  95     if IGNORE_SEQS_LACKING_GN
  96       unless gn_match
  97         outfile.puts( "ignored because no GN=: " + name )
  98         no_gn_counter += 1
  99         next
 100       end
 101     else
 102       unless gn_match
 103         outfile.puts( "no GN=: " + name )
 104       end
 105     end
 106
 107     gn =nil
 108     if gn_match
 109       gn = gn_match[1] + "_" + species
 110     else
 111       if IGNORE_SEQS_LACKING_GN
 112         puts "cannot be"
 113         exit
 114       end
 115       gn = name
 116     end
 117
 118     unless gn_to_seqs.has_key?(gn)
 119       gn_to_seqs[gn] = Msa.new
 120     end
 121     gn_to_seqs[gn].add_sequence(seq)
 122   end
 123
 124   if IGNORE_FRAGMENTS
 125     outfile.puts( "Sequences ignored because \"fragment\" in desc                : " + frag_counter.to_s )
 126   end
 127
 128   if IGNORE_SEQS_LACKING_GN
 129     outfile.puts( "Sequences ignored because no \"GN=\" in desc                  : " + no_gn_counter.to_s )
 130   end
 131   outfile.puts( "Sequences ignored because identical sequence in same species: " + same_seq_counter.to_s )
 132   outfile.puts
 133   outfile.puts
 134
 135   counter = 1
 136   gn_to_seqs.each_pair do |gene,seqs|
 137     seq = nil
 138     if seqs.get_number_of_seqs > 1
 139       outfile.puts( counter.to_s + ": " + gene )
 140       outfile.puts( seqs.to_fasta )
 141       outfile.puts
 142       outfile.puts
 143       counter += 1
 144       longest = 0
 145       longest_seq = nil
 146       for j in 0 ... seqs.get_number_of_seqs()
 147         current = seqs.get_sequence( j )
 148         if current.get_length > longest
 149           longest =  current.get_length
 150           longest_seq = current
 151         end
 152       end
 153       seq = longest_seq
 154       longest_non_unique_genes_msa.add_sequence( seq )
 155     else
 156       seq = seqs.get_sequence( 0 )
 157       unique_genes_msa.add_sequence( seq )
 158     end
 159
 160
 161     unless IGNORE_SPECIES
 162       species = species_re.match( seq.get_name )[ 1 ]
 163     else
 164       species = "XXXXX"
 165     end
 166     unless all_msa_per_species.has_key?( species )
 167       all_msa_per_species[ species ] = Msa.new
 168     end
 169     all_msa_per_species[ species ].add_sequence( seq )
 170
 171   end
 172
 173   outfile.close
 174
 175   w = FastaWriter.new
 176   w.write( unique_genes_msa, outbase + "_seqs_from_unique_genes.fasta" )
 177   w.write( longest_non_unique_genes_msa, outbase + "_longest_seqs_from_nonunique_genes.fasta" )
 178
 179   all_msa_per_species.each_pair do |s,m|
 180     w = FastaWriter.new
 181     w.write( m, outbase + "_"  + s + ".fasta" )
 182   end
 183
 184 end