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