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