2 # = lib/evo/apps/fasta_taxonomy_processor - FastaTaxonomyProcessor class
4 # Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
5 # License:: GNU Lesser General Public License (LGPL)
7 # $Id: fasta_taxonomy_processor.rb,v 1.4 2010/12/13 19:00:11 cmzmasek Exp $
10 require 'lib/evo/util/util'
11 require 'lib/evo/msa/msa_factory'
12 require 'lib/evo/msa/msa'
13 require 'lib/evo/io/msa_io'
14 require 'lib/evo/io/parser/sp_taxonomy_parser'
15 require 'lib/evo/io/parser/fasta_parser'
16 require 'lib/evo/io/writer/fasta_writer'
17 require 'lib/evo/io/writer/phylip_sequential_writer'
18 require 'lib/evo/util/command_line_arguments'
19 require 'lib/evo/apps/tseq_taxonomy_processor'
23 class FastaTaxonomyProcessor
25 PRG_NAME = "fasta_tap"
26 PRG_DATE = "2009.01.20"
27 PRG_DESC = "preprocessing of multiple sequence files in ncbi fasta format"
29 COPYRIGHT = "2009 Christian M Zmasek"
30 CONTACT = "phylosoft@gmail.com"
31 WWW = "www.phylosoft.org"
34 @tax_ids_to_sp_taxonomies = Hash.new()
39 Util.print_program_information( PRG_NAME,
48 if ARGV == nil || ARGV.length != 4
49 puts( "Usage: #{PRG_NAME}.rb <sp taxonomy file> <sequences in ncbi fasta format> <name for fasta outfile> <name for map outfile>" )
55 cla = CommandLineArguments.new( ARGV )
56 rescue ArgumentError => e
57 Util.fatal_error( PRG_NAME, "error: " + e.to_s )
59 allowed_opts = Array.new
60 disallowed = cla.validate_allowed_options_as_str( allowed_opts )
61 if ( disallowed.length > 0 )
62 Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
65 sp_taxonomy_infile = cla.get_file_name( 0 )
66 sequences_infile = cla.get_file_name( 1 )
67 sequences_outfile = cla.get_file_name( 2 )
68 mapping_outfile = cla.get_file_name( 3 )
70 Util.fatal_error_if_not_readable( PRG_NAME, sp_taxonomy_infile )
71 Util.fatal_error_if_not_readable( PRG_NAME, sequences_infile )
72 Util.fatal_error_if_not_writable( PRG_NAME, mapping_outfile )
73 Util.fatal_error_if_not_writable( PRG_NAME, sequences_outfile )
75 sp_taxonomies = SpTaxonomyParser.parse( sp_taxonomy_infile )
77 Util.print_message( PRG_NAME, "read in taxonomic data for " + sp_taxonomies.size.to_s + " species from: " + sp_taxonomy_infile )
79 fasta_parser = FastaParser.new
80 msa_fac = MsaFactory.new
82 seqs = msa_fac.create_msa_from_file( sequences_infile, fasta_parser )
84 Util.print_message( PRG_NAME, "read in " + seqs.get_number_of_seqs.to_s + " sequences from: " + sequences_infile )
86 removed = seqs.remove_redundant_sequences!( true, true )
89 Util.print_message( PRG_NAME, "going to ignore the following " + removed.size.to_s + " redundant sequences:" )
90 removed.each { | seq_name |
93 Util.print_message( PRG_NAME, "will process " + seqs.get_number_of_seqs.to_s + " non-redundant sequences" )
96 mapping_out = File.open( mapping_outfile, "a" )
98 for i in 0 ... seqs.get_number_of_seqs
99 seq = seqs.get_sequence( i )
100 seq.set_name( Util::normalize_seq_name( modify_name( seq, i, sp_taxonomies, mapping_out ), 10 ) )
105 w = FastaWriter.new()
107 w.set_max_name_length( 10 )
110 io.write_to_file( seqs, sequences_outfile, w )
111 rescue Exception => e
112 Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
116 Util.print_message( PRG_NAME, "wrote: " + mapping_outfile )
117 Util.print_message( PRG_NAME, "wrote: " + sequences_outfile )
118 Util.print_message( PRG_NAME, "OK" )
124 def modify_name( seq, i, sp_taxonomies, mapping_outfile )
128 seq_desc = seq.get_name
132 if seq_desc =~ /\[(.+)\]/
135 Util.fatal_error( PRG_NAME, "no taxonomy in [" + seq_desc + "]" )
138 matching_sp_taxonomy = nil
140 sp_taxonomies.each { |sp_taxonomy|
141 if ( sp_taxonomy.scientific_name == taxonomy_sn )
142 matching_sp_taxonomy = sp_taxonomy
146 if matching_sp_taxonomy == nil
147 Util.fatal_error( PRG_NAME, "taxonomy [" + taxonomy_sn + "] for [" + seq_desc + "] not found" )
150 new_name = i.to_s( 16 ) + "_" + matching_sp_taxonomy.code
153 if seq_desc =~ /gi\|(.+?)\|/
156 Util.fatal_error( PRG_NAME, "no gi in [" + seq_desc + "]" )
161 if seq_desc =~ /\|\s*([^|]+?)\s*\[/
165 if seq_name =~ /\[.+\]$/
166 # Redundant taxonomy information hides here.
167 seq_name = seq_name.sub(/\[.+\]$/, '')
169 if seq_name =~ /^\s*hypothetical\s+protein\s*/i
170 # Pointless information.
171 seq_name = seq_name.sub( /^\s*hypothetical\s+protein\s*/i, '' )
173 if seq_name =~ /^\s*conserved\s+hypothetical\s+protein\s*/i
174 # Pointless information.
175 seq_name = seq_name.sub( /^\s*conserved\s+hypothetical\s+protein\s*/i, '' )
179 mapping_outfile.print( new_name + "\t" +
180 TseqTaxonomyProcessor::TAXONOMY_CODE + matching_sp_taxonomy.code + "\t" +
181 TseqTaxonomyProcessor::TAXONOMY_ID + matching_sp_taxonomy.id + "\t" +
182 TseqTaxonomyProcessor::TAXONOMY_ID_TYPE + "ncbi" + "\t" +
183 TseqTaxonomyProcessor::TAXONOMY_SN + matching_sp_taxonomy.scientific_name + "\t" +
184 TseqTaxonomyProcessor::SEQ_ACCESSION + gi.to_s + "\t" +
185 TseqTaxonomyProcessor::SEQ_ACCESSION_SOURCE + "gi" + "\t" +
186 TseqTaxonomyProcessor::SEQ_NAME + seq_name + "\t" +
187 TseqTaxonomyProcessor::SEQ_MOL_SEQ + seq.get_sequence_as_string +
188 Constants::LINE_DELIMITER )
190 mapping_outfile.print( new_name + "\t" +
191 TseqTaxonomyProcessor::TAXONOMY_CODE + matching_sp_taxonomy.code + "\t" +
192 TseqTaxonomyProcessor::TAXONOMY_ID + matching_sp_taxonomy.id + "\t" +
193 TseqTaxonomyProcessor::TAXONOMY_ID_TYPE + "ncbi" + "\t" +
194 TseqTaxonomyProcessor::TAXONOMY_SN + matching_sp_taxonomy.scientific_name + "\t" +
195 TseqTaxonomyProcessor::SEQ_NAME + seq_name + "\t" +
196 TseqTaxonomyProcessor::SEQ_MOL_SEQ + seq.get_sequence_as_string +
197 Constants::LINE_DELIMITER )