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   input = ARGV[ 0 ]
  14   f = MsaFactory.new()
  15
  16   IGNORE_SEQS_LACKING_GN = false
  17
  18   msa = nil
  19
  20   begin
  21     msa = f.create_msa_from_file( input, FastaParser.new() )
  22   rescue Exception => e
  23     puts "error: " + e.to_s
  24     exit
  25   end
  26
  27   all_names = Set.new
  28   all_seqs_per_species = Hash.new
  29   gn_to_seqs = Hash.new
  30   unique_genes_msa = Msa.new
  31   longest_non_unique_genes_msa = Msa.new
  32   gn_re = /GN=(\S+)/
  33   fragment_re = /fragment/i
  34   species_re = /\[([A-Z0-9]{3,6})\]$/
  35
  36   frag_counter = 0
  37   no_gn_counter = 0
  38   same_seq_counter = 0
  39
  40   for i in 0 ... msa.get_number_of_seqs()
  41     seq = msa.get_sequence( i )
  42     name = seq.get_name
  43     if all_names.include?( name )
  44       puts "error: sequence name \"" + name + "\" is not unique (#" + i.to_s + ")"
  45       exit
  46     else
  47       all_names << name
  48     end
  49
  50     if fragment_re.match( name )
  51       puts "ignored because fragment: " + name
  52       frag_counter += 1
  53       next
  54     end
  55
  56     species = nil
  57     if species_re.match( name )
  58       s_match = species_re.match( name )
  59       species = s_match[1]
  60
  61       unless all_seqs_per_species.include?( species )
  62         all_seqs_per_species[ species ] = Set.new
  63       end
  64       all_seqs = all_seqs_per_species[ species ]
  65       mol_seq = seq.get_sequence_as_string.upcase
  66       if all_seqs.include?( mol_seq )
  67         puts "ignored because identical sequence in same species: " + name
  68         same_seq_counter += 1
  69         next
  70       else
  71         all_seqs << mol_seq
  72       end
  73     else
  74       puts "error: no species for: " + name
  75       exit
  76     end
  77
  78     gn_match = gn_re.match( name )
  79     if IGNORE_SEQS_LACKING_GN
  80       unless gn_match
  81         puts "ignored because no GN=: " + name
  82         no_gn_counter += 1
  83         next
  84       end
  85     else
  86       unless gn_match
  87         puts "no GN=: " + name
  88       end
  89     end
  90
  91     gn =nil
  92     if gn_match
  93       gn = gn_match[1] + "_" + species
  94     else
  95       if IGNORE_SEQS_LACKING_GN
  96         puts "cannot be"
  97         exit
  98       end
  99       gn = name
 100     end
 101
 102     unless gn_to_seqs.has_key?(gn)
 103       gn_to_seqs[gn] = Msa.new
 104     end
 105     gn_to_seqs[gn].add_sequence(seq)
 106   end
 107
 108   puts "Sequences ignored because \"fragment\" in desc                : " + frag_counter.to_s
 109   if IGNORE_SEQS_LACKING_GN
 110     puts "Sequences ignored because no \"GN=\" in desc                  : " + no_gn_counter.to_s
 111   end
 112   puts "Sequences ignored because identical sequence in same species: " + same_seq_counter.to_s
 113   puts
 114   puts
 115
 116   counter = 1
 117   gn_to_seqs.each_pair do |gene,seqs|
 118     if seqs.get_number_of_seqs > 1
 119       puts counter.to_s + ": " + gene
 120       puts seqs.to_fasta
 121       puts
 122       puts
 123       counter += 1
 124       longest = 0
 125       longest_seq = nil
 126       for j in 0 ... seqs.get_number_of_seqs()
 127         current = seqs.get_sequence( j )
 128         if current.get_length > longest
 129           longest =  current.get_length
 130           longest_seq = current
 131         end
 132       end
 133       longest_non_unique_genes_msa.add_sequence(longest_seq)
 134     else
 135       unique_genes_msa.add_sequence( seqs.get_sequence( 0 ) )
 136     end
 137   end
 138   w = FastaWriter.new
 139   w.write(unique_genes_msa, "seqs_from_unique_genes.fasta")
 140   w.write(longest_non_unique_genes_msa, "longest_seqs_from_nonunique_genes.fasta")
 141 end