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