+++ /dev/null
-#
-# = lib/evo/apps/domain_sequence_extractor.rb - DomainSequenceExtractor class
-#
-# Copyright:: Copyright (C) 2012 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id:Exp $
-
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-require 'lib/evo/io/parser/hmmscan_domain_extractor'
-
-module Evoruby
-
- class DomainSequenceExtractor
-
- PRG_NAME = "dsx"
- PRG_VERSION = "2.000"
- PRG_DESC = "extraction of domain sequences from hmmscan output"
- PRG_DATE = "20121001"
- COPYRIGHT = "2012 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- E_VALUE_THRESHOLD_OPTION = 'e'
- LENGTH_THRESHOLD_OPTION = 'l'
- ADD_POSITION_OPTION = 'p'
- ADD_DOMAIN_NUMBER_OPTION = 'd'
- ADD_SPECIES = 's'
- MIN_LINKER_OPT = 'ml'
- LOG_FILE_SUFFIX = '_domain_seq_extr.log'
- PASSED_SEQS_SUFFIX = '_with_passing_domains.fasta'
- FAILED_SEQS_SUFFIX = '_with_no_passing_domains.fasta'
- HELP_OPTION_1 = 'help'
- HELP_OPTION_2 = 'h'
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC ,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- ld = Constants::LINE_DELIMITER
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError
- Util.fatal_error( PRG_NAME, "error: " + $!, STDOUT )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if ( cla.get_number_of_files != 4 )
- print_help
- exit( -1 )
- end
-
- allowed_opts = Array.new
- allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
- allowed_opts.push( ADD_POSITION_OPTION )
- allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION )
- allowed_opts.push( LENGTH_THRESHOLD_OPTION )
- allowed_opts.push( ADD_SPECIES )
- allowed_opts.push( MIN_LINKER_OPT )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed,
- STDOUT )
- end
-
- domain_id = cla.get_file_name( 0 )
- hmmsearch_output = cla.get_file_name( 1 )
- fasta_sequence_file = cla.get_file_name( 2 )
- outfile = cla.get_file_name( 3 )
-
- if outfile.downcase.end_with?( ".fasta" )
- outfile = outfile[ 0 .. outfile.length - 7 ]
- elsif outfile.downcase.end_with?( ".fsa" )
- outfile = outfile[ 0 .. outfile.length - 5 ]
- end
-
-
- add_position = false
- if ( cla.is_option_set?( ADD_POSITION_OPTION ) )
- add_position = true
- end
-
- add_domain_number = false
- if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION ) )
- add_domain_number = true
- end
-
- add_species = false
- if cla.is_option_set? ADD_SPECIES
- add_species = true
- end
-
- e_value_threshold = -1.0
- if ( cla.is_option_set?( E_VALUE_THRESHOLD_OPTION ) )
- begin
- e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
- rescue ArgumentError => e
- Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- if ( e_value_threshold < 0.0 )
- Forester::Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
- end
- end
-
- length_threshold = -1
- if ( cla.is_option_set?( LENGTH_THRESHOLD_OPTION ) )
- begin
- length_threshold = cla.get_option_value_as_int( LENGTH_THRESHOLD_OPTION )
- rescue ArgumentError => e
- Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- if ( length_threshold < 0)
- Forester::Util.fatal_error( PRG_NAME, "attempt to use a negative length threshold", STDOUT )
- end
- end
-
-
- min_linker = nil
- if ( cla.is_option_set?( MIN_LINKER_OPT ) )
- begin
- min_linker = cla.get_option_value_as_int( MIN_LINKER_OPT )
- rescue ArgumentError => e
- Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- if ( !min_linker || min_linker > 100 || min_linker < -100 )
- Forester::Util.fatal_error( PRG_NAME, "unexpected value for min linker " + min_linker.to_s, STDOUT )
- end
- end
-
-
- log = String.new
-
- puts()
- puts( "Domain : " + domain_id )
- log << "Domain : " + domain_id + ld
- puts( "Hmmscan outputfile : " + hmmsearch_output )
- log << "Hmmscan outputfile : " + hmmsearch_output + ld
- puts( "Fasta sequencefile (complete sequences): " + fasta_sequence_file )
- log << "Fasta sequencefile (complete sequences): " + fasta_sequence_file + ld
- puts( "Outputfile : " + outfile + ".fasta" )
- log << "Outputfile : " + outfile + ld
- puts( "Passed sequences outfile (fasta) : " + outfile + PASSED_SEQS_SUFFIX )
- log << "Passed sequences outfile (fasta) : " + outfile + PASSED_SEQS_SUFFIX + ld
- puts( "Failed sequences outfile (fasta) : " + outfile + FAILED_SEQS_SUFFIX )
- log << "Failed sequences outfile (fasta) : " + outfile + FAILED_SEQS_SUFFIX + ld
- puts( "Logfile : " + outfile + LOG_FILE_SUFFIX )
- log << "Logfile : " + outfile + LOG_FILE_SUFFIX + ld
- if ( e_value_threshold >= 0.0 )
- puts( "E-value threshold : " + e_value_threshold.to_s )
- log << "E-value threshold : " + e_value_threshold.to_s + ld
- else
- puts( "E-value threshold : no threshold" )
- log << "E-value threshold : no threshold" + ld
- end
- if ( length_threshold > 0 )
- puts( "Length threshold : " + length_threshold.to_s )
- log << "Length threshold : " + length_threshold.to_s + ld
- else
- puts( "Length threshold : no threshold" )
- log << "Length threshold : no threshold" + ld
- end
- if ( min_linker )
- puts( "Min linker : " + min_linker.to_s )
- log << "Min linker : " + min_linker.to_s + ld
-
- end
-
-
- if ( add_position )
- puts( "Add positions (rel to complete seq) to extracted domains: true" )
- log << "Add positions (rel to complete seq) to extracted domains: true" + ld
- else
- puts( "Add positions (rel to complete seq) to extracted domains: false" )
- log << "Add positions (rel to complete seq) to extracted domains: false" + ld
- end
-
- if ( add_domain_number )
- puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): true" )
- log << "Add numbers to extracted domains (in case of more than one domain per complete seq): true" + ld
- else
- puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): false" )
- log << "Add numbers to extracted domains (in case of more than one domain per complete seq): false" + ld
- end
-
- puts
-
- domain_count = 0
- begin
- parser = HmmscanDomainExtractor.new()
- domain_count = parser.parse( domain_id,
- hmmsearch_output,
- fasta_sequence_file,
- outfile,
- outfile + PASSED_SEQS_SUFFIX,
- outfile + FAILED_SEQS_SUFFIX,
- e_value_threshold,
- length_threshold,
- add_position,
- add_domain_number,
- add_species,
- min_linker,
- log )
- rescue ArgumentError, IOError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
-
- rescue Exception => e
- puts e.backtrace
- Util.fatal_error( PRG_NAME, "unexpected exception: " + e.to_s, STDOUT )
-
- end
-
- puts
- Util.print_message( PRG_NAME, "extracted a total of " + domain_count.to_s + " domains" )
- Util.print_message( PRG_NAME, "wrote; " + outfile + ".fasta")
- Util.print_message( PRG_NAME, "wrote: " + outfile + LOG_FILE_SUFFIX )
- Util.print_message( PRG_NAME, "wrote: " + outfile + PASSED_SEQS_SUFFIX )
- Util.print_message( PRG_NAME, "wrote: " + outfile + FAILED_SEQS_SUFFIX )
-
- begin
- f = File.open( outfile + LOG_FILE_SUFFIX, 'a' )
- f.print( log )
- f.close
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- puts
- Util.print_message( PRG_NAME, "OK" )
- puts
-
- end
-
- def print_help()
- puts()
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb [options] <domain> <hmmscan outputfile> <file containing complete sequences in fasta format> <outputfile>" )
- puts()
- puts( " options: -" + E_VALUE_THRESHOLD_OPTION + "=<f>: E-value threshold, default is no threshold" )
- puts( " -" + LENGTH_THRESHOLD_OPTION + "=<i>: length threshold, default is no threshold" )
- puts( " -" + ADD_POSITION_OPTION + ": to add positions (rel to complete seq) to extracted domains" )
- puts( " -" + ADD_DOMAIN_NUMBER_OPTION + ": to add numbers to extracted domains (in case of more than one domain per complete seq) (example \"domain~2-3\")" )
- puts( " -" + ADD_SPECIES + ": to add species [in brackets]" )
- puts( " -" + MIN_LINKER_OPT + "=<i>: to extract pairs of same domains with a distance inbetween shorter than a given value" )
- puts()
- end
-
- end # class DomainSequenceExtractor
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/domains_to_forester - DomainsToForester class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: Exp $
-#
-# last modified: 06/11/2007
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/io/parser/fasta_parser'
-require 'lib/evo/sequence/protein_domain'
-require 'lib/evo/sequence/domain_structure'
-
-module Evoruby
-
- class DomainsToForester
-
- PRG_NAME = "d2f"
- PRG_DESC = "parsed hmmpfam output to forester format"
- PRG_VERSION = "1.001"
- PRG_DATE = "20120807"
- COPYRIGHT = "2012 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- E_VALUE_THRESHOLD_OPTION = "e"
- OVERWRITE_IF_SAME_FROM_TO_OPTION = "o"
- HELP_OPTION_1 = "help"
- HELP_OPTION_2 = "h"
-
- def parse( domains_list_file,
- original_seqs_file,
- outfile,
- column_delimiter,
- e_value_threshold,
- overwrite_if_same_from_to )
- Util.check_file_for_readability( domains_list_file )
- Util.check_file_for_readability( original_seqs_file )
- Util.check_file_for_writability( outfile )
-
- domain_structures = Hash.new() # protein name is key, domain structure is value
-
- f = MsaFactory.new
-
- original_seqs = f.create_msa_from_file( original_seqs_file, FastaParser.new )
- if ( original_seqs.get_number_of_seqs < 1 )
- error_msg = "\"" + original_seqs_file + "\" appears devoid of sequences in fasta-format"
- raise ArgumentError, error_msg
- end
-
- File.open( domains_list_file ) do | file |
- while line = file.gets
- if !is_ignorable?( line )
-
- a = line.split( column_delimiter )
- l = a.length
- if ( ( l < 4 ) || ( e_value_threshold >= 0.0 && l < 5 ) )
- error_msg = "unexpected format at line: " + line
- raise IOError, error_msg
- end
- protein_name = a[ 0 ]
- domain_name = a[ 1 ]
- seq_from = -1
- seq_to = -1
- ##########################################
- if domain_name =~ /RRM_\d/
- puts "ignoring " + line
- next
- end
- ##########################################
-
-
- begin
- seq_from = a[ 2 ].to_i
- rescue Exception
- error_msg = "failed to parse seq from from \"" + a[ 2 ] + "\" [line: " + line + "]"
- raise IOError, error_msg
- end
- begin
- seq_to = a[ 3 ].to_i
- rescue Exception
- error_msg = "failed to parse seq to from \"" + a[ 3 ] + "\" [line: " + line + "]"
- raise IOError, error_msg
- end
-
- e_value = -1
- if l > 4
- begin
- e_value = a[ 4 ].to_f
- rescue Exception
- error_msg = "failed to parse E-value from \"" + a[ 4 ] + "\" [line: " + line + "]"
- raise IOError, error_msg
- end
- end
-
- seq = original_seqs.get_by_name_start( protein_name )
-
- total_length = seq.get_length
-
- if ( ( ( e_value_threshold < 0.0 ) || ( e_value <= e_value_threshold ) ) )
- pd = ProteinDomain.new( domain_name, seq_from, seq_to, "", e_value )
- ds = nil
- if ( domain_structures.has_key?( protein_name ) )
- ds = domain_structures[ protein_name ]
- else
- ds = DomainStructure.new( total_length )
- domain_structures[ protein_name ] = ds
- end
- ds.add_domain( pd, overwrite_if_same_from_to )
- end
-
- end
- end
- end
-
- out = File.open( outfile, "a" )
- ds = domain_structures.sort
- for d in ds
- protein_name = d[ 0 ]
- domain_structure = d[ 1 ]
- out.print( protein_name.to_s )
- out.print( ": " )
- out.print( domain_structure.to_NHX )
- out.print( Constants::LINE_DELIMITER )
- end
-
- out.flush()
- out.close()
-
- end # parse
-
-
-
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if cla.get_number_of_files != 3
- print_help
- exit( -1 )
- end
-
- allowed_opts = Array.new
- allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
- allowed_opts.push( OVERWRITE_IF_SAME_FROM_TO_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed,
- STDOUT )
- end
-
- domains_list_file = cla.get_file_name( 0 )
- original_sequences_file = cla.get_file_name( 1 )
- outfile = cla.get_file_name( 2 )
-
-
- e_value_threshold = -1.0
- if cla.is_option_set?( E_VALUE_THRESHOLD_OPTION )
- begin
- e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- if ( e_value_threshold < 0.0 )
- Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
- end
- end
- overwrite_if_same_from_to = false
- if ( cla.is_option_set?( OVERWRITE_IF_SAME_FROM_TO_OPTION ) )
- overwrite_if_same_from_to = true
- end
-
- puts
- puts( "Domains list file : " + domains_list_file )
- puts( "Fasta sequencefile (complete sequences): " + original_sequences_file )
- puts( "Outputfile : " + outfile )
- if ( e_value_threshold >= 0.0 )
- puts( "E-value threshold : " + e_value_threshold.to_s )
- else
- puts( "E-value threshold : no threshold" )
- end
- if ( overwrite_if_same_from_to )
- puts( "Overwrite if same from and to : true" )
- else
- puts( "Overwrite if same from and to : false" )
- end
-
- puts
-
- begin
- parse( domains_list_file,
- original_sequences_file,
- outfile,
- " ",
- e_value_threshold,
- overwrite_if_same_from_to )
-
- rescue ArgumentError, IOError, StandardError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "unexpected exception: " + e.to_s, STDOUT )
- end
-
-
- puts
- Util.print_message( PRG_NAME, 'OK' )
- puts
-
- end
-
- private
-
- def print_help()
- puts
- puts( "Usage:" )
- puts
- puts( " " + PRG_NAME + ".rb [options] <domains list file (parsed hmmpfam output)> <file containing complete sequences in fasta format> <outputfile>" )
- puts()
- puts( " options: -" + E_VALUE_THRESHOLD_OPTION + "=<f> : E-value threshold, default is no threshold" )
- puts( " -" + OVERWRITE_IF_SAME_FROM_TO_OPTION + " : overwrite domain with same start and end with domain with better E-value" )
- puts
- end
-
-
-
- def is_ignorable?( line )
- return ( line !~ /[A-Za-z0-9-]/ || line =~ /^\s*#/)
- end
-
-
- end # class DomainsToForester
-
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/evo_nursery.rb - EvoNursery class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: evo_nursery.rb,v 1.11 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-
-require 'lib/evo/soft/fastme'
-require 'lib/evo/soft/tree_puzzle'
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/writer/phylip_sequential_writer'
-require 'lib/evo/io/parser/general_msa_parser'
-require 'lib/evo/io/writer/msa_writer'
-
-require 'iconv'
-
-module Evoruby
-
- class EvoNursery
- GAP_RATIO = 0.50
- GAP_RATIO_FOR_SEQS = 0.75
- MIN_LENGTH = 30
- MIN_SEQS = 4
- MAX_SEQS = 3000
- MAX_ALN_FILE_SIZE = 5000000
- MODEL = :auto
- RATES = :uniform
- FASTME_INITIAL_TREE = :GME
- ALN_NAME = '_align_'
- TREE_PUZZLE_OUTDIST = TreePuzzle::OUTDIST
- TREE_PUZZLE_OUTFILE = TreePuzzle::OUTFILE
- FASTME_OUTTREE = FastMe::OUTTREE
- FASTME_OUTPUT_D = FastMe::OUTPUT_D
-
- PRG_NAME = "evo_nursery"
- PRG_DATE = "2012.03.21"
- PRG_DESC = "pfam alignments to evolutionary trees"
- PRG_VERSION = "0.20"
- COPYRIGHT = "2009-2012 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- HELP_OPTION_1 = "help"
- HELP_OPTION_2 = "h"
-
- def run
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if RUBY_VERSION !~ /1.9/
- puts( "Your ruby version is #{RUBY_VERSION}, expected 1.9.x " )
- exit( -1 )
- end
-
- forester_home = Util.get_env_variable_value( Constants::FORESTER_HOME_ENV_VARIABLE )
- java_home = Util.get_env_variable_value( Constants::JAVA_HOME_ENV_VARIABLE )
- decorator = java_home + '/bin/java -cp ' + forester_home + '/java/forester.jar org.forester.application.decorator'
-
- if ( ARGV == nil || ARGV.length != 1 )
- help
- exit( -1 )
- end
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- help
- exit( 0 )
- end
-
- output_dir = cla.get_file_name( 0 )
-
- if output_dir !~ /\/$/
- output_dir = output_dir + '/'
- end
-
- if !File.exists?( output_dir )
- Util.fatal_error( PRG_NAME, output_dir.to_s + " does not exist", STDOUT )
- end
- ic = Iconv.new( 'UTF-8//IGNORE', 'UTF-8' )
- files = Dir.entries( "." )
- skipped = Array.new
- counter = 1
- analyzed = 0;
- begin
- files.each { |pfam_aln_file|
- if ( !File.directory?( pfam_aln_file ) &&
- pfam_aln_file !~ /^\./ &&
- pfam_aln_file !~ /.+\.tre$/ )
-
- tree_out_file = output_dir + File.basename( pfam_aln_file ) + ".xml"
-
- if File.exists?( tree_out_file )
- puts
- puts
- puts "***** skipping " + File.basename( pfam_aln_file ) + ", already exists"
- puts
- skipped.push( File.basename( pfam_aln_file ) + " [already exists]" )
- next
- end
-
- puts
- puts counter.to_s + ": " + pfam_aln_file.to_str
- counter += 1
- if File.size( pfam_aln_file ) > MAX_ALN_FILE_SIZE
- puts "***** skipping, file size: " + File.size( pfam_aln_file ).to_s
- skipped.push( File.basename( pfam_aln_file ) + " [file size: " + File.size( pfam_aln_file ).to_s + "]" )
- next
- end
-
- f = MsaFactory.new()
- msa = f.create_msa_from_file( pfam_aln_file, GeneralMsaParser.new() )
-
- if msa.get_number_of_seqs < MIN_SEQS || msa.get_number_of_seqs > MAX_SEQS
- puts "***** skipping, seqs: " + msa.get_number_of_seqs.to_s
- skipped.push( File.basename( pfam_aln_file ) + " [seqs: " + msa.get_number_of_seqs.to_s + "]" )
- next
- end
-
- msa.remove_gap_columns_w_gap_ratio!( GAP_RATIO )
-
- length = msa.get_length
- if length < MIN_LENGTH
- puts "***** skipping, length: " + length.to_s
- skipped.push( File.basename( pfam_aln_file ) + " [length: " + length.to_s + "]" )
- next
- end
-
- msa.remove_sequences_by_gap_ratio!( GAP_RATIO_FOR_SEQS )
-
- if msa.get_number_of_seqs < MIN_SEQS
- puts "***** skipping, seqs: " + msa.get_number_of_seqs.to_s
- skipped.push( File.basename( pfam_aln_file ) + " [seqs: " + msa.get_number_of_seqs.to_s + "]" )
- next
- end
-
- map_file = output_dir + File.basename( pfam_aln_file ) + ".map"
- f = File.open( map_file, 'a' )
- for i in 0 ... msa.get_number_of_seqs
- name = msa.get_sequence( i ).get_name()
- name =~ /(.+)_(.+)\/.+/
- acc = $1
- tax_code = $2
-
- mapping_str = i.to_s
- mapping_str << "\t"
- mapping_str << 'TAXONOMY_CODE:'
- mapping_str << tax_code
- mapping_str << "\t"
- mapping_str << 'SEQ_SYMBOL:'
- mapping_str << ( acc + '_' + tax_code )
- mapping_str << "\t"
- if ( acc.length < 6 )
- acc = acc + '_' + tax_code
- end
- mapping_str << 'SEQ_ACCESSION:'
- mapping_str << acc
- mapping_str << "\t"
- mapping_str << 'SEQ_ACCESSION_SOURCE:UniProtKB'
- mapping_str << "\t"
- mapping_str << 'NODE_NAME:'
- mapping_str << name
- f.print( mapping_str )
- f.print( "\n" )
- name = msa.get_sequence( i ).set_name( i.to_s )
- end
- f.close
-
- io = MsaIO.new()
- w = MsaWriter
- w = PhylipSequentialWriter.new()
- w.clean( true )
- w.set_max_name_length( 10 )
- if File.exists?( output_dir + ALN_NAME )
- File.unlink( output_dir + ALN_NAME )
- end
- io.write_to_file( msa, output_dir + ALN_NAME, w )
-
- tp = TreePuzzle.new()
- tp.run( output_dir + ALN_NAME,
- MODEL,
- RATES,
- msa.get_number_of_seqs )
-
- File.rename( output_dir + ALN_NAME, output_dir + File.basename( pfam_aln_file ) + ".aln" )
-
- fastme = FastMe.new()
- fastme.run( TREE_PUZZLE_OUTDIST, 0, FASTME_INITIAL_TREE )
-
- pfam_acc = nil
- pfam_de = nil
- File.open( pfam_aln_file ) do |file|
- while line = file.gets
- line = ic.iconv( line )
- if line =~ /^#=AC\s+(.+)/
- pfam_acc = $1
- end
- if line =~ /^#=DE\s+(.+)/
- pfam_de = $1
- end
- if pfam_acc && pfam_de
- break
- end
- end
- end
- if !pfam_acc || !pfam_de
- Util.fatal_error( PRG_NAME, "problem with " + pfam_aln_file.to_s, STDOUT )
- end
-
- puzzle_model = nil
- File.open( TREE_PUZZLE_OUTFILE ) do |file|
- while line = file.gets
- line = ic.iconv( line )
- if line =~ /^Model\s+of\s+substitution:\s+(.+)/
- puzzle_model = $1
- break
- end
- end
- end
- if !puzzle_model
- Util.fatal_error( PRG_NAME, "problem with puzzle outfile: " + TREE_PUZZLE_OUTFILE.to_s, STDOUT )
- end
-
- desc = pfam_de
- desc << ' | '
- desc << 'ML pwd estimation by TREE-PUZZLE version '
- desc << TreePuzzle::VERSION
- desc << ', model: '
- desc << puzzle_model
- desc << ', rates: '
- desc << RATES.to_s
- desc << '; tree estimation by FastME version '
- desc << FastMe::VERSION
- desc << ', initial tree: '
- desc << FASTME_INITIAL_TREE.to_s
- desc << '; aln length: '
- desc << msa.get_length.to_s
-
- cmd = decorator + " -table -p -pn=\"" + pfam_aln_file +
- "\" -pi=pfam:" + pfam_acc +
- " -pd=\"" + desc + "\" " +
- FASTME_OUTTREE + ' ' +
- map_file + ' ' + tree_out_file
-
- IO.popen( cmd , 'r+' ) do | pipe |
- pipe.close_write
- end
- analyzed += 1
-
- File.unlink( map_file )
- File.unlink(TREE_PUZZLE_OUTDIST)
- File.unlink( TREE_PUZZLE_OUTFILE )
- File.unlink( FASTME_OUTPUT_D )
- end
- }
- rescue ArgumentError, IOError, StandardError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- puts()
- puts( 'Skipped:' )
- puts()
- for i in 0 ... skipped.size
- puts i.to_s + ": " + skipped[ i ]
- end
-
- puts()
- puts( 'Skipped : ' + skipped.size.to_s + ' alignments' )
- puts( 'Analyzed: ' + analyzed.to_s + ' alignments' )
-
- puts( 'Min gap ratio for col del : ' + GAP_RATIO.to_s )
- puts( 'Min gap ratio for seq del : ' + GAP_RATIO_FOR_SEQS.to_s )
- puts( 'Minimal aln length : ' + MIN_LENGTH.to_s )
- puts( 'Minimal number of sequences: ' + MIN_SEQS.to_s )
- puts( 'Maximal number of sequences: ' + MAX_SEQS.to_s )
- puts( 'Maximal aln file size : ' + MAX_ALN_FILE_SIZE.to_s )
- puts( 'Model : ' + MODEL.to_s )
- puts( 'FastME initial tree: ' + FASTME_INITIAL_TREE.to_s )
-
- puts()
- puts( '[' + PRG_NAME + '] > OK' )
- puts()
-
- end # run
-
- private
-
- def help
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb <output dir> " )
- puts()
- end
-
-
- end # class EvoNursery
-
-end # module Evoruby
\ No newline at end of file
+++ /dev/null
-#
-# = lib/evo/apps/fasta_extractor.rb - FastaExtractor class
-#
-# Copyright:: Copyright (C) 2006-2008 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: fasta_extractor.rb,v 1.2 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-require 'lib/evo/util/util'
-require 'lib/evo/util/constants'
-require 'lib/evo/util/command_line_arguments'
-
-
-module Evoruby
-
- class FastaExtractor
-
- PRG_NAME = "fae"
- PRG_VERSION = "1.0.0"
- PRG_DESC = "extraction of nucleotide sequences from a fasta file by names from wublast search"
- PRG_DATE = "2008.08.09"
- COPYRIGHT = "2008-2009 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
- HELP_OPTION_1 = 'help'
- HELP_OPTION_2 = 'h'
-
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC ,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- ld = Constants::LINE_DELIMITER
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if ( cla.get_number_of_files != 3 )
- print_help
- exit( -1 )
- end
-
- allowed_opts = Array.new
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed,
- STDOUT )
- end
-
- input_file = cla.get_file_name( 0 )
- names_file = cla.get_file_name( 1 )
- output_file = cla.get_file_name( 2 )
-
- if !File.exist?( input_file )
- Util.fatal_error( PRG_NAME, "error: input file [#{input_file}] does not exist" )
- end
- if !File.exist?( names_file )
- Util.fatal_error( PRG_NAME, "error: names file [#{names_file}] does not exist" )
- end
- if File.exist?( output_file )
- Util.fatal_error( PRG_NAME, "error: [#{output_file }] already exists" )
- end
-
- names = extract_names_with_frames( names_file )
-
- extract_sequences( names, input_file, output_file )
-
- puts
- Util.print_message( PRG_NAME, "OK" )
- puts
-
- end
-
-
- def extract_names_with_frames( names_file )
- names = Hash.new()
- File.open( names_file ) do | file |
- while line = file.gets
- if ( !Util.is_string_empty?( line ) && !(line =~ /\s*#/ ) )
- if ( line =~ /(\S+)\s+([+|-]\d)\s+\d+\s+(\S+)/ )
- name = $1
- frame = $2
- e = $3
- names[ name ] = "[" + frame + "] [" + e + "]"
- end
- end
- end
- end
- names
- end
-
- def extract_sequences( names, fasta_file, output_file )
- output = File.open( output_file, "a" )
- matching_state = false
- counter = 0
- File.open( fasta_file ) do | file |
- while line = file.gets
- if !Util.is_string_empty?( line )
- if ( line =~ /\s*>\s*(.+)/ )
- name = $1
- if names.has_key?( name )
- matching_state = true
- counter += 1
- puts counter.to_s + ". " +name + " " + names[ name ]
- output.print( ">" + name + " " + names[ name ] )
- output.print( Evoruby::Constants::LINE_DELIMITER )
- else
- matching_state = false
- end
- elsif matching_state
- output.print( line )
- end
- end
- end
- end
- output.close()
- end
-
- def print_help()
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb <input fasta file> <names file based on blast output> <output file>" )
- puts()
- end
-
- end # class FastaExtractor
-end
\ No newline at end of file
+++ /dev/null
-#
-# = lib/evo/apps/fasta_taxonomy_processor - FastaTaxonomyProcessor class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: fasta_taxonomy_processor.rb,v 1.4 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-require 'lib/evo/util/util'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/msa/msa'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/parser/sp_taxonomy_parser'
-require 'lib/evo/io/parser/fasta_parser'
-require 'lib/evo/io/writer/fasta_writer'
-require 'lib/evo/io/writer/phylip_sequential_writer'
-require 'lib/evo/util/command_line_arguments'
-require 'lib/evo/apps/tseq_taxonomy_processor'
-
-module Evoruby
-
- class FastaTaxonomyProcessor
-
- PRG_NAME = "fasta_tap"
- PRG_DATE = "2009.01.20"
- PRG_DESC = "preprocessing of multiple sequence files in ncbi fasta format"
- PRG_VERSION = "1.00"
- COPYRIGHT = "2009 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- def initialize()
- @tax_ids_to_sp_taxonomies = Hash.new()
- end
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if ARGV == nil || ARGV.length != 4
- puts( "Usage: #{PRG_NAME}.rb <sp taxonomy file> <sequences in ncbi fasta format> <name for fasta outfile> <name for map outfile>" )
- puts()
- exit( -1 )
- end
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- allowed_opts = Array.new
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
- end
-
- sp_taxonomy_infile = cla.get_file_name( 0 )
- sequences_infile = cla.get_file_name( 1 )
- sequences_outfile = cla.get_file_name( 2 )
- mapping_outfile = cla.get_file_name( 3 )
-
- Util.fatal_error_if_not_readable( PRG_NAME, sp_taxonomy_infile )
- Util.fatal_error_if_not_readable( PRG_NAME, sequences_infile )
- Util.fatal_error_if_not_writable( PRG_NAME, mapping_outfile )
- Util.fatal_error_if_not_writable( PRG_NAME, sequences_outfile )
-
- sp_taxonomies = SpTaxonomyParser.parse( sp_taxonomy_infile )
-
- Util.print_message( PRG_NAME, "read in taxonomic data for " + sp_taxonomies.size.to_s + " species from: " + sp_taxonomy_infile )
-
- fasta_parser = FastaParser.new
- msa_fac = MsaFactory.new
-
- seqs = msa_fac.create_msa_from_file( sequences_infile, fasta_parser )
-
- Util.print_message( PRG_NAME, "read in " + seqs.get_number_of_seqs.to_s + " sequences from: " + sequences_infile )
-
- removed = seqs.remove_redundant_sequences!( true, true )
-
- if removed.size > 0
- Util.print_message( PRG_NAME, "going to ignore the following " + removed.size.to_s + " redundant sequences:" )
- removed.each { | seq_name |
- puts seq_name
- }
- Util.print_message( PRG_NAME, "will process " + seqs.get_number_of_seqs.to_s + " non-redundant sequences" )
- end
-
- mapping_out = File.open( mapping_outfile, "a" )
-
- for i in 0 ... seqs.get_number_of_seqs
- seq = seqs.get_sequence( i )
- seq.set_name( Util::normalize_seq_name( modify_name( seq, i, sp_taxonomies, mapping_out ), 10 ) )
- end
-
- io = MsaIO.new()
-
- w = FastaWriter.new()
-
- w.set_max_name_length( 10 )
- w.clean( true )
- begin
- io.write_to_file( seqs, sequences_outfile, w )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
- end
- mapping_out.close()
-
- Util.print_message( PRG_NAME, "wrote: " + mapping_outfile )
- Util.print_message( PRG_NAME, "wrote: " + sequences_outfile )
- Util.print_message( PRG_NAME, "OK" )
-
- end
-
- private
-
- def modify_name( seq, i, sp_taxonomies, mapping_outfile )
-
- #i = i + 1792
-
- seq_desc = seq.get_name
-
- taxonomy_sn = nil
-
- if seq_desc =~ /\[(.+)\]/
- taxonomy_sn = $1
- else
- Util.fatal_error( PRG_NAME, "no taxonomy in [" + seq_desc + "]" )
- end
-
- matching_sp_taxonomy = nil
-
- sp_taxonomies.each { |sp_taxonomy|
- if ( sp_taxonomy.scientific_name == taxonomy_sn )
- matching_sp_taxonomy = sp_taxonomy
- end
- }
-
- if matching_sp_taxonomy == nil
- Util.fatal_error( PRG_NAME, "taxonomy [" + taxonomy_sn + "] for [" + seq_desc + "] not found" )
- end
-
- new_name = i.to_s( 16 ) + "_" + matching_sp_taxonomy.code
-
- gi = nil
- if seq_desc =~ /gi\|(.+?)\|/
- gi = $1
- else
- Util.fatal_error( PRG_NAME, "no gi in [" + seq_desc + "]" )
- end
-
- seq_name = ""
-
- if seq_desc =~ /\|\s*([^|]+?)\s*\[/
- seq_name = $1
- end
-
- if seq_name =~ /\[.+\]$/
- # Redundant taxonomy information hides here.
- seq_name = seq_name.sub(/\[.+\]$/, '')
- end
- if seq_name =~ /^\s*hypothetical\s+protein\s*/i
- # Pointless information.
- seq_name = seq_name.sub( /^\s*hypothetical\s+protein\s*/i, '' )
- end
- if seq_name =~ /^\s*conserved\s+hypothetical\s+protein\s*/i
- # Pointless information.
- seq_name = seq_name.sub( /^\s*conserved\s+hypothetical\s+protein\s*/i, '' )
- end
-
- if gi != nil
- mapping_outfile.print( new_name + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_CODE + matching_sp_taxonomy.code + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_ID + matching_sp_taxonomy.id + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_ID_TYPE + "ncbi" + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_SN + matching_sp_taxonomy.scientific_name + "\t" +
- TseqTaxonomyProcessor::SEQ_ACCESSION + gi.to_s + "\t" +
- TseqTaxonomyProcessor::SEQ_ACCESSION_SOURCE + "gi" + "\t" +
- TseqTaxonomyProcessor::SEQ_NAME + seq_name + "\t" +
- TseqTaxonomyProcessor::SEQ_MOL_SEQ + seq.get_sequence_as_string +
- Constants::LINE_DELIMITER )
- else
- mapping_outfile.print( new_name + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_CODE + matching_sp_taxonomy.code + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_ID + matching_sp_taxonomy.id + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_ID_TYPE + "ncbi" + "\t" +
- TseqTaxonomyProcessor::TAXONOMY_SN + matching_sp_taxonomy.scientific_name + "\t" +
- TseqTaxonomyProcessor::SEQ_NAME + seq_name + "\t" +
- TseqTaxonomyProcessor::SEQ_MOL_SEQ + seq.get_sequence_as_string +
- Constants::LINE_DELIMITER )
-
- end
- new_name
- end
-
- end
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/hmmscan_parser.rb - HmmscanParser class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: hmmscan_parser.rb,v 1.5 2010/12/13 19:00:11 cmzmasek Exp $
-#
-# last modified: 11/24/2009
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-
-module Evoruby
-
- class HmmscanParser
-
- PRG_NAME = "hsp"
- PRG_VERSION = "1.0.1"
- PRG_DESC = "hmmscan parser"
- PRG_DATE = "2009.11.24"
- COPYRIGHT = "2009 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- DELIMITER_OPTION = "d"
- E_VALUE_THRESHOLD_OPTION = "e"
- IGNORE_DUF_OPTION = "i"
- PARSE_OUT_DESCRIPITION_OPTION = "a"
- HELP_OPTION_1 = "help"
- HELP_OPTION_2 = "h"
-
- def initialize
- @domain_counts = Hash.new
- end
-
- # raises ArgumentError, IOError
- def parse( inpath,
- outpath,
- column_delimiter,
- e_value_threshold,
- ignore_dufs,
- get_descriptions )
- Util.check_file_for_readability( inpath )
- Util.check_file_for_writability( outpath )
-
- outfile = File.open( outpath, "a" )
-
- query = String.new
- desc = String.new
- model = String.new
- env_from = String.new
- env_to = String.new
- i_e_value = String.new
-
- queries_count = 0
-
- nl = Constants::LINE_DELIMITER
-
- File.open( inpath ) do | file |
- while line = file.gets
- if !HmmscanParser.is_ignorable?( line ) && line =~ /^\S+\s+\S/
-
- # tn acc tlen query acc qlen Evalue score bias # of c-E i-E score bias hf ht af at ef et acc desc
- # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
- line =~ /^(\S+)\s+(\S+)\s+(\d+)\s+(\S+)\s+(\S+)\s+(\d+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\d+)\s+(\d+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\S+)\s+(.*)/
-
- model = $1
- query = $4
- i_e_value = $13.to_f
- env_from = $20.to_i
- env_to = $21.to_i
-
- if ( ( ( e_value_threshold < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
- ( !ignore_dufs || ( model !~ /^DUF\d+/ ) ) )
- count_model( model )
- outfile.print( query +
- column_delimiter )
- if ( get_descriptions )
- outfile.print( desc +
- column_delimiter )
- end
- outfile.print( model +
- column_delimiter +
- env_from.to_s +
- column_delimiter +
- env_to.to_s +
- column_delimiter +
- i_e_value.to_s )
- outfile.print( nl )
- end
- end
- end # while line = file.gets
- end
- outfile.flush()
- outfile.close()
-
- return queries_count
-
- end # def parse
-
- def count_model( model )
- if ( @domain_counts.has_key?( model ) )
- count = @domain_counts[ model ].to_i
- count += 1
- @domain_counts[ model ] = count
- else
- @domain_counts[ model ] = 1
- end
- end
-
-
- def get_domain_counts()
- return @domain_counts
- end
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if ( cla.get_number_of_files != 2 )
- print_help
- exit( -1 )
- end
-
- allowed_opts = Array.new
- allowed_opts.push( DELIMITER_OPTION )
- allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
- allowed_opts.push( IGNORE_DUF_OPTION )
- allowed_opts.push( PARSE_OUT_DESCRIPITION_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed,
- STDOUT )
- end
-
- inpath = cla.get_file_name( 0 )
- outpath = cla.get_file_name( 1 )
-
- column_delimiter = "\t"
- if ( cla.is_option_set?( DELIMITER_OPTION ) )
- begin
- column_delimiter = cla.get_option_value( DELIMITER_OPTION )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
-
- e_value_threshold = -1.0
- if ( cla.is_option_set?( E_VALUE_THRESHOLD_OPTION ) )
- begin
- e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- if ( e_value_threshold < 0.0 )
- Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
- end
- end
-
- ignore_dufs = false
- if ( cla.is_option_set?( IGNORE_DUF_OPTION ) )
- ignore_dufs = true
- end
-
- parse_descriptions = false
- if ( cla.is_option_set?( PARSE_OUT_DESCRIPITION_OPTION ) )
- parse_descriptions = true
- end
-
- puts()
- puts( "hmmpfam outputfile: " + inpath )
- puts( "outputfile : " + outpath )
- if ( e_value_threshold >= 0.0 )
- puts( "E-value threshold : " + e_value_threshold.to_s )
- else
- puts( "E-value threshold : no threshold" )
- end
- if ( parse_descriptions )
- puts( "parse descriptions: true" )
- else
- puts( "parse descriptions: false" )
- end
- if ( ignore_dufs )
- puts( "ignore DUFs : true" )
- else
- puts( "ignore DUFs : false" )
- end
- if ( column_delimiter == "\t" )
- puts( "column delimiter : TAB" )
- else
- puts( "column delimiter : " + column_delimiter )
- end
- puts()
-
- begin
- queries_count = parse( inpath,
- outpath,
- column_delimiter,
- e_value_threshold,
- ignore_dufs,
- parse_descriptions )
- rescue ArgumentError, IOError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- domain_counts = get_domain_counts()
-
- puts
- puts( "read output for a total of " + queries_count.to_s + " query sequences" )
- puts
- puts( "domain counts (considering potential E-value threshold and ignoring of DUFs):" )
- puts( "(number of different domains: " + domain_counts.length.to_s + ")" )
- puts
- puts( Util.draw_histogram( domain_counts, "#" ) )
- puts
- Util.print_message( PRG_NAME, 'OK' )
- puts
-
- end # def run()
-
- def print_help()
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb [options] <hmmscan outputfile> <outputfile>" )
- puts()
- puts( " options: -" + DELIMITER_OPTION + ": column delimiter for outputfile, default is TAB" )
- puts( " -" + E_VALUE_THRESHOLD_OPTION + ": E-value threshold, default is no threshold" )
- puts( " -" + PARSE_OUT_DESCRIPITION_OPTION + ": parse query description (in addition to query name)" )
- puts( " -" + IGNORE_DUF_OPTION + ": ignore DUFs" )
- puts()
- end
-
-
- private
-
-
- def HmmscanParser.is_ignorable?( line )
- return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
- end
-
- end # class HmmscanParser
-
-end # module Evoruby
\ No newline at end of file
+++ /dev/null
-#
-# = lib/evo/apps/msa_processor.rb - MsaProcessor class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: msa_processor.rb,v 1.33 2010/12/13 19:00:10 cmzmasek Exp $
-#
-
-require 'date'
-require 'set'
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/writer/phylip_sequential_writer'
-require 'lib/evo/io/writer/nexus_writer'
-require 'lib/evo/io/writer/fasta_writer'
-require 'lib/evo/io/parser/fasta_parser'
-require 'lib/evo/io/parser/general_msa_parser'
-require 'lib/evo/io/writer/msa_writer'
-
-module Evoruby
-
- class MsaProcessor
-
- PRG_NAME = "msa_pro"
- PRG_DATE = "2012.05.11"
- PRG_DESC = "processing of multiple sequence alignments"
- PRG_VERSION = "1.06"
- COPYRIGHT = "2008-2010 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
-
- NAME_LENGTH_DEFAULT = 10
- WIDTH_DEFAULT_FASTA = 60
- INPUT_TYPE_OPTION = "i"
- OUTPUT_TYPE_OPTION = "o"
- MAXIMAL_NAME_LENGTH_OPTION = "n"
- WIDTH_OPTION = "w"
- CLEAN_UP_SEQ_OPTION = "c"
- REM_RED_OPTION = "rem_red"
- REMOVE_GAP_COLUMNS_OPTION = "rgc"
- REMOVE_GAP_ONLY_COLUMNS = "rgoc"
- REMOVE_COLUMNS_GAP_RATIO_OPTION = "rr"
- REMOVE_ALL_GAP_CHARACTERS_OPTION = "rg"
- REMOVE_ALL_SEQUENCES_LISTED_OPTION = "r"
- KEEP_ONLY_SEQUENCES_LISTED_OPTION = "k"
-
- KEEP_MATCHING_SEQUENCES_OPTION = "mk"
- REMOVE_MATCHING_SEQUENCES_OPTION = "mr"
-
- TRIM_OPTION = "t"
- REMOVE_SEQS_GAP_RATIO_OPTION = "rsgr"
- REMOVE_SEQS_NON_GAP_LENGTH_OPTION = "rsl"
- SPLIT = "split"
- LOG_SUFFIX = "_msa_pro.log"
- HELP_OPTION_1 = "help"
- HELP_OPTION_2 = "h"
-
-
- def initialize()
- @input_format_set = false
- @output_format_set = false
- @fasta_input = false
- @phylip_input = true
- @name_length = NAME_LENGTH_DEFAULT
- @name_length_set = false
- @width = WIDTH_DEFAULT_FASTA # fasta only
- @pi_output = true
- @fasta_output = false
- @nexus_output = false
- @clean = false # phylip only
- @rgc = false
- @rgoc = false
- @rg = false # fasta only
- @rem_red = false
- @rgr = -1
- @rsgr = -1
- @rsl = -1
- @remove_matching = nil
- @keep_matching = nil
-
- @seqs_name_file = nil
- @remove_seqs = false
- @keep_seqs = false
- @trim = false
- @split = -1
- @first = -1
- @last = -1
- end
-
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if ( ARGV == nil || ARGV.length < 1 )
- Util.print_message( PRG_NAME, "Illegal number of arguments" )
- print_help
- exit( -1 )
- end
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "Error: " + e.to_s, STDOUT )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if ( cla.get_number_of_files != 2 || ARGV.length < 2 )
- Util.print_message( PRG_NAME, "Illegal number of arguments" )
- print_help
- exit( -1 )
- end
-
- allowed_opts = Array.new
- allowed_opts.push( INPUT_TYPE_OPTION )
- allowed_opts.push( OUTPUT_TYPE_OPTION )
- allowed_opts.push( MAXIMAL_NAME_LENGTH_OPTION )
- allowed_opts.push( WIDTH_OPTION )
- allowed_opts.push( CLEAN_UP_SEQ_OPTION )
- allowed_opts.push( REMOVE_GAP_COLUMNS_OPTION )
- allowed_opts.push( REMOVE_GAP_ONLY_COLUMNS )
- allowed_opts.push( REMOVE_COLUMNS_GAP_RATIO_OPTION )
- allowed_opts.push( REMOVE_ALL_GAP_CHARACTERS_OPTION )
- allowed_opts.push( REMOVE_ALL_SEQUENCES_LISTED_OPTION )
- allowed_opts.push( KEEP_ONLY_SEQUENCES_LISTED_OPTION )
- allowed_opts.push( TRIM_OPTION )
- allowed_opts.push( REMOVE_SEQS_GAP_RATIO_OPTION )
- allowed_opts.push( REMOVE_SEQS_NON_GAP_LENGTH_OPTION )
- allowed_opts.push( SPLIT )
- allowed_opts.push( REM_RED_OPTION )
- allowed_opts.push( KEEP_MATCHING_SEQUENCES_OPTION )
- allowed_opts.push( REMOVE_MATCHING_SEQUENCES_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed )
- end
-
- input = cla.get_file_name( 0 )
- output = cla.get_file_name( 1 )
-
- analyze_command_line( cla )
-
- begin
- Util.check_file_for_readability( input )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- begin
- Util.check_file_for_writability( output )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if ( @rg )
- set_pi_output( false )
- set_fasta_output( true )
- set_nexus_output( false )
- end
-
- if ( !@input_format_set )
- fasta_like = false
- begin
- fasta_like = Util.looks_like_fasta?( input )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- @fasta_input = fasta_like
- @phylip_input = !fasta_like
- if ( !@output_format_set )
- @fasta_output = fasta_like
- @pi_output = !fasta_like
- @nexus_output = false
- end
- end
-
- ld = Constants::LINE_DELIMITER
- log = PRG_NAME + " " + PRG_VERSION + " [" + PRG_DATE + "]" + " LOG" + ld
- now = DateTime.now
- log << "Date/time: " + now.to_s + ld
-
- puts()
- puts( "Input alignment : " + input )
- log << "Input alignment : " + input + ld
- puts( "Output alignment : " + output )
- log << "Output alignment : " + output + ld
- if ( @phylip_input )
- puts( "Input is : Phylip, or something like it" )
- log << "Input is : Phylip, or something like it" + ld
- elsif ( @fasta_input )
- puts( "Input is : Fasta" )
- log << "Input is : Fasta" + ld
- end
- if( @rgr >= 0 )
- puts( "Max col gap ratio: " + @rgr.to_s )
- log << "Max col gap ratio: " + @rgr.to_s + ld
- elsif ( @rgc )
- puts( "Remove gap colums" )
- log << "Remove gap colums" + ld
- elsif( @rgoc )
- puts( "Remove gap only colums" )
- log << "Remove gap only colums" + ld
- end
- if ( @clean )
- puts( "Clean up : true" )
- log << "Clean up : true" + ld
- end
-
- if ( @pi_output )
- puts( "Output is : Phylip interleaved" )
- log << "Output is : Phylip interleaved" + ld
- elsif ( @fasta_output )
- puts( "Output is : Fasta" )
- log << "Output is : Fasta" + ld
- if ( @width )
- puts( "Width : " + @width.to_s )
- log << "Width : " + @width.to_s + ld
- end
- if ( @rg )
- puts( "Remove all gap characters (alignment is destroyed)" )
- log << "Remove all gap characters (alignment is destroyed)" + ld
- end
- elsif ( @nexus_output )
- puts( "Output is : Nexus" )
- log << "Output is : Nexus" + ld
- end
- if ( @name_length_set || !@fasta_output )
- puts( "Max name length : " + @name_length.to_s )
- log << "Max name length : " + @name_length.to_s + ld
- end
- if( @rsgr >= 0 )
- puts( "Remove sequences for which the gap ratio > " + @rsgr.to_s )
- log << "Remove sequences for which the gap ratio > " + @rsgr.to_s + ld
- end
- if( @rsl >= 0 )
- puts( "Remove sequences with less than " + @rsl.to_s + " non-gap characters" )
- log << "Remove sequences with less than " + @rsl.to_s + " non-gap characters" + ld
- end
- if ( @remove_seqs )
- puts( "Remove sequences listed in: " + @seqs_name_file )
- log << "Remove sequences listed in: " + @seqs_name_file + ld
- elsif ( @keep_seqs )
- puts( "Keep only sequences listed in: " + @seqs_name_file )
- log << "Keep only sequences listed in: " + @seqs_name_file + ld
- end
- if ( @trim )
- puts( "Keep only columns from: "+ @first.to_s + " to " + @last.to_s )
- log << "Keep only columns from: "+ @first.to_s + " to " + @last.to_s + ld
- end
- if ( @rem_red )
- puts( "Remove redundant sequences: true" )
- log << "Remove redundant sequences: true" + ld
- end
- if ( @split > 0 )
- puts( "Split : " + @split.to_s )
- log << "Split : " + @split.to_s + ld
- end
- puts()
-
- f = MsaFactory.new()
-
- msa = nil
-
- begin
- if ( @phylip_input )
- msa = f.create_msa_from_file( input, GeneralMsaParser.new() )
- elsif ( @fasta_input )
- msa = f.create_msa_from_file( input, FastaParser.new() )
- end
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- if ( msa.is_aligned() )
- Util.print_message( PRG_NAME, "Length of original alignment : " + msa.get_length.to_s )
- log << "Length of original alignment : " + msa.get_length.to_s + ld
- gp = msa.calculate_gap_proportion
- Util.print_message( PRG_NAME, "Gap-proportion of original alignment : " + gp.to_s )
- log << "Gap-proportion of original alignment : " + gp.to_s + ld
- else
- Util.print_message( PRG_NAME, "the input is not aligned" )
- log << "The input is not aligned" + ld
- end
-
- all_names = Set.new()
- for i in 0 ... msa.get_number_of_seqs()
- current_name = msa.get_sequence( i ).get_name
- if all_names.include?( current_name )
- Util.print_warning_message( PRG_NAME, "sequence name [" + current_name + "] is not unique" )
- else
- all_names.add( current_name )
- end
- end
-
- begin
-
- if ( @remove_seqs || @keep_seqs )
- names = Util.file2array( @seqs_name_file, true )
- if ( names == nil || names.length() < 1 )
- error_msg = "file \"" + @seqs_name_file.to_s + "\" appears empty"
- Util.fatal_error( PRG_NAME, error_msg )
- end
-
- if ( @remove_seqs )
- c = 0
- for i in 0 ... names.length()
- to_delete = msa.find_by_name( names[ i ], true, false )
- if ( to_delete.length() < 1 )
- error_msg = "sequence name \"" + names[ i ] + "\" not found"
- Util.fatal_error( PRG_NAME, error_msg )
- elsif ( to_delete.length() > 1 )
- error_msg = "sequence name \"" + names[ i ] + "\" is not unique"
- Util.fatal_error( PRG_NAME, error_msg )
- else
- msa.remove_sequence!( to_delete[ 0 ] )
- c += 1
- end
- end
- Util.print_message( PRG_NAME, "Removed " + c.to_s + " sequences" )
- log << "Removed " + c.to_s + " sequences" + ld
- elsif ( @keep_seqs )
- msa_new = Msa.new()
- r = 0
- k = 0
- for j in 0 ... msa.get_number_of_seqs()
- if ( names.include?( msa.get_sequence( j ).get_name() ) )
- msa_new.add_sequence( msa.get_sequence( j ) )
- k += 1
- else
- r += 1
- end
- end
- msa = msa_new
- Util.print_message( PRG_NAME, "Kept " + k.to_s + " sequences" )
- log << "Kept " + k.to_s + " sequences" + ld
- Util.print_message( PRG_NAME, "Removed " + r.to_s + " sequences" )
- log << "removed " + r.to_s + " sequences" + ld
- end
- end
-
- if ( @trim )
- msa.trim!( @first, @last )
- end
- if( @rgr >= 0 )
- msa.remove_gap_columns_w_gap_ratio!( @rgr )
- elsif ( @rgc )
- msa.remove_gap_columns!()
- elsif( @rgoc )
- msa.remove_gap_only_columns!()
- end
- if( @rsgr >= 0 )
- n = msa.get_number_of_seqs()
- removed = msa.remove_sequences_by_gap_ratio!( @rsgr )
- k = msa.get_number_of_seqs()
- r = n - k
- Util.print_message( PRG_NAME, "Kept " + k.to_s + " sequences" )
- log << "Kept " + k.to_s + " sequences" + ld
- Util.print_message( PRG_NAME, "Removed " + r.to_s + " sequences" )
- log << "Removed " + r.to_s + " sequences:" + ld
- removed.each { | seq_name |
- log << " " + seq_name + ld
- }
- end
- if( @rsl >= 0 )
- n = msa.get_number_of_seqs()
- removed = msa.remove_sequences_by_non_gap_length!( @rsl )
- k = msa.get_number_of_seqs()
- r = n - k
- Util.print_message( PRG_NAME, "Kept " + k.to_s + " sequences" )
- log << "Kept " + k.to_s + " sequences" + ld
- Util.print_message( PRG_NAME, "Removed " + r.to_s + " sequences" )
- log << "Removed " + r.to_s + " sequences:" + ld
- removed.each { | seq_name |
- log << " " + seq_name + ld
- }
- end
- if ( @keep_matching )
- n = msa.get_number_of_seqs
- to_be_removed = Set.new
- for ii in 0 ... n
- seq = msa.get_sequence( ii )
- if !seq.get_name.downcase.index( @keep_matching.downcase )
- to_be_removed.add( ii )
- end
- end
- to_be_removed_ary = to_be_removed.to_a.sort.reverse
- to_be_removed_ary.each { | index |
- msa.remove_sequence!( index )
- }
- # msa = sort( msa )
- end
- if ( @remove_matching )
- n = msa.get_number_of_seqs
- to_be_removed = Set.new
- for iii in 0 ... n
-
- seq = msa.get_sequence( iii )
-
- if seq.get_name.downcase.index( @remove_matching.downcase )
- to_be_removed.add( iii )
- end
- end
- to_be_removed_ary = to_be_removed.to_a.sort.reverse
- to_be_removed_ary.each { | index |
- msa.remove_sequence!( index )
- }
- msa = sort( msa )
- end
-
-
-
- if ( @split > 0 )
- begin
- msas = msa.split( @split, true )
- io = MsaIO.new()
- w = MsaWriter
- if ( @pi_output )
- w = PhylipSequentialWriter.new()
- w.clean( @clean )
- w.set_max_name_length( @name_length )
- elsif( @fasta_output )
- w = FastaWriter.new()
- w.set_line_width( @width )
- if ( @rg )
- w.remove_gap_chars( true )
- Util.print_warning_message( PRG_NAME, "removing gap character, the output is likely to become unaligned" )
- log << "removing gap character, the output is likely to become unaligned" + ld
- end
- w.clean( @clean )
- if ( @name_length_set )
- w.set_max_name_length( @name_length )
- end
- elsif( @nexus_output )
- w = NexusWriter.new()
- w.clean( @clean )
- w.set_max_name_length( @name_length )
- end
- i = 0
- for m in msas
- i = i + 1
- io.write_to_file( m, output + "_" + i.to_s, w )
- end
- Util.print_message( PRG_NAME, "wrote " + msas.length.to_s + " files" )
- log << "wrote " + msas.length.to_s + " files" + ld
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- end
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- if ( @split <= 0 )
-
- unless ( @rg )
- if ( msa.is_aligned() )
- Util.print_message( PRG_NAME, "Length of processed alignment : " + msa.get_length.to_s )
- log << "Length of processed alignment : " + msa.get_length.to_s + ld
- gp = msa.calculate_gap_proportion
- Util.print_message( PRG_NAME, "Gap-proportion of processed alignment: " + gp.to_s )
- log << "Gap-proportion of processed alignment: " + gp.to_s + ld
- else
- Util.print_warning_message( PRG_NAME, "output is not aligned" )
- log << "output is not aligned" + ld
- end
- end
-
- if @rem_red
- removed = msa.remove_redundant_sequences!( true, false )
- if removed.size > 0
- identicals = msa.get_identical_seqs_detected
- log << "the following " + identicals.size.to_s + " sequences are identical:" + ld
- identicals.each { | s |
- log << s + ld
- }
- log << "ignoring the following " + removed.size.to_s + " redundant sequences:" + ld
- removed.each { | seq_name |
- log << seq_name + ld
- }
- Util.print_message( PRG_NAME, "will store " + msa.get_number_of_seqs.to_s + " non-redundant sequences" )
- log << "will store " + msa.get_number_of_seqs.to_s + " non-redundant sequences" + ld
- end
- end
-
- io = MsaIO.new()
-
- w = MsaWriter
-
- if ( @pi_output )
- w = PhylipSequentialWriter.new()
- w.clean( @clean )
- w.set_max_name_length( @name_length )
- elsif( @fasta_output )
- w = FastaWriter.new()
- w.set_line_width( @width )
- if ( @rg )
- w.remove_gap_chars( true )
- Util.print_warning_message( PRG_NAME, "removing gap characters, the output is likely to become unaligned" )
- log << "removing gap character, the output is likely to become unaligned" + ld
- end
- w.clean( @clean )
- if ( @name_length_set )
- w.set_max_name_length( @name_length )
- end
- elsif( @nexus_output )
- w = NexusWriter.new()
- w.clean( @clean )
- w.set_max_name_length( @name_length )
- end
-
-
- begin
- io.write_to_file( msa, output, w )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- begin
- f = File.open( output + LOG_SUFFIX, 'a' )
- f.print( log )
- f.close
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
-
- end
- Util.print_message( PRG_NAME, "OK" )
- puts
- end
-
-
- private
-
- def sort( msa )
- names = Set.new
- for i in 0 ... msa.get_number_of_seqs
- name = msa.get_sequence( i ).get_name
- names.add( name )
- end
- sorted_ary = names.to_a.sort
- new_msa = Msa.new
- sorted_ary.each { | seq_name |
- seq = msa.get_sequence( msa.find_by_name( seq_name, true, false )[ 0 ] )
- new_msa.add_sequence( seq )
- }
- new_msa
- end
-
- def set_fasta_input( fi = true )
- @fasta_input = fi
- @input_format_set = true
- end
- def set_phylip_input( pi = true )
- @phylip_input = pi
- @input_format_set = true
- end
- def set_name_length( i )
- @name_length = i
- @name_length_set = true
- end
- def set_width( i )
- @width = i
- end
- def set_fasta_output( fo = true )
- @fasta_output = fo
- @output_format_set = true
- end
- def set_pi_output( pso = true )
- @pi_output = pso
- @output_format_set = true
- end
- def set_nexus_output( nexus = true )
- @nexus_output = nexus
- @output_format_set = true
- end
- def set_clean( c = true )
- @clean = c
- end
- def set_remove_gap_columns( rgc = true )
- @rgc = rgc
- end
- def set_remove_gap_only_columns( rgoc = true )
- @rgoc = rgoc
- end
- def set_remove_gaps( rg = true )
- @rg = rg
- end
- def set_remove_gap_ratio( rgr )
- @rgr = rgr
- end
- def set_remove_seqs_gap_ratio( rsgr )
- @rsgr = rsgr
- end
- def set_remove_seqs_min_non_gap_length( rsl )
- @rsl = rsl
- end
- def set_remove_seqs( file )
- @seqs_name_file = file
- @remove_seqs = true
- @keep_seqs = false
- end
- def set_keep_seqs( file )
- @seqs_name_file = file
- @keep_seqs = true
- @remove_seqs = false
- end
- def set_trim( first, last )
- @trim = true
- @first = first
- @last = last
- end
- def set_remove_matching( remove )
- @remove_matching = remove
- end
- def set_keep_matching( keep )
- @keep_matching = keep
- end
- def set_rem_red( rr )
- @rem_red = rr
- end
-
-
-
- def set_split( s )
- if ( s > 0 )
- @split = s
- @clean = false # phylip only
- @rgc = false
- @rgoc = false
- @rg = false # fasta only
- @rgr = -1
- @rsgr = -1
- @rsl = -1
- @seqs_name_file = nil
- @remove_seqs = false
- @keep_seqs = false
- @trim = false
- @first = -1
- @last = -1
- end
- end
- def analyze_command_line( cla )
- if ( cla.is_option_set?( INPUT_TYPE_OPTION ) )
- begin
- type = cla.get_option_value( INPUT_TYPE_OPTION )
- if ( type == "p" )
- set_phylip_input( true )
- set_fasta_input( false )
- elsif ( type == "f" )
- set_fasta_input( true )
- set_phylip_input( false )
- end
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( OUTPUT_TYPE_OPTION ) )
- begin
- type = cla.get_option_value( OUTPUT_TYPE_OPTION )
- if ( type == "p" )
- set_pi_output( true )
- set_fasta_output( false )
- set_nexus_output( false )
- elsif ( type == "f" )
- set_pi_output( false )
- set_fasta_output( true )
- set_nexus_output( false )
- elsif ( type == "n" )
- set_pi_output( false )
- set_fasta_output( false )
- set_nexus_output( true )
- end
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( MAXIMAL_NAME_LENGTH_OPTION ) )
- begin
- l = cla.get_option_value_as_int( MAXIMAL_NAME_LENGTH_OPTION )
- set_name_length( l )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( WIDTH_OPTION ) )
- begin
- w = cla.get_option_value_as_int( WIDTH_OPTION )
- set_width( w )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( CLEAN_UP_SEQ_OPTION ) )
- set_clean( true )
- end
- if ( cla.is_option_set?( REMOVE_GAP_COLUMNS_OPTION ) )
- set_remove_gap_columns( true )
- end
- if ( cla.is_option_set?( REM_RED_OPTION ) )
- set_rem_red( true )
- end
- if ( cla.is_option_set?( REMOVE_GAP_ONLY_COLUMNS ) )
- set_remove_gap_only_columns( true )
- end
- if ( cla.is_option_set?( REMOVE_ALL_GAP_CHARACTERS_OPTION ) )
- set_remove_gaps( true )
- end
- if ( cla.is_option_set?( REMOVE_COLUMNS_GAP_RATIO_OPTION ) )
- begin
- f = cla.get_option_value_as_float( REMOVE_COLUMNS_GAP_RATIO_OPTION )
- set_remove_gap_ratio( f )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( REMOVE_ALL_SEQUENCES_LISTED_OPTION ) )
- begin
- s = cla.get_option_value( REMOVE_ALL_SEQUENCES_LISTED_OPTION )
- set_remove_seqs( s )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( KEEP_ONLY_SEQUENCES_LISTED_OPTION ) )
- begin
- s = cla.get_option_value( KEEP_ONLY_SEQUENCES_LISTED_OPTION )
- set_keep_seqs( s )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( TRIM_OPTION ) )
- begin
- s = cla.get_option_value( TRIM_OPTION )
- if ( s =~ /(\d+)-(\d+)/ )
- set_trim( $1.to_i(), $2.to_i() )
- else
- puts( "illegal argument" )
- print_help
- exit( -1 )
- end
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( REMOVE_SEQS_GAP_RATIO_OPTION ) )
- begin
- f = cla.get_option_value_as_float( REMOVE_SEQS_GAP_RATIO_OPTION )
- set_remove_seqs_gap_ratio( f )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( REMOVE_SEQS_NON_GAP_LENGTH_OPTION ) )
- begin
- f = cla.get_option_value_as_int( REMOVE_SEQS_NON_GAP_LENGTH_OPTION )
- set_remove_seqs_min_non_gap_length( f )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( SPLIT ) )
- begin
- s = cla.get_option_value_as_int( SPLIT )
- set_split( s )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
-
- end
- if ( cla.is_option_set?( REMOVE_MATCHING_SEQUENCES_OPTION ) )
- begin
- s = cla.get_option_value( REMOVE_MATCHING_SEQUENCES_OPTION )
- set_remove_matching( s )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
- if ( cla.is_option_set?( KEEP_MATCHING_SEQUENCES_OPTION ) )
- begin
- s = cla.get_option_value( KEEP_MATCHING_SEQUENCES_OPTION )
- set_keep_matching( s )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
- end
- end
-
-
- end
-
- def print_help()
- puts()
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb [options] <input alignment> <output>" )
- puts()
- puts( " options: -" + INPUT_TYPE_OPTION + "=<input type>: f for fasta, p for phylip selex type" )
- puts( " -" + OUTPUT_TYPE_OPTION + "=<output type>: f for fasta, n for nexus, p for phylip sequential (default)" )
- puts( " -" + MAXIMAL_NAME_LENGTH_OPTION + "=<n>: n=maximal name length (default for phylip 10, for fasta: unlimited )" )
- puts( " -" + WIDTH_OPTION + "=<n>: n=width (fasta output only, default is 60)" )
- puts( " -" + CLEAN_UP_SEQ_OPTION + ": clean up sequences" )
- puts( " -" + REMOVE_GAP_COLUMNS_OPTION + ": remove gap columns" )
- puts( " -" + REMOVE_GAP_ONLY_COLUMNS + ": remove gap-only columns" )
- puts( " -" + REMOVE_COLUMNS_GAP_RATIO_OPTION + "=<n>: remove columns for which ( seqs with gap / number of sequences > n )" )
- puts( " -" + REMOVE_ALL_GAP_CHARACTERS_OPTION + ": remove all gap characters (destroys alignment, fasta output only)" )
- puts( " -" + REMOVE_ALL_SEQUENCES_LISTED_OPTION + "=<file>: remove all sequences listed in file" )
- puts( " -" + KEEP_ONLY_SEQUENCES_LISTED_OPTION + "=<file>: keep only sequences listed in file" )
- puts( " -" + TRIM_OPTION + "=<first>-<last>: remove columns before first and after last" )
- puts( " -" + REMOVE_SEQS_GAP_RATIO_OPTION + "=<n>: remove sequences for which the gap ratio > n (after column operations)" )
- puts( " -" + REMOVE_SEQS_NON_GAP_LENGTH_OPTION + "=<n> remove sequences with less than n non-gap characters (after column operations)" )
- puts( " -" + REMOVE_MATCHING_SEQUENCES_OPTION + "=<s> remove all sequences with names containing s" )
- puts( " -" + KEEP_MATCHING_SEQUENCES_OPTION + "=<s> keep only sequences with names containing s" )
- puts( " -" + SPLIT + "=<n> split a fasta file into n files of equal number of sequences (expect for " )
- puts( " last one), cannot be used with other options" )
- puts( " -" + REM_RED_OPTION + ": remove redundant sequences" )
- puts()
- end
-
-
-
-
-
- end # class MsaProcessor
-
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/multi_sequence_extractor.rb - MultiSequenceExtractor class
-#
-# Copyright:: Copyright (C) 2006-2008 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: multi_sequence_extractor.rb,v 1.10 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/msa/msa'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/parser/fasta_parser'
-require 'lib/evo/io/writer/fasta_writer'
-require 'lib/evo/util/command_line_arguments'
-
-
-
-module Evoruby
-
- class MultiSequenceExtractor
-
- PRG_NAME = "mse"
- PRG_VERSION = "1.02"
- PRG_DESC = "extraction of sequences by name from multiple multi-sequence ('fasta') files"
- PRG_DATE = "2012.07.20"
- COPYRIGHT = "2008-2012 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
- HELP_OPTION_1 = 'help'
- HELP_OPTION_2 = 'h'
-
- EXT_OPTION = 'e'
- EXTRACT_LINKERS_OPTION = 'l'
- LOG_SUFFIX = ".mse_log"
- LINKERS_SUFFIX = ".linkers"
- FASTA_SUFFIX = ".fasta"
- FASTA_WITH_NORMALIZED_IDS_SUFFIX = ".ni.fasta"
- NORMALIZED_IDS_MAP_SUFFIX = ".nim"
- PROTEINS_LIST_FILE_SEPARATOR = "\t"
-
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC ,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- ld = Constants::LINE_DELIMITER
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if ( cla.get_number_of_files != 4 && cla.get_number_of_files != 5 )
- print_help
- exit( -1 )
- end
-
- allowed_opts = Array.new
- allowed_opts.push(EXT_OPTION)
- allowed_opts.push(EXTRACT_LINKERS_OPTION)
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed,
- STDOUT )
- end
-
- seq_names_files_suffix = cla.get_file_name( 0 )
- input_dir = cla.get_file_name( 1 )
- out_dir = cla.get_file_name( 2 )
- out_dir_doms = cla.get_file_name( 3 )
- mapping_file = nil
-
- if ( cla.get_number_of_files == 5 )
- mapping_file = cla.get_file_name( 4 )
- begin
- Util.check_file_for_readability( mapping_file )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- end
-
- extension = 0
- if cla.is_option_set?(EXT_OPTION)
- extension = cla.get_option_value_as_int(EXT_OPTION)
- if extension < 0
- extension = 0
- end
- end
-
- extract_linkers = false
- if cla.is_option_set?(EXTRACT_LINKERS_OPTION)
- extract_linkers = true
- end
-
- if !File.exist?( input_dir )
- Util.fatal_error( PRG_NAME, "error: input directory [#{input_dir}] does not exist" )
- end
- if !File.exist?( out_dir )
- Util.fatal_error( PRG_NAME, "error: output directory [#{out_dir}] does not exist" )
- end
- if !File.exist?( out_dir_doms )
- Util.fatal_error( PRG_NAME, "error: output directory [#{out_dir_doms}] does not exist" )
- end
- if !File.directory?( input_dir )
- Util.fatal_error( PRG_NAME, "error: [#{input_dir}] is not a directory" )
- end
- if !File.directory?( out_dir )
- Util.fatal_error( PRG_NAME, "error: [#{out_dir}] is not a directory" )
- end
- if !File.directory?( out_dir_doms )
- Util.fatal_error( PRG_NAME, "error: [#{out_dir_doms}] is not a directory" )
- end
-
-
- log = String.new
-
- log << "Program : " + PRG_NAME + ld
- log << "Version : " + PRG_VERSION + ld
- log << "Program date : " + PRG_DATE + ld
-
- puts()
- puts( "Sequence names files suffix: " + seq_names_files_suffix )
- log << "Sequence names files suffix: " + seq_names_files_suffix + ld
- puts( "Input dir : " + input_dir )
- log << "Input dir : " + input_dir + ld
- puts( "Output dir : " + out_dir )
- log << "Output dir : " + out_dir + ld
- puts( "Output dir domains : " + out_dir_doms )
- log << "Output dir domains : " + out_dir_doms + ld
- if ( mapping_file != nil )
- puts( "Mapping file : " + mapping_file )
- log << "Mapping file : " + mapping_file + ld
- end
- if extension > 0
- puts( "Extension : " + extension.to_s )
- log << "Extension : " + extension.to_s + ld
- end
- if extract_linkers
- puts( "Extract linkers : true" )
- log << "Extract linkers : true" + ld
- end
- log << "Date : " + Time.now.to_s + ld
- puts
-
- if ( mapping_file != nil )
- species_codes_to_paths = extract_mappings( mapping_file )
- end
-
- input_files = obtain_inputfiles( input_dir, seq_names_files_suffix )
-
- counter = 0
-
- input_files.each { |input_file|
- counter += 1
- puts
- puts
- puts counter.to_s + "/" + input_files.size.to_s
- read_seq_family_file( input_file,
- seq_names_files_suffix,
- input_dir,
- species_codes_to_paths,
- log,
- out_dir,
- out_dir_doms,
- mapping_file,
- extension,
- extract_linkers )
- }
- puts
- Util.print_message( PRG_NAME, "OK" )
- puts
-
- end
-
-
- def read_seq_family_file( input_file,
- seq_names_files_suffix,
- input_dir,
- species_codes_to_paths,
- log,
- out_dir,
- out_dir_doms,
- mapping_file,
- extension,
- extract_linkers )
-
- begin
- Util.check_file_for_readability( input_file )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- basename = File.basename( input_file, seq_names_files_suffix )
- out_file_path_fasta_file = out_dir + Constants::FILE_SEPARATOR + basename + FASTA_SUFFIX
- out_file_path_normalized_ids_fasta_file = out_dir + Constants::FILE_SEPARATOR + basename + FASTA_WITH_NORMALIZED_IDS_SUFFIX
- out_file_path_ids_map = out_dir + Constants::FILE_SEPARATOR + basename + NORMALIZED_IDS_MAP_SUFFIX
- doms_out_file_path_fasta_file = out_dir_doms + Constants::FILE_SEPARATOR + basename + "_domains" + FASTA_SUFFIX
- doms_ext_out_file_path_fasta_file = nil
- if extension > 0
- doms_ext_out_file_path_fasta_file = out_dir_doms + Constants::FILE_SEPARATOR + basename + "_domains_ext_" + extension.to_s + FASTA_SUFFIX
- end
- begin
- Util.check_file_for_writability( out_file_path_fasta_file )
- Util.check_file_for_writability( out_file_path_normalized_ids_fasta_file )
- Util.check_file_for_writability( out_file_path_ids_map )
- Util.check_file_for_writability( doms_out_file_path_fasta_file )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- ids_map_writer = nil
- begin
- ids_map_writer = File.open( out_file_path_ids_map, 'a' )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- current_species = ""
- current_msa = nil
- new_msa = Msa.new
- new_msa_normalized_ids = Msa.new
- new_msa_domains = Msa.new
- new_msa_domains_extended = Msa.new
- per_species_counter = 0
- linkers = ""
-
- puts basename
-
- File.open( input_file ) do | file |
- while line = file.gets
- line.strip!
- if !Util.is_string_empty?( line ) && !(line =~ /\s*#/ )
- values = line.split( PROTEINS_LIST_FILE_SEPARATOR )
- mod_line = nil
- if ( values.length < 2 )
- Util.fatal_error( PRG_NAME, "unexpected format: " + line )
- end
- species = values[ 0 ]
- if species == "BRADI" || species == "ASPNG" || species == "SCLSC" || species == "PTEVA" || species == "EIMTE"
- next
- end
- seq_name = values[ 1 ]
- domain_ranges = nil
- if ( values.length > 3 )
- domain_ranges_block = values[ 3 ]
- domain_ranges = domain_ranges_block.split( "/" )
- end
- if ( species != current_species )
- current_species = species
- my_file = input_dir + Constants::FILE_SEPARATOR + current_species
-
- if ( !File.exist?( my_file ) )
- if species_codes_to_paths == nil
- Util.fatal_error( PRG_NAME, "error: [#{my_file}] not found and no mapping file provided" )
- elsif ( !species_codes_to_paths.has_key?( current_species ) )
- Util.fatal_error( PRG_NAME, "error: species [#{current_species}] not found in mapping file [#{mapping_file}]" )
- end
- my_file = species_codes_to_paths[ current_species ]
- end
- my_path = File.expand_path( my_file )
- my_readlink = my_path
- if ( File.symlink?( my_path ) )
- my_readlink = File.readlink( my_path )
- end
- current_msa = read_fasta_file( my_file )
-
- if ( per_species_counter > 0 )
- print_counts( per_species_counter, log, Constants::LINE_DELIMITER )
- per_species_counter = 0
- end
- puts " " + current_species + " [" + my_readlink + "]"
- log << current_species + " [" + my_readlink + "]" + Constants::LINE_DELIMITER
- end
- puts " " + seq_name
- log << " " + seq_name + Constants::LINE_DELIMITER
- per_species_counter = per_species_counter + 1
- seq = nil
-
- if current_msa.find_by_name_start( seq_name, true ).size > 0
- begin
- seq = current_msa.get_by_name_start( seq_name, true ).copy
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- else
- # Not found, try finding by partial match.
- begin
- seq = current_msa.get_by_name( seq_name, true, true )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- end
-
- normalized_id = per_species_counter.to_s( 16 ).upcase +
- "_" + current_species
-
- per_species_counter.to_i
-
- ids_map_writer.write( normalized_id + ": " + seq.get_name + Constants::LINE_DELIMITER )
-
- orig_name = nil
- if seq != nil
- orig_name = seq.get_name
- seq.set_name( seq.get_name + " [" + current_species + "]" )
- new_msa.add_sequence( seq )
- else
- Util.fatal_error( PRG_NAME, "unexected error: seq is nil" )
- end
-
- if domain_ranges != nil
- first = true
- prev_to = -1
-
- domain_ranges.each { |range|
- if range != nil && range.length > 0
- s = range.split("-")
- from = s[ 0 ].to_i
- to = s[ 1 ].to_i
- new_msa_domains.add_sequence( Sequence.new( orig_name +
- " [" + from.to_s +
- "-" + to.to_s +
- "] [" + basename + "] [" +
- current_species + "]",
- seq.get_sequence_as_string[from..to] ) )
- if extension > 0
- from_e = from - extension
- if from_e < 0
- from_e = 0
- end
- to_e = to + extension
- if to_e > seq.get_sequence_as_string.length - 1
- to_e = seq.get_sequence_as_string.length - 1
- end
- new_msa_domains_extended.add_sequence( Sequence.new( orig_name +
- " [" + from.to_s +
- "-" + to.to_s +
- "] [extended by " +
- extension.to_s +
- "] [" + basename + "] [" +
- current_species + "]",
- seq.get_sequence_as_string[ from_e..to_e ] ) )
- end # extension > 0
- if extract_linkers
- if first
- first = false
- f = 0
- t = from - 1
- if extension > 0
- f = from - extension
- end
- mod_line = line + "\t[" + get_linker_sequence( f, t, seq ) + "|"
- else
- mod_line += get_linker_sequence( prev_to + 1, from - 1, seq ) + "|"
- end
- prev_to = to
- end
- end # range != nil && range.length > 0
- }
- if extract_linkers && prev_to > 0
- f = prev_to + 1
- t = seq.get_sequence_as_string.length - 1
- if extension > 0
- t = prev_to + extension
- end
- mod_line += get_linker_sequence( f, t, seq ) + "]"
- end
- end
-
- new_msa_normalized_ids.add_sequence( Sequence.new( normalized_id, seq.get_sequence_as_string ) )
- if mod_line
- linkers << mod_line + Constants::LINE_DELIMITER
- end
- end # !Util.is_string_empty?( line ) && !(line =~ /\s*#/ )
- end # while line = file.gets
-
- end
-
- ids_map_writer.close
-
- if ( per_species_counter > 0 )
- print_counts( per_species_counter, log, Constants::LINE_DELIMITER )
- end
-
- io = MsaIO.new()
-
- fasta_writer = FastaWriter.new()
- fasta_writer.remove_gap_chars
- fasta_writer.clean
-
- begin
- io.write_to_file( new_msa, out_file_path_fasta_file, fasta_writer )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if new_msa_domains != nil
- begin
- io.write_to_file( new_msa_domains, doms_out_file_path_fasta_file, fasta_writer )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- end
-
- if extension > 0 && new_msa_domains_extended != nil
- begin
- io.write_to_file( new_msa_domains_extended, doms_ext_out_file_path_fasta_file, fasta_writer )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- end
-
- begin
- io.write_to_file( new_msa_normalized_ids, out_file_path_normalized_ids_fasta_file, fasta_writer )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- begin
- f = File.open( out_dir + Constants::FILE_SEPARATOR + basename + LOG_SUFFIX , 'a' )
- f.print( log )
- f.close
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if extract_linkers && linkers != nil
- begin
- f = File.open( out_dir + Constants::FILE_SEPARATOR + basename + LINKERS_SUFFIX , 'a' )
- f.print( linkers )
- f.close
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- end
-
-
- end
-
-
- def get_linker_sequence( from, to, seq )
- if from < 0
- from = 0
- end
- if to > seq.get_sequence_as_string.length - 1
- to = seq.get_sequence_as_string.length - 1
- end
- if from > to
- return ""
- else
- return from.to_s + "-" + to.to_s + ":" + seq.get_sequence_as_string[ from..to ]
- end
- end
-
- def obtain_inputfiles( input_dir, seq_names_files_suffix )
- input_files = Array.new()
- Dir.foreach( input_dir ) { |file_name|
- if file_name.index( seq_names_files_suffix ) == ( file_name.size - seq_names_files_suffix.size )
- input_files.push( input_dir + Constants::FILE_SEPARATOR + file_name )
- end
- }
- input_files
- end
-
- def extract_mappings( mapping_file )
- species_code_to_path = Hash.new()
- File.open( mapping_file ) do | file |
- while line = file.gets
- if ( !Util.is_string_empty?( line ) && !(line =~ /\s*#/ ) )
- if ( line =~ /(\S+)\s+(\S+)/ )
- species = $1
- path = $2
- if ( species_code_to_path.has_key?( species ) )
- Util.fatal_error( PRG_NAME, "error: species code [#{species}] is not unique" )
- end
- if ( species_code_to_path.has_value?( path ) )
- Util.fatal_error( PRG_NAME, "error: path [#{path}] is not unique" )
- end
- if ( !File.exist?( path ) )
- Util.fatal_error( PRG_NAME, "error: file [#{path}] does not exist" )
- end
- if ( !File.file?( path ) )
- Util.fatal_error( PRG_NAME, "error: [#{path}] is not a regular file" )
- end
- if ( !File.readable?( path ) )
- Util.fatal_error( PRG_NAME, "error: file [#{path}] is not readable" )
- end
- if ( File.size( path ) < 10000 )
- Util.fatal_error( PRG_NAME, "error: file [#{path}] appears too small" )
- end
- if ( !Util.looks_like_fasta?( path ) )
- Util.fatal_error( PRG_NAME, "error: file [#{path}] does not appear to be a fasta file" )
- end
- species_code_to_path[ species ] = path
- puts species + " -> " + path
- end
- end
- end
- end
- species_code_to_path
- end
-
- def print_counts( per_species_counter, log, ld )
- puts " [sum: " + per_species_counter.to_s + "]"
- log << " [sum: " + per_species_counter.to_s + "]" + ld
- end
-
- def read_fasta_file( input )
- f = MsaFactory.new()
- msa = nil
- begin
- msa = f.create_msa_from_file( input, FastaParser.new() )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- msa
- end
-
- def print_help()
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb <sequence names file suffix> <input dir containing sequence names files " +
- "and possibly genome multiple-sequence ('fasta') files> <output directory for sequences> <output directory for domains> [mapping file for " +
- "genome multiple-sequence ('fasta') files not in input dir]" )
- puts()
- puts( " option: -" + EXT_OPTION + "=<int>: to extend extracted domains" )
- puts( " -" + EXTRACT_LINKERS_OPTION + " : to extract linkers" )
- puts()
- puts( " " + "Example: \"mse.rb .prot . protein_seqs domain_seqs ../genome_locations.txt\"" )
- puts()
- end
-
- end # class MultiSequenceExtractor
-end
\ No newline at end of file
+++ /dev/null
-#
-# = lib/evo/apps/ - class
-#
-# Copyright:: Copyright (C) 2009 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: new_tap.rb,v 1.4 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-require 'lib/evo/util/util'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/msa/msa'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/parser/fasta_parser'
-require 'lib/evo/io/parser/general_msa_parser'
-require 'lib/evo/io/writer/fasta_writer'
-require 'lib/evo/io/writer/phylip_sequential_writer'
-require 'lib/evo/util/command_line_arguments'
-
-module Evoruby
-
- class TaxonomyProcessor
-
- PRG_NAME = ""
- PRG_DATE = "2009.10.09"
- PRG_DESC = "replacement of labels in multiple sequence files"
- PRG_VERSION = "1.00"
- COPYRIGHT = "2009 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- REMOVE_REDUNDANT_SEQS_OPTION = "rr"
-
- def initialize()
- @taxonomies = Hash.new()
- end
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if ( ARGV == nil || ( ARGV.length != 3 && ARGV.length != 4 ) )
- puts( "Usage: #{PRG_NAME}.rb <input sequences> <output sequences> <output map>" )
- puts()
- puts( " options: -" + REMOVE_REDUNDANT_SEQS_OPTION + ": to remove redundant sequences" )
- puts()
- exit( -1 )
- end
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- input = cla.get_file_name( 0 )
- output = cla.get_file_name( 1 )
- map_file = cla.get_file_name( 2 )
-
- allowed_opts = Array.new
- allowed_opts.push( REMOVE_REDUNDANT_SEQS_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
- end
-
-
- remove_redudant = false
- if ( cla.is_option_set?( REMOVE_REDUNDANT_SEQS_OPTION ) )
- remove_redudant = true
- end
-
- if ( File.exists?( output ) )
- Util.fatal_error( PRG_NAME, "outfile [" + output + "] already exists" )
- end
- if ( File.exists?( map_file ) )
- Util.fatal_error( PRG_NAME, "map file [" + map_file + "] already exists" )
- end
- if ( !File.exists?( input) )
- Util.fatal_error( PRG_NAME, "infile [" + input + "] does not exist" )
- end
-
- fasta_like = Util.looks_like_fasta?( input )
-
- puts()
- puts( "Input alignment : " + input )
- puts( "Output alignment: " + output )
- puts( "Output map : " + map_file )
- if ( fasta_like )
- puts( "Format : Fasta" )
- else
- puts( "Format : Phylip like" )
- end
- puts()
-
- species_map = Hash.new
-
- f = MsaFactory.new()
- begin
- if ( fasta_like )
- msa = f.create_msa_from_file( input, FastaParser.new() )
- else
- msa = f.create_msa_from_file( input, GeneralMsaParser.new() )
- end
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "failed to read file: " + e.to_s )
- end
-
- if ( msa == nil || msa.get_number_of_seqs() < 1 )
- Util.fatal_error( PRG_NAME, "failed to read MSA" )
- end
- begin
- Util.check_file_for_writability( map_file )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_, STDOUT )
- end
-
- if ( remove_redudant )
- removed = msa.remove_redundant_sequences!( true )
- if removed.size > 0
- Util.print_message( PRG_NAME, "going to ignore the following " + removed.size.to_s + " redundant sequences:" )
- removed.each { | seq_name |
- puts seq_name
- }
- Util.print_message( PRG_NAME, "will process " + msa.get_number_of_seqs.to_s + " non redundant sequences" )
- end
- end
-
- lf = File.open( map_file, "a" )
- for i in 0 ... msa.get_number_of_seqs
- seq = msa.get_sequence( i )
- end
-
- io = MsaIO.new()
- w = nil
- if ( fasta_like )
- w = FastaWriter.new()
- else
- w = PhylipSequentialWriter.new()
- end
- w.set_max_name_length( 10 )
- w.clean( true )
- begin
- io.write_to_file( msa, output, w )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
- end
- lf.close()
- if ( @taxonomies.length > 0 )
- Util.print_message( PRG_NAME, "number of unique taxonomies: " + @taxonomies.length.to_s )
- end
- Util.print_message( PRG_NAME, "wrote: " + map_file )
- Util.print_message( PRG_NAME, "wrote: " + output )
- Util.print_message( PRG_NAME, "OK" )
- end
-
- end # class
-
-end # module Evoruby
\ No newline at end of file
+++ /dev/null
-#!/usr/local/bin/ruby -w
-#
-# = lib/evo/apps/phylogenies_decorator
-#
-# Copyright:: Copyright (C) 2006-2008 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# decoration of phylogenies with sequence/species names and domain architectures
-#
-# $Id: phylogenies_decorator.rb,v 1.34 2010/12/13 19:00:11 cmzmasek Exp $
-#
-# Environment variable FORESTER_HOME needs to point to the appropriate
-# directory (e.g. setenv FORESTER_HOME $HOME/SOFTWARE_DEV/ECLIPSE_WORKSPACE/forester-atv/)
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-
-require 'date'
-
-module Evoruby
-
- class PhylogeniesDecorator
-
- DECORATOR_OPTIONS_SEQ_NAMES = '-r=1 -mdn'
- # -mdn is a hidden expert option to rename e.g. "6_ORYLA3" to "6_[3]_ORYLA"
- #DECORATOR_OPTIONS_SEQ_NAMES = '-sn -r=1'
- DECORATOR_OPTIONS_DOMAINS = '-r=1'
- IDS_MAPFILE_SUFFIX = '.nim'
- DOMAINS_MAPFILE_SUFFIX = '.dff'
- SLEEP_TIME = 0.1
- REMOVE_NI = true
- TMP_FILE = '___PD___'
- LOG_FILE = '00_phylogenies_decorator.log'
- FORESTER_HOME = ENV[Constants::FORESTER_HOME_ENV_VARIABLE]
- JAVA_HOME = ENV[Constants::JAVA_HOME_ENV_VARIABLE]
-
- PRG_NAME = "phylogenies_decorator"
- PRG_DATE = "2008.09.02"
- PRG_DESC = "decoration of phylogenies with sequence/species names and domain architectures"
- PRG_VERSION = "1.0.1"
- COPYRIGHT = "2008-2009 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- IDS_ONLY_OPTION = "n"
- DOMAINS_ONLY_OPTION = "d"
- HELP_OPTION_1 = "help"
- HELP_OPTION_2 = "h"
-
- NL = Constants::LINE_DELIMITER
-
- def run
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if ( ARGV == nil || ARGV.length > 3 || ARGV.length < 2 )
- print_help
- exit( -1 )
- end
-
- if FORESTER_HOME == nil || FORESTER_HOME.length < 1
- Util.fatal_error( PRG_NAME, "apparently environment variable #{Constants::FORESTER_HOME_ENV_VARIABLE} has not been set" )
- end
- if JAVA_HOME == nil || JAVA_HOME.length < 1
- Util.fatal_error( PRG_NAME, "apparently environment variable #{Constants::JAVA_HOME_ENV_VARIABLE} has not been set" )
- end
-
- if !File.exist?( FORESTER_HOME )
- Util.fatal_error( PRG_NAME, '[' + FORESTER_HOME + '] does not exist' )
- end
- if !File.exist?( JAVA_HOME )
- Util.fatal_error( PRG_NAME, '[' + JAVA_HOME + '] does not exist' )
- end
-
- decorator = JAVA_HOME + '/bin/java -cp ' + FORESTER_HOME + '/java/forester.jar org.forester.application.decorator'
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- if ( cla.is_option_set?( HELP_OPTION_1 ) ||
- cla.is_option_set?( HELP_OPTION_2 ) )
- print_help
- exit( 0 )
- end
-
- if File.exist?( LOG_FILE )
- Util.fatal_error( PRG_NAME, 'logfile [' + LOG_FILE + '] already exists' )
- end
-
- allowed_opts = Array.new
- allowed_opts.push( IDS_ONLY_OPTION )
- allowed_opts.push( DOMAINS_ONLY_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
- end
-
- ids_only = false
- domains_only = false
-
- in_suffix = cla.get_file_name( 0 )
- out_suffix = cla.get_file_name( 1 )
-
- if cla.is_option_set?( IDS_ONLY_OPTION )
- ids_only = true
- end
- if cla.is_option_set?( DOMAINS_ONLY_OPTION )
- domains_only = true
- end
-
- if ( ids_only && domains_only )
- Util.fatal_error( PRG_NAME, 'attempt to use ids only and domains only at the same time' )
- end
-
- log = String.new
-
- now = DateTime.now
- log << "Program : " + PRG_NAME + NL
- log << "Version : " + PRG_VERSION + NL
- log << "Program date : " + PRG_DATE + NL
- log << "Options for seq names: " + DECORATOR_OPTIONS_SEQ_NAMES + NL
- log << "Options for domains : " + DECORATOR_OPTIONS_DOMAINS + NL
- log << "FORESTER_HOME : " + FORESTER_HOME + NL
- log << "JAVA_HOME : " + JAVA_HOME + NL + NL
- log << "Date/time: " + now.to_s + NL
- log << "Directory: " + Dir.getwd + NL + NL
-
- Util.print_message( PRG_NAME, 'input suffix : ' + in_suffix )
- Util.print_message( PRG_NAME, 'output suffix : ' + out_suffix )
-
- log << 'input suffix : ' + in_suffix + NL
- log << 'output suffix : ' + out_suffix + NL
-
- if ( File.exists?( TMP_FILE ) )
- File.delete( TMP_FILE )
- end
-
- files = Dir.entries( "." )
-
- counter = 0
-
- files.each { | phylogeny_file |
- if ( !File.directory?( phylogeny_file ) &&
- phylogeny_file !~ /^\./ &&
- phylogeny_file !~ /^00/ &&
- phylogeny_file !~ /#{out_suffix}$/ &&
- phylogeny_file =~ /#{in_suffix}$/ )
- begin
- Util.check_file_for_readability( phylogeny_file )
- rescue ArgumentError
- Util.fatal_error( PRG_NAME, 'can not read from: ' + phylogeny_file + ': '+ $! )
- end
-
- counter += 1
-
- outfile = phylogeny_file.sub( /#{in_suffix}$/, out_suffix )
-
- if REMOVE_NI
- outfile = outfile.sub( /_ni_/, '_' )
- end
-
- if File.exists?( outfile )
- msg = counter.to_s + ': ' + phylogeny_file + ' -> ' + outfile +
- ' : already exists, skipping'
- Util.print_message( PRG_NAME, msg )
- log << msg + NL
- next
- end
-
- Util.print_message( PRG_NAME, counter.to_s + ': ' + phylogeny_file + ' -> ' + outfile )
- log << counter.to_s + ': ' + phylogeny_file + ' -> ' + outfile + NL
-
- phylogeny_id = get_id( phylogeny_file )
-
- ids_mapfile_name = nil
- domains_mapfile_name = nil
-
- if ids_only
- ids_mapfile_name = get_file( files, phylogeny_id, IDS_MAPFILE_SUFFIX )
- elsif domains_only
- domains_mapfile_name = get_file( files, phylogeny_id, DOMAINS_MAPFILE_SUFFIX )
- else
- ids_mapfile_name = get_file( files, phylogeny_id, IDS_MAPFILE_SUFFIX )
- domains_mapfile_name = get_file( files, phylogeny_id, DOMAINS_MAPFILE_SUFFIX )
- end
-
- if domains_mapfile_name != nil
- begin
- Util.check_file_for_readability( domains_mapfile_name )
- rescue ArgumentError
- Util.fatal_error( PRG_NAME, 'failed to read from [#{domains_mapfile_name}]: ' + $! )
- end
- end
-
- if ids_mapfile_name != nil
- begin
- Util.check_file_for_readability( ids_mapfile_name )
- rescue ArgumentError
- Util.fatal_error( PRG_NAME, 'failed to read from [#{ids_mapfile_name}]: ' + $! )
- end
- end
-
- if domains_mapfile_name != nil
- if ids_mapfile_name != nil
- my_outfile = TMP_FILE
- else
- my_outfile = outfile
- end
- cmd = decorator + ' ' + DECORATOR_OPTIONS_DOMAINS + ' ' +
- '-f=d ' + phylogeny_file + ' ' +
- domains_mapfile_name + ' ' + my_outfile
- execute_cmd( cmd, log )
- end
-
- if ids_mapfile_name != nil
- if domains_mapfile_name != nil
- my_infile = TMP_FILE
- else
- my_infile = phylogeny_file
- end
- cmd = decorator + ' ' + DECORATOR_OPTIONS_SEQ_NAMES + ' ' +
- '-f=s ' + my_infile + ' ' +
- ids_mapfile_name + ' ' + outfile
- execute_cmd( cmd, log )
- end
-
- if ( File.exists?( TMP_FILE ) )
- File.delete( TMP_FILE )
- end
- end
- }
- open( LOG_FILE, 'w' ) do | f |
- f.write( log )
- end
- puts
- Util.print_message( PRG_NAME, 'OK' )
- puts
- end # def run
-
- def execute_cmd( cmd, log )
- log << 'excuting ' + cmd + NL
- IO.popen( cmd , 'r+' ) do | pipe |
- pipe.close_write
- log << pipe.read + NL + NL
- end
- sleep( SLEEP_TIME )
- end
-
-
- def get_id( phylogeny_file_name )
- phylogeny_file_name =~ /^([^_]+)/
- $1
- end
-
- def get_file( files_in_dir, phylogeny_id, suffix_pattern )
- matching_files = Array.new
- files_in_dir.each { | file |
-
- if ( !File.directory?( file ) &&
- file !~ /^\./ &&
- file !~ /^00/ &&
- file =~ /^#{phylogeny_id}.*#{suffix_pattern}$/ )
- matching_files << file
- end
- }
- if matching_files.length < 1
- Util.fatal_error( PRG_NAME, 'no file matching [' + phylogeny_id +
- '_] [' + suffix_pattern + '] present in current directory' )
- elsif matching_files.length > 1
- Util.fatal_error( PRG_NAME, 'more than one file matching [' + phylogeny_id +
- '_] [' + suffix_pattern + '] present in current directory' )
- end
- matching_files[ 0 ]
- end
-
- def print_help()
- puts( "Usage:" )
- puts()
- puts( " " + PRG_NAME + ".rb [options] <suffix of intrees to be decorated> <suffix for decorated outtrees> " )
- puts()
- puts( " options: -" + IDS_ONLY_OPTION + ": decorate with sequence/species names only" )
- puts( " -" + DOMAINS_ONLY_OPTION + ": decorate with domain structures" )
- puts()
- end
- end # class PhylogenyiesDecorator
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/phylogeny_factory - PhylogenyFactory class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: phylogeny_factory.rb,v 1.32 2010/12/13 19:00:11 cmzmasek Exp $
-
-require 'lib/evo/util/constants'
-require 'lib/evo/util/util'
-require 'lib/evo/util/command_line_arguments'
-
-require 'set'
-require 'date'
-
-module Evoruby
-
- class PhylogenyFactory
-
- PRG_NAME = "phylogeny_factory"
- PRG_DATE = "2010.05.26"
- PRG_DESC = "automated phylogeny reconstruction using queing system"
- PRG_VERSION = "1.1"
- COPYRIGHT = "2010 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- USE_JOB_SUBMISSION_SYSTEM_OPTION = 's'
- LOG_FILE = '00_phylogeny_factory.log'
- TEMPLATE_FILE = '00_phylogeny_factory.template'
- PBS_O_WORKDIR = '$PBS_O_WORKDIR/'
- MIN_LENGTH_DEFAULT = 40
- WALLTIME = '100:00:00'
- QUEUE = 'default'
-
- TMP_CMD_FILE_SUFFIX = '_QSUB'
-
- HMM = 'HMM'
- RSL = 'RSL'
-
- OPTION_OPEN = '%['
- OPTION_CLOSE = ']%'
-
- WAIT = 1.0
-
- NL = Constants::LINE_DELIMITER
-
- def run
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
- allowed_opts = Array.new
- allowed_opts.push( USE_JOB_SUBMISSION_SYSTEM_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME,
- "unknown option(s): " + disallowed,
- STDOUT )
- end
-
- if File.exists?( LOG_FILE )
- puts( '[' + PRG_NAME + '] > log file [' + LOG_FILE + '] already exists' )
- exit( -1 )
- end
-
- if !File.exists?( TEMPLATE_FILE )
- puts( '[' + PRG_NAME + '] > template file [' + TEMPLATE_FILE + '] not found' )
- exit( -1 )
- end
-
- use_job_submission_system = false
- if cla.is_option_set?( USE_JOB_SUBMISSION_SYSTEM_OPTION )
- use_job_submission_system = true
- end
-
- log = String.new
-
- now = DateTime.now
- log << "Program : " + PRG_NAME + NL
- log << "Version : " + PRG_VERSION + NL
- log << "Program date: " + PRG_DATE + NL + NL
- log << "Date/time : " + now.to_s + NL
- log << "Directory : " + Dir.getwd + NL + NL
-
- puts( '[' + PRG_NAME + '] > reading ' + TEMPLATE_FILE )
-
- paths = Hash.new # path placeholder -> full path
- min_lengths = Hash.new # alignment id -> minimal length
- hmms = Hash.new # alignment id -> hmm
- options = Hash.new # option placeholder -> option
- ids = Set.new
-
- commands = Array.new
-
- log << "////////////////////////////////////////////////////////////////// #{NL}"
- log << "Template file [" + TEMPLATE_FILE + "]:#{NL}"
-
- command = String.new
-
- open( TEMPLATE_FILE ).each { | line |
- log << line
- if ( line =~ /^#/ )
- elsif ( line =~ /^\$\s*(\S+)\s*=\s*(\S+)/ )
- paths[ $1 ] = $2
- puts( '[' + PRG_NAME + '] > paths : ' + $1 + ' => ' + $2 )
-
- elsif ( line =~ /^%\s*#{HMM}\s*(\S+)\s*=\s*(\S+)/ )
- hmms[ $1 ] = $2
- puts( '[' + PRG_NAME + '] > hmms : ' + $1 + ' => ' + $2 )
-
- elsif ( line =~ /^%\s*#{RSL}\s*(\S+)\s*=\s*(\S+)/ )
- min_lengths[ $1 ] = $2
- puts( '[' + PRG_NAME + '] > min lengths: ' + $1 + ' => ' + $2 )
-
- elsif ( line =~ /^%\s*(\S+)\s*=\s*(\S+)/ )
- options[ $1 ] = $2
- puts( '[' + PRG_NAME + '] > options : ' + $1 + ' => ' + $2 )
-
- elsif ( line =~ /^>\s*(.+)/ )
- command = command + $1 + ";#{NL}"
-
- elsif ( line =~ /^-/ )
- commands << prepare( command, paths )
- command = String.new
- end
- }
- log << "\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ #{NL}#{NL}"
-
- files = Dir.entries( "." )
-
- files.each { | file |
- if ( !File.directory?( file ) &&
- file !~ /^\./ &&
- file !~ /#{TEMPLATE_FILE}/ &&
- file !~ /.bck$/ &&
- file !~ /.log$/ &&
- file !~ /nohup/ &&
- file !~ /^00/ )
- aln_name = file.to_str
- id = get_id( aln_name )
- if !ids.include?( id )
- ids.add( id )
- end
- puts( '[' + PRG_NAME + '] > file [id] : ' + aln_name + ' [' + id + ']' )
- commands.each do | cmd |
-
- cmd = subst_hmm( cmd, aln_name, hmms )
- cmd = subst_min_length( cmd, aln_name, min_lengths )
- cmd = subst_options( cmd, options )
- if use_job_submission_system
- cmd = subst_aln_name( cmd, PBS_O_WORKDIR + aln_name )
- else
- cmd = subst_aln_name( cmd, aln_name )
- end
-
- if ( cmd =~ /%/ )
- cmd =~ /(%.*?%)/
- problem = $1
- puts( '[' + PRG_NAME + '] > WARNING : [' + id + '] command still contains placeholder: ' + problem )
- log << "WARNING: command still contains placeholder: " + cmd + NL
- else
- tmp_cmd_file = file.to_str[ 0..4 ] + TMP_CMD_FILE_SUFFIX
- if ( File.exists?( tmp_cmd_file ) )
- File.delete( tmp_cmd_file )
- end
- if use_job_submission_system
- open( tmp_cmd_file, 'w' ) do |f|
- f.write( cmd )
- end
- end
-
- log << cmd + NL
-
- if use_job_submission_system
- IO.popen( 'qsub -q ' + QUEUE + ' -l walltime=' + WALLTIME + ' ' + tmp_cmd_file , 'r+' ) do | pipe |
- pipe.close_write
- end
- else
- spawn( 'nohup ' + cmd + ' &', STDERR => "/dev/null" )
- end
-
- sleep( WAIT )
- if ( File.exists?( tmp_cmd_file ) )
- File.delete( tmp_cmd_file )
- end
- end
- end
- end
- }
-
- open( LOG_FILE, 'w' ) do | f |
- f.write( log )
- end
-
- puts()
- puts( '[' + PRG_NAME + '] > OK' )
- puts()
-
- end # def run
-
- def prepare( command, paths )
- paths.each_pair{ | name, full |
- command = command.gsub( name, full )
- }
- command
- end
-
- def subst_options( command, options )
- opt_placeholders = command.scan( /%\[\S+\]%/ )
- opt_placeholders.each { | opt_placeholder |
- opt_placeholder = opt_placeholder.gsub( OPTION_OPEN , '' )
- opt_placeholder = opt_placeholder.gsub( OPTION_CLOSE, '' )
- opt_value = options[ opt_placeholder ]
- if ( opt_value != nil && opt_value.size > 0 )
- command = command.gsub( OPTION_OPEN + opt_placeholder + OPTION_CLOSE, opt_value )
- end
- }
- command
- end
-
- def subst_aln_name( command, aln_name )
- command = command.gsub( '$', aln_name )
- command
- end
-
- def subst_hmm( command, aln_name, hmms )
- id = get_id( aln_name )
- hmm = hmms[ id ]
- if ( hmm != nil && hmm.size > 0 )
- command = command.gsub( OPTION_OPEN + HMM + OPTION_CLOSE, hmm )
- end
- command
- end
-
- def subst_min_length( command, aln_name, min_lengths )
- id = get_id( aln_name )
- min_length = min_lengths[ id ]
- if ( min_length != nil && min_length.size > 0 )
- command = command.gsub( OPTION_OPEN + RSL + OPTION_CLOSE, min_length )
- else
- command = command.gsub( OPTION_OPEN + RSL + OPTION_CLOSE, MIN_LENGTH_DEFAULT.to_s )
- end
- command
- end
-
- def get_id( aln_name )
- aln_name =~ /^([^_]+)/
- $1
- end
-
- end # class PhylogenyFactory
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/taxonomy_processor - TaxonomyProcessor class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: taxonomy_processor.rb,v 1.26 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-require 'lib/evo/util/util'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/msa/msa'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/parser/fasta_parser'
-require 'lib/evo/io/parser/general_msa_parser'
-require 'lib/evo/io/writer/fasta_writer'
-require 'lib/evo/io/writer/phylip_sequential_writer'
-require 'lib/evo/util/command_line_arguments'
-
-module Evoruby
-
- class TaxonomyProcessor
-
- PRG_NAME = "tap"
- PRG_DATE = "2012.09.27"
- PRG_DESC = "replacement of species names in multiple sequence files"
- PRG_VERSION = "1.02"
- COPYRIGHT = "2012 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- SIMPLE = true
-
- EXTRACT_TAXONOMY_OPTION = "t"
-
- def initialize()
- @taxonomies = Hash.new()
- end
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if ( ARGV == nil || ( ARGV.length != 1 && ARGV.length != 3 && ARGV.length != 4 && ARGV.length != 5 && ARGV.length != 6 ) )
- puts( "Usage: #{PRG_NAME}.rb [options] [input map file] <input sequences> [output sequences] [output id list]" )
- puts()
- puts( " options: -" + EXTRACT_TAXONOMY_OPTION + ": to extract taxonomy information from bracketed expression" )
- puts()
- exit( -1 )
- end
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
-
-
- mapfile = nil
- input = nil
- output = nil
- list_file = nil
-
-
-
- if cla.get_number_of_files == 4
- mapfile = cla.get_file_name( 0 )
- input = cla.get_file_name( 1 )
- output = cla.get_file_name( 2 )
- list_file = cla.get_file_name( 3 )
- elsif cla.get_number_of_files == 3
- input = cla.get_file_name( 0 )
- output = cla.get_file_name( 1 )
- list_file = cla.get_file_name( 2 )
- elsif cla.get_number_of_files == 1
- input = cla.get_file_name( 0 )
- i = nil
- if input.downcase.end_with?( ".fasta" )
- i = input[ 0 .. input.length - 7 ]
- elsif input.downcase.end_with?( ".fsa" )
- i = input[ 0 .. input.length - 5 ]
- else
- i = input
- end
- output = i + "_ni.fasta"
- list_file = i + ".nim"
- end
-
-
- allowed_opts = Array.new
- allowed_opts.push( EXTRACT_TAXONOMY_OPTION )
-
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
- end
-
- extract_taxonomy = false
- if ( cla.is_option_set?( EXTRACT_TAXONOMY_OPTION ) )
- extract_taxonomy = true
- end
-
- if ( File.exists?( output ) )
- Util.fatal_error( PRG_NAME, "outfile [" + output + "] already exists" )
- end
- if ( File.exists?( list_file ) )
- Util.fatal_error( PRG_NAME, "list file [" + list_file + "] already exists" )
- end
- if ( !File.exists?( input) )
- Util.fatal_error( PRG_NAME, "infile [" + input + "] does not exist" )
- end
- if ( mapfile != nil && !File.exists?( mapfile ) )
- Util.fatal_error( PRG_NAME, "mapfile [" + mapfile + "] does not exist" )
- end
-
- fasta_like = Util.looks_like_fasta?( input )
-
- puts()
- if mapfile != nil
- puts( "Map file : " + mapfile )
- end
- puts( "Input alignment : " + input )
- puts( "Output alignment: " + output )
- puts( "Name list : " + list_file )
- if ( fasta_like )
- puts( "Format : Fasta" )
- else
- puts( "Format : Phylip like" )
- end
- if ( extract_taxonomy )
- puts( "Extract taxonomy: true" )
- end
- puts()
-
- species_map = Hash.new
- if mapfile != nil
- File.open( mapfile ) do | file |
- while line = file.gets
- if ( line =~/(.+)#(.+)/ || line =~/(.+)\s+(.+)/ )
- species_map[ $1 ] = $2
- Util.print_message( PRG_NAME, "mapping: " + $1 + ' => ' + $2 )
- end
- end
- end
- end
-
- f = MsaFactory.new()
- begin
- if ( fasta_like )
- msa = f.create_msa_from_file( input, FastaParser.new() )
- else
- msa = f.create_msa_from_file( input, GeneralMsaParser.new() )
- end
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "failed to read file: " + e.to_s )
- end
-
- if ( msa == nil || msa.get_number_of_seqs() < 1 )
- Util.fatal_error( PRG_NAME, "failed to read MSA" )
- end
- begin
- Util.check_file_for_writability( list_file )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_, STDOUT )
- end
-
- #removed = msa.remove_redundant_sequences!( true )
- #if removed.size > 0
- # Util.print_message( PRG_NAME, "going to ignore the following " + removed.size.to_s + " redundant sequences:" )
- # removed.each { | seq_name |
- # puts seq_name
- # }
- # Util.print_message( PRG_NAME, "will process " + msa.get_number_of_seqs.to_s + " non redundant sequences" )
- #end
-
- lf = File.open( list_file, "a" )
- for i in 0 ... msa.get_number_of_seqs
- seq = msa.get_sequence( i )
- seq.set_name( Util::normalize_seq_name( modify_name( seq.get_name(), i, lf, species_map, extract_taxonomy ), 10 ) )
- end
-
- io = MsaIO.new()
- w = nil
- if ( fasta_like )
- w = FastaWriter.new()
- else
- w = PhylipSequentialWriter.new()
- end
- w.set_max_name_length( 10 )
- w.clean( true )
- begin
- io.write_to_file( msa, output, w )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
- end
- lf.close()
- if ( @taxonomies.length > 0 )
- Util.print_message( PRG_NAME, "number of unique taxonomies: " + @taxonomies.length.to_s )
- end
- Util.print_message( PRG_NAME, "wrote: " + list_file )
- Util.print_message( PRG_NAME, "wrote: " + output )
- Util.print_message( PRG_NAME, "OK" )
- end
-
- private
-
- def modify_name( desc, counter, file, species_map, extract_taxonomy )
- new_desc = nil
- my_species = nil
- # if desc =~ /^>?\s*\S{1,10}_([0-9A-Z]{3,5})/
- if desc =~ /^>?\s*\S{1,10}_([A-Z]{3,5})/
- new_desc = counter.to_s( 16 ) + "_" + $1
- elsif SIMPLE
- new_desc = counter.to_s( 16 )
- elsif extract_taxonomy
- if ( desc.count( "[" ) != desc.count( "]" ) )
- Util.fatal_error( PRG_NAME, "illegal bracket count in: " + desc )
- end
- species = nil
- species_map.each_key do | key |
- if desc =~ /[\b|_]#{key}\b/ # Added boundaries to prevent e.g. RAT matching ARATH.
- species = species_map[ key ]
- new_desc = counter.to_s( 16 ) + "_" + species
- break
- end
- end
- if species == nil
- if desc =~/.*\[(\S{3,}?)\]/
- species = $1
- species.strip!
- species.upcase!
- species.gsub!( /\s+/, " " )
- species.gsub!( /-/, "" )
- species.gsub!( /\)/, "" )
- species.gsub!( /\(/, "" )
- species.gsub!( /\'/, "" )
- if species =~ /\S+\s\S+/ || species =~ /\S{3,5}/
- if species =~ /(\S+)\s(\S+)/
- code = $1[ 0..2 ] + $2[ 0..1 ]
- elsif species =~ /\S{3,5}/
- code = species
- elsif species.count( " " ) > 2
- species =~ /(\S+)\s+(\S+)\s+(\S+)$/
- third_last = $1
- second_last = $2
- last = $3
- code = code[ 0 ] + third_last[ 0 ] + second_last[ 0 ] + last[ 0 ] + last[ last.size - 1 ]
- elsif species.count( " " ) > 1
- species =~ /(\S+)\s+(\S+)$/
- second_last = $1
- last = $2
- code = code[ 0..1 ] + second_last[ 0 ] + last[ 0 ] + last[ last.size - 1 ]
- end
- new_desc = counter.to_s( 16 ) + "_" + code
- if @taxonomies.has_key?( code )
- if ( !@taxonomies.has_value?( species ) )
- Util.fatal_error( PRG_NAME, "code [#{code}] is not unique in [#{desc}]" )
- end
- else
- if ( @taxonomies.has_value?( species ) )
- Util.fatal_error( PRG_NAME, "genome [#{species}] is not unique in [#{desc}]" )
- else
- @taxonomies[ code ] = species
- end
- end
- else
- Util.fatal_error( PRG_NAME, "illegal format [#{species}] in: " + desc )
- end
- else
- Util.fatal_error( PRG_NAME, "illegal format in: " + desc )
- end
- end
- else
- species = nil
- my_species = nil
- species_map.each_key do | key |
- if desc =~ /#{key}/
- species = species_map[ key ]
- species = species.gsub( /\s+/, "" )
- species = species.gsub( /_/, " " )
- my_species = species
- if species =~ /(\S+)\s+(\S+)/
- species = $1[0..2] + $2[0..1]
- end
- species = species.gsub( /\s+/, "" )
- species = species.slice(0, 5)
- species.upcase!
- break
- end
- end
- if species == nil
- Util.fatal_error( PRG_NAME, "species not found in: " + desc )
- else
- new_desc = counter.to_s( 16 ) + "_" + species
- end
- end
- if new_desc == nil
- Util.fatal_error( PRG_NAME, "failed to extract species from: " + desc )
- end
- if my_species != nil
- file.print( new_desc + ": " + desc + " [" + my_species + "]" + "\n" )
- else
- file.print( new_desc + ": " + desc + "\n" )
- end
- new_desc
- end
-
- end # class TaxonomyProcessor
-
-end # module Evoruby
+++ /dev/null
-#
-# = lib/evo/apps/tseq_taxonomy_processor - TseqTaxonomyProcessor class
-#
-# Copyright:: Copyright (C) 2006-2007 Christian M. Zmasek
-# License:: GNU Lesser General Public License (LGPL)
-#
-# $Id: tseq_taxonomy_processor.rb,v 1.6 2010/12/13 19:00:11 cmzmasek Exp $
-
-
-require 'lib/evo/util/util'
-require 'lib/evo/msa/msa_factory'
-require 'lib/evo/msa/msa'
-require 'lib/evo/io/msa_io'
-require 'lib/evo/io/parser/sp_taxonomy_parser'
-require 'lib/evo/io/parser/ncbi_tseq_parser'
-require 'lib/evo/io/writer/fasta_writer'
-require 'lib/evo/io/writer/phylip_sequential_writer'
-require 'lib/evo/util/command_line_arguments'
-
-module Evoruby
-
- class TseqTaxonomyProcessor
-
- PRG_NAME = "tseq_tap"
- PRG_DATE = "2009.01.06"
- PRG_DESC = "preprocessing of multiple sequence files in ncbi tseq xml format"
- PRG_VERSION = "1.02"
- COPYRIGHT = "2009 Christian M Zmasek"
- CONTACT = "phylosoft@gmail.com"
- WWW = "www.phylosoft.org"
-
- TAXONOMY_CODE = "TAXONOMY_CODE:"
- TAXONOMY_ID = "TAXONOMY_ID:"
- TAXONOMY_ID_TYPE = "TAXONOMY_ID_TYPE:"
- TAXONOMY_SN = "TAXONOMY_SN:"
- TAXONOMY_CN = "TAXONOMY_CN:"
- SEQ_ACCESSION = "SEQ_ACCESSION:"
- SEQ_ACCESSION_SOURCE = "SEQ_ACCESSION_SOURCE:"
- SEQ_SECONDARY_ACCESSION = "SEQ_SECONDARY_ACCESSION:"
- SEQ_SYMBOL = "SEQ_SYMBOL:"
- SEQ_NAME = "SEQ_NAME:"
- SEQ_MOL_SEQ = "SEQ_MOL_SEQ:"
-
- def initialize()
- @tax_ids_to_sp_taxonomies = Hash.new()
- end
-
- def run()
-
- Util.print_program_information( PRG_NAME,
- PRG_VERSION,
- PRG_DESC,
- PRG_DATE,
- COPYRIGHT,
- CONTACT,
- WWW,
- STDOUT )
-
- if ARGV == nil || ARGV.length != 4
- puts( "Usage: #{PRG_NAME}.rb <sp taxonomy file> <sequences in tseq xml format> <name for fasta outfile> <name for map outfile>" )
- puts()
-
- exit( -1 )
- end
-
- begin
- cla = CommandLineArguments.new( ARGV )
- rescue ArgumentError => e
- Util.fatal_error( PRG_NAME, "error: " + e.to_s )
- end
- allowed_opts = Array.new
- disallowed = cla.validate_allowed_options_as_str( allowed_opts )
- if ( disallowed.length > 0 )
- Util.fatal_error( PRG_NAME, "unknown option(s): " + disallowed )
- end
-
- sp_taxonomy_infile = cla.get_file_name( 0 )
- sequences_infile = cla.get_file_name( 1 )
- sequences_outfile = cla.get_file_name( 2 )
- mapping_outfile = cla.get_file_name( 3 )
-
- Util.fatal_error_if_not_readable( PRG_NAME, sp_taxonomy_infile )
- Util.fatal_error_if_not_readable( PRG_NAME, sequences_infile )
- Util.fatal_error_if_not_writable( PRG_NAME, mapping_outfile )
- Util.fatal_error_if_not_writable( PRG_NAME, sequences_outfile )
-
- sp_taxonomies = SpTaxonomyParser.parse( sp_taxonomy_infile )
-
- Util.print_message( PRG_NAME, "read in taxonomic data for " + sp_taxonomies.size.to_s + " species from: " + sp_taxonomy_infile )
-
- tseq_parser = NcbiTSeqParser.new
- msa_fac = MsaFactory.new
-
- seqs = msa_fac.create_msa_from_file( sequences_infile, tseq_parser )
-
- Util.print_message( PRG_NAME, "read in " + seqs.get_number_of_seqs.to_s + " sequences from: " + sequences_infile )
-
- removed = seqs.remove_redundant_sequences!( true, true )
-
- if removed.size > 0
- Util.print_message( PRG_NAME, "going to ignore the following " + removed.size.to_s + " redundant sequences:" )
- removed.each { | seq_name |
- puts seq_name
- }
- Util.print_message( PRG_NAME, "will process " + seqs.get_number_of_seqs.to_s + " non-redundant sequences" )
- end
-
- mapping_out = File.open( mapping_outfile, "a" )
-
- for i in 0 ... seqs.get_number_of_seqs
- seq = seqs.get_sequence( i )
- if seq.get_taxonomy == nil
- Util.fatal_error( PRG_NAME, "sequence [" + seq.get_name + "] has no taxonomy information" )
- end
- seq.set_name( Util::normalize_seq_name( modify_name( seq, i, sp_taxonomies, mapping_out ), 10 ) )
- end
-
- io = MsaIO.new()
-
- w = FastaWriter.new()
-
- w.set_max_name_length( 10 )
- w.clean( true )
- begin
- io.write_to_file( seqs, sequences_outfile, w )
- rescue Exception => e
- Util.fatal_error( PRG_NAME, "failed to write file: " + e.to_s )
- end
- mapping_out.close()
-
- Util.print_message( PRG_NAME, "wrote: " + mapping_outfile )
- Util.print_message( PRG_NAME, "wrote: " + sequences_outfile )
- Util.print_message( PRG_NAME, "OK" )
-
- end
-
- private
-
- def modify_name( seq, i, sp_taxonomies, mapping_outfile )
-
- tax_id = seq.get_taxonomy.get_id
- matching_sp_taxonomy = nil
-
- if @tax_ids_to_sp_taxonomies.has_key?( tax_id )
- # This is so that a second lookup will be much faster.
- matching_sp_taxonomy = @tax_ids_to_sp_taxonomies[ tax_id ]
- else
- sp_taxonomies.each { |sp_taxonomy|
- if ( sp_taxonomy.id == tax_id )
- if matching_sp_taxonomy != nil
- Util.fatal_error( PRG_NAME, "taxonomy id [" + tax_id.to_s + "] is not unique" )
- end
- matching_sp_taxonomy = sp_taxonomy
- @tax_ids_to_sp_taxonomies[ tax_id ] = sp_taxonomy
- end
- }
- end
- if matching_sp_taxonomy == nil
- Util.fatal_error( PRG_NAME, "taxonomy id [" + tax_id.to_s + "] for [" + seq.get_taxonomy.get_name + "] not found" )
- end
-
- new_name = i.to_s( 16 ) + "_" + matching_sp_taxonomy.code
-
- seq_name = seq.get_name
- if seq_name =~ /\[.+\]$/
- # Redundant taxonomy information hides here.
- seq_name = seq_name.sub(/\[.+\]$/, '')
- end
- if seq_name =~ /^\s*hypothetical\s+protein\s*/i
- # Pointless information.
- seq_name = seq_name.sub( /^\s*hypothetical\s+protein\s*/i, '' )
- end
-
- mapping_outfile.print( new_name + "\t" +
- TAXONOMY_CODE + matching_sp_taxonomy.code + "\t" +
- TAXONOMY_ID + tax_id + "\t" +
- TAXONOMY_ID_TYPE + seq.get_taxonomy.get_id_source + "\t" +
- TAXONOMY_SN + matching_sp_taxonomy.scientific_name + "\t" +
- SEQ_ACCESSION + seq.get_accession + "\t" +
- SEQ_ACCESSION_SOURCE + seq.get_accession_source + "\t" +
- SEQ_SYMBOL + seq.get_symbol + "\t" +
- SEQ_NAME + seq_name + "\t" +
- SEQ_MOL_SEQ + seq.get_sequence_as_string +
- Constants::LINE_DELIMITER )
- new_name
- end
-
- end
-
-end # module Evoruby
\ No newline at end of file
git://source.jalview.org / jalview.git / commitdiff