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 = Set.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
  35   frag_counter = 0
  36   no_gn_counter = 0
  37
  38   for i in 0 ... msa.get_number_of_seqs()
  39     seq = msa.get_sequence( i )
  40     name = seq.get_name
  41     if all_names.include?( name )
  42       puts "error: " + name + " is not unique (#" + i.to_s + ")"
  43       exit
  44     else
  45       all_names << name
  46     end
  47     mol_seq = seq.get_sequence_as_string.upcase
  48     if all_seqs.include?( mol_seq )
  49       puts "error: sequence of " + name + " is not unique (#" +  i.to_s + ")"
  50       exit
  51     else
  52       all_seqs << mol_seq
  53     end
  54
  55     if fragment_re.match( name )
  56       puts "ignored because fragment: " + name
  57       frag_counter += 1
  58       next
  59     end
  60
  61     gn_match = gn_re.match( name )
  62     if IGNORE_SEQS_LACKING_GN
  63       unless gn_match
  64         puts "ignored because no GN=: " + name
  65         no_gn_counter += 1
  66         next
  67       end
  68     else
  69       unless gn_match
  70         puts "no GN=: " + name
  71       end
  72     end
  73
  74     gn =nil
  75     if gn_match
  76       gn = gn_match[1]
  77     else
  78       if IGNORE_SEQS_LACKING_GN
  79         puts "cannot be"
  80         exit
  81       end
  82       gn = name
  83     end
  84
  85     unless gn_to_seqs.has_key?(gn)
  86       gn_to_seqs[gn] = Msa.new
  87     end
  88     gn_to_seqs[gn].add_sequence(seq)
  89   end
  90
  91   puts "Sequences ignored because \"fragment\" in desc: " + frag_counter.to_s
  92   puts "Sequences ignored because no \"GN=\" in desc  : " + no_gn_counter.to_s
  93   puts
  94   puts
  95
  96   counter = 1
  97   gn_to_seqs.each_pair do |gene,seqs|
  98     if seqs.get_number_of_seqs > 1
  99       puts counter.to_s + ": " + gene
 100       puts seqs.to_fasta
 101       puts
 102       puts
 103       counter += 1
 104       longest = 0
 105       longest_seq = nil
 106       for j in 0 ... seqs.get_number_of_seqs()
 107         current = seqs.get_sequence( j )
 108         if current.get_length > longest
 109           longest =  current.get_length
 110           longest_seq = current
 111         end
 112       end
 113       longest_non_unique_genes_msa.add_sequence(longest_seq)
 114     else
 115       unique_genes_msa.add_sequence( seqs.get_sequence( 0 ) )
 116     end
 117   end
 118   w = FastaWriter.new
 119   w.write(unique_genes_msa, "seqs_from_unique_genes.fasta")
 120   w.write(longest_non_unique_genes_msa, "longest_seqs_from_nonunique_genes.fasta")
 121 end