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   msa = nil
  17
  18   begin
  19     msa = f.create_msa_from_file( input, FastaParser.new() )
  20   rescue Exception => e
  21     puts "error: " + e.to_s
  22     exit
  23   end
  24
  25   gn_to_seqs = Hash.new
  26   unique_genes_msa = Msa.new
  27   longest_non_unique_genes_msa = Msa.new
  28   gn_re = /GN=(\S+)/
  29   fragment_re = /fragment/i
  30
  31   frag_counter = 0
  32   no_gn_counter = 0
  33
  34   for i in 0 ... msa.get_number_of_seqs()
  35     seq = msa.get_sequence( i )
  36     name = seq.get_name
  37     if fragment_re.match( name )
  38       puts "ignored because fragment: " + name
  39       frag_counter += 1
  40       next
  41     end
  42     gn_match = gn_re.match( name )
  43     unless gn_match
  44       puts "ignored because no GN=: " + name
  45       no_gn_counter += 1
  46       next
  47     end
  48     gn = gn_match[1]
  49     unless gn_to_seqs.has_key?(gn)
  50       gn_to_seqs[gn] = Msa.new
  51     end
  52     gn_to_seqs[gn].add_sequence(seq)
  53   end
  54
  55   puts "Sequeunces ignored because \"fragment\" in desc: " + frag_counter.to_s
  56   puts "Sequeunces ignored because no \"GN=\" in desc  : " + no_gn_counter.to_s
  57   puts
  58   puts
  59
  60   counter = 1
  61   gn_to_seqs.each_pair do |gene,seqs|
  62     if seqs.get_number_of_seqs > 1
  63       puts counter.to_s + ": " + gene
  64       puts seqs.to_fasta
  65       puts
  66       puts
  67       counter += 1
  68       longest = 0
  69       longest_seq = nil
  70       for j in 0 ... seqs.get_number_of_seqs()
  71         current = seqs.get_sequence( j )
  72         if current.get_length > longest
  73           longest =  current.get_length
  74           longest_seq = current
  75         end
  76       end
  77       longest_non_unique_genes_msa.add_sequence(longest_seq)
  78     else
  79       unique_genes_msa.add_sequence( seqs.get_sequence( 0 ) )
  80     end
  81   end
  82   w = FastaWriter.new
  83   w.write(unique_genes_msa, "seqs_from_unique_genes.fasta")
  84   w.write(longest_non_unique_genes_msa, "longest_seqs_from_nonunique_genes.fasta")
  85 end