+++ /dev/null
-# $Id: forester.pm,v 1.26 2010/12/13 19:00:22 cmzmasek Exp $
-#
-# FORESTER -- software libraries and applications
-# for evolutionary biology research and applications.
-#
-# Copyright (C) 2007-2009 Christian M. Zmasek
-# Copyright (C) 2007-2009 Burnham Institute for Medical Research
-# All rights reserved
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License, or (at your option) any later version.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
-#
-# Contact: phylosoft @ gmail . com
-# WWW: www.phylosoft.org/forester
-#
-#
-#
-
-
-package forester;
-use strict;
-require Exporter;
-
-our $VERSION = 1.000;
-
-our @ISA = qw( Exporter );
-
-our @EXPORT = qw( executeConsense
- executePhyloPl
- executePuzzleDQO
- executePuzzleDQObootstrapped
- pfam2phylipMatchOnly
- startsWithSWISS_PROTname
- isPfamSequenceLine
- isPfamCommentLine
- containsPfamNamedSequence
- isRFline
- executeProtpars
- setModelForPuzzle
- setRateHeterogeneityOptionForPuzzle
- setParameterEstimatesOptionForPuzzle
- executePuzzleBootstrapped
- executePuzzle
- executeFastme
- executeNeighbor
- executeFitch
- executeBionj
- executeWeighbor
- executePhyml
- executeHmmfetch
- addDistsToQueryToPWDfile
- testForTextFilePresence
- exitWithWarning
- dieWithUnexpectedError
- addSlashAtEndIfNotPresent
- $LENGTH_OF_NAME
- $MIN_NUMBER_OF_AA
- $TREMBL_ACDEOS_FILE
- $SWISSPROT_ACDEOS_FILE
- $SPECIES_NAMES_FILE
- $SPECIES_TREE_FILE_DEFAULT
- $MULTIPLE_TREES_FILE_SUFFIX
- $LOG_FILE_SUFFIX
- $ALIGN_FILE_SUFFIX
- $TREE_FILE_SUFFIX
- $ADDITION_FOR_RIO_ANNOT_TREE
- $SUFFIX_PWD
- $SUFFIX_BOOT_STRP_POS
- $MULTIPLE_PWD_FILE_SUFFIX
- $SUFFIX_PWD_NOT_BOOTS
- $SUFFIX_HMM
- $MATRIX_FOR_PWD
- $RIO_PWD_DIRECTORY
- $RIO_BSP_DIRECTORY
- $RIO_NBD_DIRECTORY
- $RIO_ALN_DIRECTORY
- $RIO_HMM_DIRECTORY
- $PFAM_FULL_DIRECTORY
- $PFAM_SEED_DIRECTORY
- $PRIOR_FILE_DIR
- $PFAM_HMM_DB
- $FORESTER_JAR
- $SEQBOOT
- $NEIGHBOR
- $PROTPARS
- $CONSENSE
- $PROML
- $PHYLIP_VERSION
- $PUZZLE
- $PUZZLE_VERSION
- $FASTME
- $FASTME_VERSION
- $BIONJ
- $BIONJ_VERSION
- $WEIGHBOR
- $WEIGHBOR_VERSION
- $RAXML
- $RAXML_VERSION
- $PHYML
- $PHYML_VERSION
- $HMMALIGN
- $HMMSEARCH
- $HMMBUILD
- $HMMFETCH
- $SFE
- $HMMCALIBRATE
- $P7EXTRACT
- $MULTIFETCH
- $BOOTSTRAP_CZ
- $BOOTSTRAP_CZ_PL
- $SUPPORT_TRANSFER
- $SUPPORT_STATISTICS
- $NEWICK_TO_PHYLOXML
- $PHYLO_PL
- $RIO_PL
- $DORIO
- $PUZZLE_DQO
- $BOOTSTRAPS
- $PATH_TO_FORESTER
- $JAVA
- $NODE_LIST
- $RIO_SLAVE_DRIVER
- $RIO_SLAVE
- $TEMP_DIR_DEFAULT
- $EXPASY_SPROT_SEARCH_DE
- $EXPASY_SPROT_SEARCH_AC
- );
-
-
-
-
-# =============================================================================
-# =============================================================================
-#
-# THESE VARIABLES ARE ENVIRONMENT DEPENDENT, AND NEED TO BE SET ACCORDINGLY
-# BY THE USER
-# -------------------------------------------------------------------------
-#
-
-# For using just "phylo_pl.pl", only the following variables need to be set
-# $JAVA
-# $FORESTER_JAR
-# $TEMP_DIR_DEFAULT
-# $SEQBOOT
-# $CONSENSE
-# $PUZZLE
-# $FASTME
-# $NEIGHBOR
-# $FITCH
-# $BIONJ
-# $WEIGHBOR
-# $PHYML
-# $PROTPARS
-
-# Software directory:
-# ---------------------
-
-our $SOFTWARE_DIR = "/home/czmasek/SOFTWARE/";
-
-
-# Java virtual machine:
-# ---------------------
-our $JAVA = $SOFTWARE_DIR."JAVA/jdk1.6.0_03/bin/java";
-
-
-# Where all the temporary files can be created:
-# ---------------------------------------------
-our $TEMP_DIR_DEFAULT = "/tmp/";
-
-
-# Programs from Joe Felsenstein's PHYLIP package:
-# -----------------------------------------------
-our $SEQBOOT = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/seqboot";
-our $NEIGHBOR = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/neighbor";
-our $PROTPARS = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/protpars";
-our $PROML = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/proml";
-our $FITCH = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/fitch";
-our $CONSENSE = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/consense";
-our $PHYLIP_VERSION = "3.68";
-
-# TREE-PUZZLE:
-# ------------
-our $PUZZLE = $SOFTWARE_DIR."TREE_PUZZLE/tree-puzzle-5.2/src/puzzle";
-our $PUZZLE_VERSION = "5.2";
-
-# FASTME:
-# -----------------------------------------------------
-our $FASTME = $SOFTWARE_DIR."FASTME/fastme2.0/fastme";
-our $FASTME_VERSION = "2.0";
-
-# BIONJ:
-# -----------------------------------------------------
-our $BIONJ = $SOFTWARE_DIR."BIONJ/bionj";
-our $BIONJ_VERSION = "[1997]";
-
-# WEIGHBOR:
-# -----------------------------------------------------
-our $WEIGHBOR = $SOFTWARE_DIR."WEIGHBOR/Weighbor/weighbor";
-our $WEIGHBOR_VERSION = "1.2.1";
-
-# PHYML:
-# -----------------------------------------------------
-our $PHYML = $SOFTWARE_DIR."PHYML/phyml_v2.4.4/exe/phyml_linux";
-our $PHYML_VERSION = "2.4.4";
-
-# RAXML:
-# -----------------------------------------------------
-our $RAXML = $SOFTWARE_DIR."RAXML/RAxML-7.0.4/raxmlHPC";
-our $RAXML_VERSION = "7.0.4";
-
-
-# forester.jar. This jar file is currently available at: http://www.phylosoft.org
-# -------------------------------------------------------------------------------
-
-our $FORESTER_JAR = $SOFTWARE_DIR."FORESTER/DEV/forester/forester/java/forester.jar";
-
-
-
-# End of variables which need to be set by the user for using "phylo_pl.pl".
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-# Tool from forester.jar to transfer support values:
-# -------------------------------------------------
-our $SUPPORT_TRANSFER = $JAVA." -cp $FORESTER_JAR org.forester.application.support_transfer";
-
-
-
-# Tool from forester.jar for simple statistics for support values:
-# ----------------------------------------------------------------
-our $SUPPORT_STATISTICS = $JAVA." -cp $FORESTER_JAR org.forester.application.support_statistics";
-
-
-# Tool from forester.jar to transfer nh to phyloXML:
-# -------------------------------------------------
-our $NEWICK_TO_PHYLOXML = $JAVA." -cp $FORESTER_JAR org.forester.application.phyloxml_converter";
-
-
-
-# FORESTER itself (currently not needed for "phylo_pl.pl"):
-# ---------------------------------------------------------
-our $PATH_TO_FORESTER = "";
-
-
-# Pfam data (not needed for phylo_pl.pl):
-# --------------------------------------
-our $PFAM_FULL_DIRECTORY = "/path/to/Pfam/Full/";
-our $PFAM_SEED_DIRECTORY = "/path/to/Pfam/Seed/";
-our $PFAM_HMM_DB = "/path/to/Pfam/Pfam_ls"; # Need to run "hmmindex" on this
- # to produce .ssi file.
- # Then, for example
- # "setenv HMMERDB /home/rio/pfam-6.6/"
-
-
-$PATH_TO_FORESTER = &addSlashAtEndIfNotPresent( $PATH_TO_FORESTER );
-
-
-# Description lines and species from SWISS-PROT and TrEMBL (not needed for phylo_pl.pl):
-# -------------------------------------------------------------------------------------
-our $TREMBL_ACDEOS_FILE = $PATH_TO_FORESTER."data/trembl22_ACDEOS_1-6";
-
-our $SWISSPROT_ACDEOS_FILE = $PATH_TO_FORESTER."data/sp40_ACDEOS_1-6";
-
-
-
-# Names of species which can be analyzed and analyzed
-# against (must also be in tree $SPECIES_TREE_FILE_DEFAULT).
-# By using a list with less species, RIO analyses become faster
-# but lose phylogenetic resolution.
-# For many purposes, list "tree_of_life_bin_1-6_species_list"
-# in "data/species/" might be sufficient:
-# (not needed for phylo_pl.pl)
-# --------------------------------------------------------------
-our $SPECIES_NAMES_FILE = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6_species_list";
-
-
-
-# A default species tree in NHX format.
-# For many purposes, tree "tree_of_life_bin_1-6.nhx"
-# in "data/species/" might be fine:
-# (not needed for phylo_pl.pl)
-# --------------------------------------------------
-our $SPECIES_TREE_FILE_DEFAULT = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6.nhx";
-
-
-
-# Data for using precalculated distances:
-# (not needed for phylo_pl.pl)
-# ---------------------------------------
-our $MATRIX_FOR_PWD = 2; # The matrix which has been used for the pwd in $RIO_PWD_DIRECTORY.
- # 0=JTT, 1=PAM, 2=BLOSUM 62, 3=mtREV24, 5=VT, 6=WAG.
-
-our $RIO_PWD_DIRECTORY = $PATH_TO_FORESTER."example_data/"; # all must end with "/"
-our $RIO_BSP_DIRECTORY = $PATH_TO_FORESTER."example_data/";
-our $RIO_NBD_DIRECTORY = $PATH_TO_FORESTER."example_data/";
-our $RIO_ALN_DIRECTORY = $PATH_TO_FORESTER."example_data/";
-our $RIO_HMM_DIRECTORY = $PATH_TO_FORESTER."example_data/";
-
-
-
-#
-# End of variables which need to be set by the user.
-#
-# =============================================================================
-# =============================================================================
-
-
-
-
-
-$TEMP_DIR_DEFAULT = &addSlashAtEndIfNotPresent( $TEMP_DIR_DEFAULT );
-$PFAM_FULL_DIRECTORY = &addSlashAtEndIfNotPresent( $PFAM_FULL_DIRECTORY );
-$PFAM_SEED_DIRECTORY = &addSlashAtEndIfNotPresent( $PFAM_SEED_DIRECTORY );
-
-
-
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# These variables should normally not be changed:
-#
-
-our $PRIOR_FILE_DIR = $PATH_TO_FORESTER."data/priors_for_hmmbuild/";
- # Directory containing dirichlet prior
- # files needed for certain aligments
- # by hmmbuild (e.g. Collagen).
-
-
-
-
-
-# TREE-PUZZLE:
-our $PUZZLE_DQO = $PATH_TO_FORESTER."puzzle_dqo/src/puzzle";
-
-# HMMER:
-our $HMMALIGN = $PATH_TO_FORESTER."hmmer/binaries/hmmalign";
-our $HMMSEARCH = $PATH_TO_FORESTER."hmmer/binaries/hmmsearch";
-our $HMMBUILD = $PATH_TO_FORESTER."hmmer/binaries/hmmbuild";
-our $HMMFETCH = $PATH_TO_FORESTER."hmmer/binaries/hmmfetch";
-our $SFE = $PATH_TO_FORESTER."hmmer/binaries/sfetch";
-our $HMMCALIBRATE = $PATH_TO_FORESTER."hmmer/binaries/hmmcalibrate";
-
-our $P7EXTRACT = $PATH_TO_FORESTER."perl/p7extract.pl";
-our $MULTIFETCH = $PATH_TO_FORESTER."perl/multifetch.pl";
-
-
-# RIO/FORESTER:
-our $BOOTSTRAP_CZ = $PATH_TO_FORESTER."C/bootstrap_cz";
-our $BOOTSTRAP_CZ_PL = $PATH_TO_FORESTER."perl/bootstrap_cz.pl";
-#our $SUPPORT_TRANSFER = $JAVA." -cp $PATH_TO_FORESTER"."java forester.tools.transfersBranchLenghts";
-#our $SUPPORT_TRANSFER = $JAVA." -cp /home/czmasek/SOFTWARE/FORESTER/forester3/forester.jar org.forester.tools.SupportTransfer";
-
-our $PHYLO_PL = $PATH_TO_FORESTER."perl/phylo_pl.pl";
-our $RIO_PL = $PATH_TO_FORESTER."perl/rio.pl";
-our $DORIO = $JAVA." -cp $PATH_TO_FORESTER"."java forester.tools.DoRIO";
-# parallel RIO:
-our $RIO_SLAVE_DRIVER = $PATH_TO_FORESTER."perl/rio_slave_driver.pl";
-our $RIO_SLAVE = $PATH_TO_FORESTER."perl/rio_slave.pl";
-our $NODE_LIST = $PATH_TO_FORESTER."data/node_list.dat";
-
-#
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
-our $BOOTSTRAPS = 100;
-our $MIN_NUMBER_OF_AA = 20; # After removal of gaps, if less, gaps are not removed.
-our $LENGTH_OF_NAME = 10;
-
-
-
-
-our $MULTIPLE_TREES_FILE_SUFFIX = ".mlt";
-our $LOG_FILE_SUFFIX = ".log";
-our $ALIGN_FILE_SUFFIX = ".aln";
-our $TREE_FILE_SUFFIX = ".nhx";
-our $ADDITION_FOR_RIO_ANNOT_TREE = ".rio";
-our $SUFFIX_PWD = ".pwd";
-our $MULTIPLE_PWD_FILE_SUFFIX = ".mpwd";
-our $SUFFIX_BOOT_STRP_POS = ".bsp";
-our $SUFFIX_PWD_NOT_BOOTS = ".nbd";
-our $SUFFIX_HMM = ".hmm";
-
-our $EXPASY_SPROT_SEARCH_DE = "http://www.expasy.org/cgi-bin/sprot-search-de?";
-our $EXPASY_SPROT_SEARCH_AC = "http://www.expasy.org/cgi-bin/sprot-search-ac?";
-
-
-
-# One argument: input multiple trees file
-# Last modified: 07/05/01
-sub executeConsense {
- my $in = $_[ 0 ];
-
- &testForTextFilePresence( $in );
-
- system( "$CONSENSE >/dev/null 2>&1 << !
-$in
-Y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$CONSENSE $in\"" );
-
- return;
-}
-
-
-
-# Four arguments:
-# 1. options ("-" is not necessary)
-# 2. alignment or pwd file
-# 3. outfile
-# 4. temp dir
-# Last modified: 07/05/01
-sub executePhyloPl {
-
- my $opts = $_[ 0 ];
- my $B = $_[ 1 ];
- my $C = $_[ 2 ];
- my $D = $_[ 3 ];
-
- &testForTextFilePresence( $B );
-
- $opts = "-".$opts;
-
- system( "$PHYLO_PL $opts $B $C $D" )
- && &dieWithUnexpectedError( "Could not execute \"$PHYLO_PL $opts $B $C $D\"" );
-
-} ## executePhyloPl
-
-
-
-
-# Two arguments:
-# 1. Name of inputfile
-# 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
-# 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
-sub executePuzzleDQO {
- my $in = $_[ 0 ];
- my $matrix_option = $_[ 1 ];
- my $mat = "";
-
- &testForTextFilePresence( $in );
-
- $mat = setModelForPuzzle( $matrix_option );
-
- system( "$PUZZLE_DQO $in >/dev/null 2>&1 << !$mat
-y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$PUZZLE_DQO\"" );
-
- return;
-
-} ## executePuzzleDQO
-
-
-
-
-# Two arguments:
-# 1. Name of inputfile
-# 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
-# 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
-# Last modified: 01/28/02
-sub executePuzzleDQObootstrapped {
- my $in = $_[ 0 ];
- my $matrix_option = $_[ 1 ];
-
-
- my $l = 0;
- my $slen = 0;
- my $counter = 0;
- my $mat = "";
- my $a = "";
- my @a = ();
-
- &testForTextFilePresence( $in );
-
- open( GRP, "<$in" ) || &dieWithUnexpectedError( "Cannot open file \"$in\"" );
- while( <GRP> ) {
- if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
- $counter++;
- }
- }
- close( GRP );
-
- $l = `cat $in | wc -l`;
- $slen = $l / $counter;
-
- system( "split -$slen $in $in.splt." )
- && &dieWithUnexpectedError( "Could not execute \"split -$slen $in $in.splt.\"" );
-
- @a = <$in.splt.*>;
-
- $mat = setModelForPuzzle( $matrix_option );
-
- foreach $a ( @a ) {
-
- system( "$PUZZLE_DQO $a >/dev/null 2>&1 << !$mat
-y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$PUZZLE_DQO $a\"" );
-
- system( "cat $a.dist >> $in.dist" )
- && &dieWithUnexpectedError( "Could not execute \"cat outdist >> $in.dist\"" );
-
- unlink( $a, $a.".dist" );
- }
-
- return;
-
-} ## executePuzzleDQObootstrapped
-
-
-
-# Transfers a Pfam (SELEX) alignment to a
-# PHYLIP sequential style alignment.
-# It only writes "match columns" as indicated by the
-# "# RF" line ('x' means match).
-#
-# Three arguments:
-# 1. infile name
-# 2. outfile name
-# 3. 1 to NOT ensure that match states contain only 'A'-'Z' or '-'
-#
-# Returns the number of match states (=length of output alignment),
-# the length of the input alignment,
-# the number of seqs in the input alignment
-#
-# Last modified: 07/07/01
-#
-sub pfam2phylipMatchOnly {
-
- my $infile = $_[ 0 ];
- my $outfile = $_[ 1 ];
- my $ne = $_[ 2 ];
- my @seq_name = ();
- my @seq_array = ();
- my $return_line = "";
- my $seq = "";
- my $x = 0;
- my $y = 0;
- my $i = 0;
- my $x_offset = 0;
- my $max_x = 0;
- my $rf_y = 0;
- my $number_colum = 0;
- my $not_ensure = 0;
- my $saw_rf_line = 0;
-
- if ( $ne && $ne == 1 ) {
- $not_ensure = 1;
- }
-
- &testForTextFilePresence( $infile );
-
- open( INPP, "$infile" ) || &dieWithUnexpectedError( "Cannot open file \"$infile\"" );
-
- # This reads in the first block. It reads in the seq names.
- while ( 1 ) {
- if ( &isPfamSequenceLine( $return_line ) ) {
- $return_line =~ /^(\S+)\s+(\S+)/;
- $seq_name[ $y ] = substr( $1, 0, $LENGTH_OF_NAME );
- $seq = $2;
- for ( $x = 0; $x < length( $seq ); $x++ ) {
- $seq_array[ $x ][ $y ] = substr( $seq, $x, 1 );
- }
- $y++;
- }
- elsif ( &isRFline( $return_line ) ) {
- $saw_rf_line = 1;
- $return_line =~ /\s+(\S+)\s*$/;
- $seq = $1;
- $x_offset = length( $seq );
- $rf_y = $y;
- for ( $x = 0; $x < $x_offset; $x++ ) {
- $seq_array[ $x ][ $rf_y ] = substr( $seq, $x, 1 );
- }
- last;
- }
-
- $return_line = <INPP>;
-
- if ( !$return_line ) {
- &dieWithUnexpectedError( "Alignment not in expected format (no RF line)" );
- }
- }
-
- if ( $saw_rf_line != 1 ) {
- &dieWithUnexpectedError( "Alignment not in expected format (no RF line)" );
- }
-
- $y = 0;
- $max_x = 0;
-
- # This reads all blocks after the 1st one.
- while ( $return_line = <INPP> ) {
- if ( &isPfamSequenceLine( $return_line ) ) {
- $return_line =~ /^\S+\s+(\S+)/;
- $seq = $1;
- for ( $x = 0; $x < length( $seq ); $x++ ) {
- $seq_array[ $x + $x_offset ][ $y % $rf_y ] = substr( $seq, $x, 1 );
- }
- $y++;
- }
- elsif ( &isRFline( $return_line ) ) {
- if ( $y != $rf_y ) {
- &dieWithUnexpectedError( "Alignment not in expected format" );
- }
-
- $return_line =~ /\s+(\S+)\s*$/;
- $seq = $1;
- $max_x = length( $seq );
-
- for ( $x = 0; $x < length( $seq ); $x++ ) {
- $seq_array[ $x + $x_offset ][ $rf_y ] = substr( $seq, $x, 1 );
- }
-
- $y = 0;
- $x_offset = $x_offset + $max_x;
- $max_x = 0;
- }
- }
-
- close( INPP );
-
- # Counts the match states, and hence the number of aa in the alignment:
- for ( $x = 0; $x < $x_offset; $x++ ) {
- if ( !$seq_array[ $x ][ $rf_y ] ) {
- &dieWithUnexpectedError( "Alignment not in expected format" );
- }
- if ( $seq_array[ $x ][ $rf_y ] eq 'x' ) {
- $number_colum++;
- }
- }
-
- # Writes the file:
-
- open( OUTPP, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
- print OUTPP "$rf_y $number_colum\n";
- for ( $y = 0; $y < $rf_y; $y++ ) {
- print OUTPP "$seq_name[ $y ]";
- for ( $i = 0; $i < ( $LENGTH_OF_NAME - length( $seq_name[ $y ] ) ); $i++ ) {
- print OUTPP " ";
- }
- for ( $x = 0; $x < $x_offset; $x++ ) {
- if ( $seq_array[ $x ][ $rf_y ] eq 'x' ) {
- if ( !$seq_array[ $x ][ $y ] ) {
- &dieWithUnexpectedError( "Alignment not in expected format" );
- }
- if ( $not_ensure != 1 && $seq_array[ $x ][ $y ] !~ /[A-Z]|-/ ) {
- &dieWithUnexpectedError( "Alignment not in expected format (match states must only contain 'A'-'Z' or '-')" );
- }
- print OUTPP "$seq_array[ $x ][ $y ]";
- }
- }
- print OUTPP "\n";
- }
- close( OUTPP );
-
- return $number_colum, $x_offset, $rf_y;
-
-} ## pfam2phylipMatchOnly
-
-
-
-# Returns whether the argument (a String)
-# starts with a SWISS-PROT name (SEQN_SPECI).
-# Last modified: 06/21/01
-sub startsWithSWISS_PROTname {
- return ( $_[ 0 ] =~ /^[A-Z0-9]{1,4}_[A-Z0-9]{1,5}/ );
-}
-
-
-
-# Returns whether the argument starts with XXX.. XXXXX.. and the first
-# character is not a "#".
-# Last modified: 06/21/01
-sub isPfamSequenceLine {
- return( !&isPfamCommentLine( $_[ 0 ] )
- && &containsPfamNamedSequence( $_[ 0 ] ) );
-}
-
-
-
-# Returns whether the argument does start with a "#".
-# Last modified: 06/21/01
-sub isPfamCommentLine {
- return ( $_[ 0 ] =~ /^#/ );
-}
-
-
-
-# Returns whether the argument starts with XXX XXXXX.
-# Last modified: 06/21/01
-sub containsPfamNamedSequence {
- return ( $_[ 0 ] =~ /^\S+\s+\S+/ );
-}
-
-
-# Returns whether the argument starts with XXX XXXXX.
-# Last modified: 06/21/01
-sub isRFline {
- return ( $_[ 0 ] =~ /^#.*RF/ );
-}
-
-
-
-# Three arguments:
-# 1. pairwise distance file
-# 2. number of bootstraps
-# 3. initial tree: BME, GME or NJ
-# Last modified: 2008/12/31
-sub executeFastme {
- my $inpwd = $_[ 0 ];
- my $bs = $_[ 1 ];
- my $init_opt = $_[ 2 ];
-
- &testForTextFilePresence( $inpwd );
- my $command = "";
- if ( $bs > 0 ) {
- $command = "$FASTME -b $init_opt -i $inpwd -n $bs -s b";
- }
- else {
- $command = "$FASTME -b $init_opt -i $inpwd -s b";
- }
- print $command;
-
- system( $command );
-
-
-} ## executeFastme
-
-
-# Four arguments:
-# 1. pairwise distance file
-# 2. number of bootstraps
-# 3. seed for random number generator
-# 4. lower-triangular data matrix? 1: yes; no, otherwise
-sub executeNeighbor {
- my $inpwd = $_[ 0 ];
- my $bs = $_[ 1 ];
- my $s = $_[ 2 ];
- my $l = $_[ 3 ];
- my $multi = "";
- my $lower = "";
-
- &testForTextFilePresence( $inpwd );
-
- if ( $bs >= 2 ) {
- $multi = "
-M
-$bs
-$s";
- }
- if ( $l == 1 ) {
- $lower = "
-L";
- }
-
- system( "$NEIGHBOR >/dev/null 2>&1 << !
-$inpwd$multi$lower
-2
-3
-Y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$NEIGHBOR $inpwd$multi$lower\"" );
-
-} ## executeNeighbor
-
-
-# Seven arguments:
-# 1. pairwise distance file
-# 2. number of bootstraps
-# 3. seed for random number generator
-# 4. number of jumbles for input order
-# 5. lower-triangular data matrix? 1: yes; no, otherwise
-# 6. FM for Fitch-Margoliash, ME for ME# 6.
-# 7. 1 to use globale rearragements
-sub executeFitch {
- my $inpwd = $_[ 0 ];
- my $bs = $_[ 1 ];
- my $s = $_[ 2 ];
- my $j = $_[ 3 ];
- my $l = $_[ 4 ];
- my $m = $_[ 5 ];
- my $use_global_rearr = $_[ 6 ];
- my $jumble = "";
- my $multi = "";
- my $lower = "";
- my $method = "";
-
- my $global = "";
- if ( $use_global_rearr == 1 ) {
- $global = "
-G";
- }
-
- &testForTextFilePresence( $inpwd );
-
- if ( $m eq "FM" ) {
- $method = "";
- }
- elsif ( $m eq "ME" ) {
- $method = "
-D";
- }
- else {
- &dieWithUnexpectedError( "method for FITCH must be either FM or ME" );
- }
-
- if ( $j >= 1 ) {
- $jumble = "
-J
-$s
-$j";
- }
-
- if ( $bs >= 2 ) {
- $multi = "
-M
-$bs
-$s";
- }
- if ( $l == 1 ) {
- $lower = "
-L";
- }
-
- # jumble must be set BEFORE multi!
- system( "$FITCH 2>&1 << !
-$inpwd$method$global$jumble$multi$lower
-3
-Y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$FITCH $inpwd$method$global$jumble$multi$lower\"" );
- # 3: Do NOT print out tree
-
-} ## executeFitch
-
-
-
-# Two arguments:
-# 1. pairwise distance file
-# 2. outfile
-sub executeBionj {
- my $inpwd = $_[ 0 ];
- my $out = $_[ 1 ];
-
- &testForTextFilePresence( $inpwd );
- my $command = "$BIONJ $inpwd $out";
-
- system( $command )
- && &dieWithUnexpectedError( $command );
-
-}
-
-# Four arguments:
-# 1. (effective) sequence length
-# 2. (effective) number of bases
-# 3. pairwise distance file
-# 4. outfile
-sub executeWeighbor {
- my $L = $_[ 0 ];
- my $b = $_[ 1 ];
- my $i = $_[ 2 ];
- my $o = $_[ 3 ];
-
- &testForTextFilePresence( $i );
- my $command = "$WEIGHBOR -L $L -b $b -i $i -o $o";
-
- system( $command )
- && &dieWithUnexpectedError( $command );
-
-}
-
-# Six arguments:
-# 1. DNA or Amino-Acids sequence filename (PHYLIP format)
-# 2. number of data sets to analyse (ex:3)
-# 3. Model: JTT | MtREV | Dayhoff | WAG | VT | DCMut | Blosum62 (Amino-Acids)
-# 4. number of relative substitution rate categories (ex:4), positive integer
-# 5. starting tree filename (Newick format), your tree filename | BIONJ for a distance-based tree
-# 6. 1 to estimate proportion of invariable sites, otherwise, fixed proportion "0.0" is used
-# PHYML produces several results files :
-# <sequence file name>_phyml_lk.txt : likelihood value(s)
-# <sequence file name>_phyml_tree.txt : inferred tree(s)
-# <sequence file name>_phyml_stat.txt : detailed execution stats
-sub executePhyml {
- my $sequences = $_[ 0 ];
- my $data_sets = $_[ 1 ];
- my $model = $_[ 2 ];
- my $nb_categ = $_[ 3 ];
- my $tree = $_[ 4 ];
- my $estimate_invar_sites = $_[ 5 ];
-
- if ( $data_sets < 1 ) {
- $data_sets = 1
- }
-
- my $invar = "0.0"; # proportion of invariable sites,
- # a fixed value (ex:0.0) | e to get the maximum likelihood estimate
- if ( $estimate_invar_sites == 1 ) {
- $invar = "e";
- }
-
- my $data_type = "1"; # 0 = DNA | 1 = Amino-Acids
- my $format = "i"; # i = interleaved sequence format | s = sequential
- my $bootstrap_sets = "0"; # number of bootstrap data sets to generate (ex:2)
- # only works with one data set to analyse
-
- my $alpha = "e"; # gamma distribution parameter,
- # a fixed value (ex:1.0) | e to get the maximum likelihood estimate
-
- my $opt_topology = "y"; # optimise tree topology ? y | n
- my $opt_lengths = "y"; # optimise branch lengths and rate parameters ? y | n
-
- if ( $data_sets > 1 ) {
- # No need to calc branch lengths for bootstrapped analysis
- $opt_lengths = "n";
- }
-
- &testForTextFilePresence( $sequences );
- my $command = "$PHYML $sequences $data_type $format $data_sets $bootstrap_sets $model $invar $nb_categ $alpha $tree $opt_topology $opt_lengths";
-
- print( "\n$command\n");
-
- system( $command )
- && &dieWithUnexpectedError( $command );
-
-}
-
-
-
-
-# Four arguments:
-# 1. name of alignment file (in correct format!)
-# 2. number of bootstraps
-# 3. jumbles: 0: do not jumble; >=1 number of jumbles
-# 4. seed for random number generator
-sub executeProtpars {
- my $align = $_[ 0 ];
- my $bs = $_[ 1 ];
- my $rand = $_[ 2 ];
- my $s = $_[ 3 ];
- my $jumble = "";
- my $multi = "";
-
-
- &testForTextFilePresence( $align );
-
- if ( $bs > 1 && $rand < 1 ) {
- $rand = 1;
- }
-
- if ( $rand >= 1 ) {
- $jumble = "
-J
-$s
-$rand";
- }
-
- if ( $bs > 1 ) {
- $multi = "
-M
-D
-$bs";
- }
-
- system( "$PROTPARS 2>&1 << !
-$align$jumble$multi
-Y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$PROTPARS $align$jumble$multi\"" );
- # 3: Do NOT print out tree
-
-
- return;
-
-} ## executeProtpars
-
-
-
-# "Model of substitution" order for DQO TREE-PUZZLE 5.0:
-# Auto
-# m -> Dayhoff (Dayhoff et al. 1978)
-# m -> JTT (Jones et al. 1992)
-# m -> mtREV24 (Adachi-Hasegawa 1996)
-# m -> BLOSUM62 (Henikoff-Henikoff 92)
-# m -> VT (Mueller-Vingron 2000)
-# m -> WAG (Whelan-Goldman 2000)
-# m -> Auto
-# One argument:
-# matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
-# 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
-# Last modified: 07/07/01
-sub setModelForPuzzle {
- my $matrix_option = $_[ 0 ];
- my $matr = "";
-
- if ( $matrix_option == 0 ) { # JTT
- $matr = "
-m
-m";
- }
- elsif ( $matrix_option == 2 ) { # BLOSUM 62
- $matr = "
-m
-m
-m
-m";
- }
- elsif ( $matrix_option == 3 ) { # mtREV24
- $matr = "
-m
-m
-m";
- }
- elsif ( $matrix_option == 5 ) { # VT
- $matr = "
-m
-m
-m
-m
-m";
- }
- elsif ( $matrix_option == 6 ) { # WAG
- $matr = "
-m
-m
-m
-m
-m
-m";
- }
- elsif ( $matrix_option == 7 ) { # auto
- $matr = "";
- }
- else { # PAM
- $matr = "
-m"
- }
-
- return $matr;
-
-} ## setModelForPuzzle
-
-# One argument:
-# Model of rate heterogeneity:
-# 1 for "8 Gamma distributed rates"
-# 2 for "Two rates (1 invariable + 1 variable)"
-# 3 for "Mixed (1 invariable + 8 Gamma rates)"
-# otherwise: Uniform rate
-# Last modified: 09/08/03
-sub setRateHeterogeneityOptionForPuzzle {
- my $rate_heterogeneity_option = $_[ 0 ];
- my $opt = "";
-
- if ( $rate_heterogeneity_option == 1 ) {
- $opt = "
-w";
- }
- elsif ( $rate_heterogeneity_option == 2 ) {
- $opt = "
-w
-w";
- }
- elsif ( $rate_heterogeneity_option == 3 ) {
- $opt = "
-w
-w
-w";
- }
- else {
- $opt = "";
- }
-
- return $opt;
-} ## setRateHeterogeneityOptionForPuzzle
-
-
-# One argument:
-# Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
-# Last modified: 09/08/03
-sub setParameterEstimatesOptionForPuzzle {
- my $parameter_estimates_option = $_[ 0 ];
- my $opt = "";
-
- if ( $parameter_estimates_option == 1 ) {
- $opt = "
-e";
- }
- else {
- $opt = "";
- }
-
- return $opt;
-} ## setParameterEstimatesOptionForPuzzle
-
-
-
-# three/four/five arguments:
-# 1. Name of inputfile
-# 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
-# 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
-# 3. Number of sequences in alignment
-# 4. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
-# 5. Model of rate heterogeneity:
-# 1 for "8 Gamma distributed rates"
-# 2 for "Two rates (1 invariable + 1 variable)"
-# 3 for "Mixed (1 invariable + 8 Gamma rates)"
-# otherwise: Uniform rate
-sub executePuzzleBootstrapped {
- my $in = $_[ 0 ];
- my $matrix_option = $_[ 1 ];
- my $number_of_seqs = $_[ 2 ];
- my $parameter_estimates_option = $_[ 3 ];
- my $rate_heterogeneity_option = $_[ 4 ];
-
- my $l = 0;
- my $slen = 0;
- my $counter = 0;
- my $mat = "";
- my $est = "";
- my $rate = "";
- my $a = "";
- my @a = ();
-
- &testForTextFilePresence( $in );
-
- open( GRP, "<$in" ) || die "\n\n$0: Unexpected error: Cannot open file <<$in>>: $!";
- while( <GRP> ) {
- if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
- $counter++;
- }
- }
- close( GRP );
-
- $l = `cat $in | wc -l`;
- $slen = $l / $counter;
-
- system( "split --suffix-length=4 -$slen $in $in.splt." )
- && die "\n\n$0: executePuzzleDQObootstrapped: Could not execute \"split --suffix-length=4 -$slen $in $in.splt.\": $!";
-
- @a = <$in.splt.*>;
-
- $mat = setModelForPuzzle( $matrix_option );
- if ( $parameter_estimates_option ) {
- $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
- }
- if ( $rate_heterogeneity_option ) {
- $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
- }
-
- my $k="";
- if ( $number_of_seqs <= 257 ) {
- $k = "k";
- }
-
- foreach $a ( @a ) {
- print "-".$a."\n";
- system( "$PUZZLE $a << !
-$k
-k
-k$mat$est$rate
-y
-!" )
- && die "$0: Could not execute \"$PUZZLE $a\"";
-
- system( "cat $a.dist >> $in.dist" )
- && die "$0: Could not execute \"cat outdist >> $in.dist\"";
-
- unlink( $a, $a.".dist", $a.".puzzle" );
- }
-
- return;
-
-} ## executePuzzleBootstrapped
-
-
-
-
-
-# three/four/five arguments:
-# 1. Name of inputfile
-# 2. Matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
-# 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
-# 3. Number of sequences in alignment
-# 4. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
-# 5. Model of rate heterogeneity:
-# 1 for "8 Gamma distributed rates"
-# 2 for "Two rates (1 invariable + 1 variable)"
-# 3 for "Mixed (1 invariable + 8 Gamma rates)"
-# otherwise: Uniform rate
-sub executePuzzle {
- my $in = $_[ 0 ];
- my $matrix_option = $_[ 1 ];
- my $number_of_seqs = $_[ 2 ];
- my $parameter_estimates_option = $_[ 3 ];
- my $rate_heterogeneity_option = $_[ 4 ];
- my $mat = "";
- my $est = "";
- my $rate = "";
-
- &testForTextFilePresence( $in );
-
- $mat = &setModelForPuzzle( $matrix_option );
- if ( $parameter_estimates_option ) {
- $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
- }
- if ( $rate_heterogeneity_option ) {
- $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
- }
-
- my $k="";
- if ( $number_of_seqs <= 257 ) {
- $k = "k";
- }
-
-
- system( "$PUZZLE $in << !
-$k
-k
-k$mat$est$rate
-y
-!" )
- && die "$0: Could not execute \"$PUZZLE\"";
-
- return;
-
-} ## executePuzzle
-
-
-
-
-# Preparation of the pwd file
-sub addDistsToQueryToPWDfile {
- my $pwd_file = $_[ 0 ];
- my $disttoquery_file = $_[ 1 ];
- my $outfile = $_[ 2 ];
- my $name_of_query = $_[ 3 ];
- my $name_of_query_ = "";
- my $return_line_pwd = "";
- my $return_line_dq = "";
- my $num_of_sqs = 0;
- my $block = 0;
- my $name_from_pwd = "X";
- my $name_from_dq = "Y";
- my @dists_to_query = ();
- my $i = 0;
-
- &testForTextFilePresence( $pwd_file );
- &testForTextFilePresence( $disttoquery_file );
-
- $name_of_query_ = $name_of_query;
- for ( my $j = 0; $j <= ( $LENGTH_OF_NAME - length( $name_of_query ) - 1 ); ++$j ) {
- $name_of_query_ .= " ";
- }
-
- open( OUT_AD, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
- open( IN_PWD, "$pwd_file" ) || &dieWithUnexpectedError( "Cannot open file \"$pwd_file\"" );
- open( IN_DQ, "$disttoquery_file" ) || &dieWithUnexpectedError( "Cannot open file \"$disttoquery_file\"" );
-
- W: while ( $return_line_pwd = <IN_PWD> ) {
-
-
- if ( $return_line_pwd =~ /^\s*(\d+)\s*$/ ) {
- $num_of_sqs = $1;
- $num_of_sqs++;
- if ( $block > 0 ) {
- print OUT_AD "$name_of_query_ ";
- for ( my $j = 0; $j < $i; ++$j ) {
- print OUT_AD "$dists_to_query[ $j ] ";
- }
- print OUT_AD "0.0\n";
- }
- print OUT_AD " $num_of_sqs\n";
- $block++;
- @dists_to_query = ();
- $i = 0;
- }
-
- if ( $block == 1
- && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
- $name_from_pwd = $1;
-
- if ( !defined( $return_line_dq = <IN_DQ> ) ) {
- &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
- }
-
- if ( $return_line_dq !~ /\S/ ) {
- if ( !defined( $return_line_dq = <IN_DQ> ) ) {
- &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
- }
- }
- $return_line_dq =~ /^\s*(\S+)\s+(\S+)/;
- $name_from_dq = $1;
- $dists_to_query[ $i++ ] = $2;
-
-
- if ( $name_from_pwd ne $name_from_dq ) {
- &dieWithUnexpectedError( "Order of sequence names in \"$pwd_file\" and \"$disttoquery_file\" is not the same" );
- }
- print OUT_AD $return_line_pwd;
-
- }
- elsif ( $block > 1
- && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
- $name_from_pwd = $1;
- if ( !defined( $return_line_dq = <IN_DQ> ) ) {
- &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
- }
- if ( $return_line_dq !~ /\S/ ) {
- if ( !defined( $return_line_dq = <IN_DQ>) ) {
- &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
- }
- }
- $return_line_dq =~ /^\s*\S+\s+(\S+)/;
- $dists_to_query[ $i++ ] = $1;
- print OUT_AD $return_line_pwd;
- }
- }
- print OUT_AD "$name_of_query_ ";
- for ( my $j = 0; $j < $i; ++$j ) {
- print OUT_AD "$dists_to_query[ $j ] ";
- }
- print OUT_AD "0.0\n";
-
- close( OUT_AD );
- close( IN_PWD );
- close( IN_DQ );
- return $block;
-
-} ## addDistsToQueryToPWDfile
-
-
-
-
-# Three arguments:
-# 1. HMMER model db
-# 2. name of HMM
-# 3. outputfile name
-# Last modified: 02/27/01
-sub executeHmmfetch {
-
- my $db = $_[ 0 ];
- my $name = $_[ 1 ];
- my $outfile = $_[ 2 ];
-
- system( "$HMMFETCH $db $name > $outfile" )
- && &dieWithUnexpectedError( "Could not execute \"$HMMFETCH $db $name > $outfile\"" );
- return;
-
-} ## executeHmmfetch
-
-
-
-# Checks wether a file is present, not empty and a plain textfile.
-# One argument: name of file.
-# Last modified: 07/07/01
-sub testForTextFilePresence {
- my $file = $_[ 0 ];
- unless ( ( -s $file ) && ( -f $file ) && ( -T $file ) ) {
- dieWithUnexpectedError( "File \"$file\" does not exist, is empty, or is not a plain textfile" );
- }
-} ## testForTextFilePresence
-
-
-# Last modified: 02/21/03
-sub addSlashAtEndIfNotPresent {
- my $filename = $_[ 0 ];
- $filename =~ s/\s+//g;
- unless ( $filename =~ /\/$/ ) {
- $filename = $filename."/";
- }
- return $filename;
-} ## addSlashAtEndIfNotPresent
-
-
-
-# Last modified: 02/15/02
-sub exitWithWarning {
-
- my $text = $_[ 0 ];
- if ( defined( $_[ 1 ] ) && $_[ 1 ] == 1 ) {
- print( "<H4 class=\"error\">user error</H4>\n" );
- print( "<P>\n" );
- print( "<B>$text</B>\n" );
- print( "</P>\n" );
- print( "<P>  </P>\n" );
- }
- else {
- print( "\n\n$text\n\n" );
- }
-
- exit( 0 );
-
-} ## exit_with_warning
-
-
-
-# Last modified: 02/15/02
-sub dieWithUnexpectedError {
-
- my $text = $_[ 0 ];
-
- die( "\n\n$0:\nUnexpected error (should not have happened):\n$text\n$!\n\n" );
-
-} ## dieWithUnexpectedError
-
-
-
-1;
+++ /dev/null
-#!/usr/bin/perl -W
-#
-# $Id: phylo_pl.pl,v 1.32 2010/12/13 19:00:22 cmzmasek Exp $
-#
-# FORESTER -- software libraries and applications
-# for evolutionary biology research and applications.
-#
-# Copyright (C) 2008-2009 Christian M. Zmasek
-# Copyright (C) 2008-2009 Burnham Institute for Medical Research
-# All rights reserved
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License, or (at your option) any later version.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
-#
-# Contact: phylosoft @ gmail . com
-# WWW: www.phylosoft.org/forester
-#
-#
-#
-# Requirements phylo_pl is part of the FORESTER libraries.
-# ------------ Many of its global variables are set via forester.pm.
-#
-#
-# Note. Use xt.pl (for Pfam alignments) or mt.pl (for other alignments)
-# to run phylo_pl.pl on whole directories of alignments files.
-#
-#
-#
-#
-# =========================
-#
-# METHOD ORDER (IMPORTANT!)
-# 1. FastME
-# 2. phylip NJ
-# 3. phylip fitch FM
-# 4. phylip fitch ME
-# 5. BIONJ
-# 6. Weighbor
-# 7. Raxml
-# 8. phyml
-# 9. phylip proml
-# 10. phylip protpars
-# 11. all
-#==========================
-
-use strict;
-
-use FindBin;
-use lib $FindBin::Bin;
-use forester;
-
-my $VERSION = "1.0.1";
-my $LAST_MODIFIED = "2009.10.02";
-
-my $RAXML_MODEL_BASE = "PROTGAMMA";
-my $RAXML_ALGORITHM = "a";
-
-my $TEMP_DIR_DEFAULT = "/tmp/phylo_pl_"; # Where all the infiles, outfiles, etc will be created.
-
-my $bootstraps = 100; # 0,1: do not bootstrap. Default: 100
-my $matrix = 5; # 0 = JTT
- # 1 = PAM
- # 2 = BLOSUM 62
- # 3 = mtREV24
- # 5 = VT - default
- # 6 = WAG
- # 7 = auto in puzzle
- # 8 = DCMut in PHYML, VT in TREE-PUZZLE
-my $rate_heterogeneity = 0; # 0 = Uniform rate (default)
- # 1 = 8 Gamma distributed rates
- # 2 = Two rates (1 invariable + 1 variable)
- # 3 = Mixed (1 invariable + 8 Gamma rates)
-my $seed = 9; # Seed for random number generators. Default: 9
-my $keep_multiple_trees = 0; # 0: delete multiple tree file
- # 1: do not delete multiple tree file
-my $exact_parameter_est = 0; # 0: no; 1: yes
-
-my $phyml_rel_substitution_rate_cat = 4;
-
-my $jumbles = 2;
-my $use_fastme = 0; # 0: no; 1: yes
-my $use_phylip_nj = 0; # 0: no; 1: yes
-my $use_phylip_fitch_fm = 0; # 0: no; 1: yes
-my $use_phylip_fitch_me = 0; # 0: no; 1: yes
-my $use_bionj = 0; # 0: no; 1: yes
-my $use_weighbor = 0; # 0: no; 1: yes
-my $use_raxml = 0; # 0: no; 1: yes
-my $use_phyml = 0; # 0: no; 1: yes
-my $use_proml = 0; # 0: no; 1: yes
-my $use_protpars = 0; # 0: no; 1: yes
-my $use_global_rearr = 0; # 0: no; 1: yes
-my $estimate_invar_sites = 0; # 0: no; 1: yes
-
-my $fastme_init_tree_opt = "NJ";
-
-my %seqnames = (); # number => seqname
-my %numbers = (); # seqname => number
-my $options = "";
-my $infile = "";
-my $pwdfile = "";
-my $outfile = "";
-my $logfile = "";
-my $multipwdfile = "";
-my $distancefile = "";
-my $log = "";
-my $ii = 0;
-my $temp_dir = "";
-my $current_dir = "";
-my @out = ();
-my $number_of_seqs = 0;
-my $number_of_aa = 0;
-
-my $use_pwd_based_methods = 0;
-
-print( "\n");
-print( "phylo_pl $VERSION ($LAST_MODIFIED)\n" );
-print( "__________________________________\n");
-print( "\n\n");
-
-
-
-unless ( @ARGV == 2 || @ARGV == 3 || @ARGV == 4 || @ARGV == 5 ) {
- &printUsage();
- exit ( -1 );
-}
-
-
-
-# Analyzes the options:
-# ---------------------
-
-if ( $ARGV[ 0 ] =~ /^-.+/ ) {
-
- unless ( @ARGV > 2 ) {
- &printUsage();
- exit ( -1 );
- }
- $options = $ARGV[ 0 ];
-
- if ( @ARGV != 3 && @ARGV != 4 ) {
- &printUsage();
- exit ( -1 );
- }
- $infile = $ARGV[ 1 ];
- $outfile = $ARGV[ 2 ];
- if ( @ARGV == 4 ) {
- $temp_dir = $ARGV[ 3 ];
- }
- if ( $options =~ /B(\d+)/ ) {
- $bootstraps = $1;
- if ( $bootstraps <= 1 ) {
- $bootstraps = 0;
- }
- elsif ( $bootstraps <= 9 ) {
- $bootstraps = 0;
- print "\n\nphylo_pl: WARNING: Bootstrap number must be devisable by 10,\nno bootstrapping.\n\n";
- }
- elsif ( $bootstraps % 10 != 0 ) {
- $bootstraps = $bootstraps - $bootstraps % 10; # to ensure $bootstraps % 10 == 0
- print "\n\nphylo_pl: WARNING: Bootstrap number must be devisable by 10,\nset to $bootstraps.\n\n";
- }
- }
- if ( $options =~ /n/ ) {
- $use_phylip_nj = 1;
- }
- if ( $options =~ /q@(\d)/ ) {
- $use_fastme = 1;
- my $opt = $1;
- if ( $opt == 1 ) {
- $fastme_init_tree_opt = "GME";
- }
- elsif ( $opt == 2 ) {
- $fastme_init_tree_opt = "BME";
- }
- elsif ( $opt == 3 ) {
- $fastme_init_tree_opt = "NJ";
- }
- else {
- &printUsage();
- exit ( -1 );
- }
- }
- if ( $options =~ /f/ ) {
- $use_phylip_fitch_fm = 1;
- }
- if ( $options =~ /e/ ) {
- $use_phylip_fitch_me = 1;
- }
- if ( $options =~ /b/ ) {
- $use_bionj = 1;
- }
- if ( $options =~ /w/ ) {
- $use_weighbor = 1;
- }
- if ( $options =~ /x/ ) {
- $use_raxml = 1;
- }
- if ( $options =~ /y/ ) {
- $use_phyml = 1;
- }
- if ( $options =~ /o/ ) {
- $use_proml = 1;
- }
- if ( $options =~ /p/ ) {
- $use_protpars = 1;
- }
- if ( $options =~ /G/ ) {
- $use_global_rearr = 1;
- }
- if ( $options =~ /I/ ) {
- $estimate_invar_sites = 1;
- }
- if ( $options =~ /j(\d+)/ ) {
- $jumbles = $1;
- if ( $jumbles < 1 ) {
- $jumbles = 0;
- }
- }
- if ( $options =~ /r(\d+)/ ) {
- $phyml_rel_substitution_rate_cat = $1;
- if ( $phyml_rel_substitution_rate_cat < 1 ) {
- $phyml_rel_substitution_rate_cat = 1;
- }
- }
- if ( $options =~ /J/ ) {
- $matrix = 0; # JTT
- }
- if ( $options =~ /P/ ) {
- $matrix = 1; # PAM
- }
- if ( $options =~ /L/ ) {
- $matrix = 2; # Blosum 62
- }
- if ( $options =~ /M/ ) {
- $matrix = 3; # mtREV24
- }
- if ( $options =~ /W/ ) {
- $matrix = 6; # WAG
- }
- if ( $options =~ /A/ ) {
- $matrix = 7; # auto
- }
- if ( $options =~ /D/ ) {
- $matrix = 8; # DCMut in PHYML and RAXML, VT in PUZZLE
- }
- if ( $options =~ /S(\d+)/ ) {
- $seed = $1;
- }
- if ( $options =~ /X/ ) {
- $keep_multiple_trees = 1;
- }
- if ( $options =~ /E/ ) {
- $exact_parameter_est = 1;
- }
- if ( $options =~ /g/ ) {
- $rate_heterogeneity = 1;
- }
- if ( $options =~ /t/ ) {
- $rate_heterogeneity = 2;
- }
- if ( $options =~ /m/ ) {
- $rate_heterogeneity = 3;
- }
-}
-else {
- unless ( @ARGV == 2 || @ARGV == 3 ) {
- &printUsage();
- exit ( -1 );
- }
- $infile = $ARGV[ 0 ];
- $outfile = $ARGV[ 1 ];
- if ( @ARGV == 3 ) {
- $temp_dir = $ARGV[ 2 ];
- }
-}
-
-if ( $use_fastme != 1 &&
- $use_phylip_nj != 1 &&
- $use_phylip_fitch_fm != 1 &&
- $use_phylip_fitch_me != 1 &&
- $use_bionj != 1 &&
- $use_weighbor != 1 &&
- $use_raxml != 1 &&
- $use_phyml != 1 &&
- $use_proml != 1 &&
- $use_protpars != 1 ) {
-
- $use_fastme = 1;
- $use_phylip_nj = 1;
- $use_phylip_fitch_fm = 1;
- $use_phylip_fitch_me = 1;
- $use_bionj = 1;
- $use_raxml = 1;
- $use_weighbor = 1;
- $use_phyml = 1;
- $use_proml = 1;
- $use_protpars = 1;
-}
-
-
-if ( $use_fastme == 1 ||
- $use_phylip_nj == 1 ||
- $use_phylip_fitch_fm == 1 ||
- $use_phylip_fitch_me == 1 ||
- $use_bionj == 1 ||
- $use_weighbor == 1 ) {
- $use_pwd_based_methods = 1;
-}
-else {
- $use_pwd_based_methods = 0;
-}
-
-$current_dir = `pwd`;
-$current_dir =~ s/\s//;
-
-if ( $outfile !~ /^\// ) {
- # outfile is not absolute path.
- $outfile = $current_dir."/".$outfile;
-}
-
-
-
-# TREE-PUZZLE sets the option in this way:
-# If two rates or mixed, exact parameter estimates are used.
-if ( $rate_heterogeneity == 2
-|| $rate_heterogeneity == 3 ) {
- $exact_parameter_est = 1
-}
-
-
-if ( $outfile =~ /\.xml$/i ) {
- $outfile =~ s/\.xml//i;
-}
-elsif ( $outfile =~ /\.aln$/i ) {
- $outfile =~ s/\.aln//i;
-}
-elsif ( $outfile =~ /\.fasta$/i ) {
- $outfile =~ s/\.fasta//i;
-}
-elsif ( $outfile =~ /\.fas$/i ) {
- $outfile =~ s/\.fas//i;
-}
-elsif ( $outfile =~ /\.seqs$/i ) {
- $outfile =~ s/\.seqs//i;
-}
-
-
-$logfile = $outfile.$LOG_FILE_SUFFIX;
-$multipwdfile = $outfile.$MULTIPLE_PWD_FILE_SUFFIX;
-$distancefile = $outfile.$SUFFIX_PWD_NOT_BOOTS;
-
-&dieIfFileExists( $logfile );
-&dieIfFileExists( $multipwdfile );
-&dieIfFileExists( $distancefile );
-
-
-my $fastme_outtree = $outfile."_fme.xml";
-my $phylip_nj_outtree = $outfile."_pnj.xml";
-my $phylip_fm_outtree = $outfile."_pfm.xml";
-my $phylip_me_outtree = $outfile."_pme.xml";
-my $bionj_outtree = $outfile."_bionj.xml";
-my $weighbor_outtree = $outfile."_weigh.xml";
-my $raxml_outtree = $outfile."_raxml.xml";
-my $phyml_outtree = $outfile."_phyml.xml";
-my $proml_outtree = $outfile."_proml.xml";
-my $protpars_outtree = $outfile."_ppp.xml";
-my $all_outtree = $outfile."_comb.xml";
-
-my $multitreefile_fastme = $outfile."_fme".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_phylip_nj = $outfile."_pnj".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_phylip_fm = $outfile."_pfm".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_phylip_me = $outfile."_pme".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_bionj = $outfile."_bionj".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_weighbor = $outfile."_weigh".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_raxml = $outfile."_raxml".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_phyml = $outfile."_phyml".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_proml = $outfile."_proml".$MULTIPLE_TREES_FILE_SUFFIX;
-my $multitreefile_protpars = $outfile."_ppp".$MULTIPLE_TREES_FILE_SUFFIX;
-
-if ( $use_fastme == 1 ) {
- &dieIfFileExists( $fastme_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_fastme );
- }
-}
-if( $use_phylip_nj == 1 ) {
- &dieIfFileExists( $phylip_nj_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_phylip_nj );
- }
-}
-if( $use_phylip_fitch_fm == 1 ) {
- &dieIfFileExists( $phylip_fm_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_phylip_fm );
- }
-}
-if( $use_phylip_fitch_me == 1 ) {
- &dieIfFileExists( $phylip_me_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_phylip_me );
- }
-}
-if( $use_bionj == 1 ) {
- &dieIfFileExists( $bionj_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_bionj );
- }
-}
-if( $use_weighbor == 1 ) {
- &dieIfFileExists( $weighbor_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_weighbor );
- }
-}
-if( $use_raxml == 1 ) {
- &dieIfFileExists( $raxml_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_raxml );
- }
-}
-if( $use_phyml == 1 ) {
- &dieIfFileExists( $phyml_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_phyml );
- }
-}
-if( $use_proml == 1 ) {
- &dieIfFileExists( $proml_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_proml );
- }
-}
-if( $use_protpars == 1 ) {
- &dieIfFileExists( $protpars_outtree );
- if ( $keep_multiple_trees == 1 && $bootstraps > 1 ) {
- &dieIfFileExists( $multitreefile_protpars );
- }
-}
-if ( $bootstraps > 1 ) {
- &dieIfFileExists( $all_outtree );
-}
-if ( $infile ne "" ) {
- unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
- die "\n\nphylo_pl: Input alignment file \"$infile\" does not exist, is empty, or is not a plain textfile.\n\n";
- }
-}
-
-
-
-
-# Prints out the options:
-# -----------------------
-
-
-$log = "\n$0 logfile:\n";
-$log = $log."Version: $VERSION\n\n";
-
-
-
-if ( $infile ne "" ) {
- $log = $log."Input : $infile\n";
-}
-
-if ( $keep_multiple_trees == 1 && $bootstraps >= 2 ) {
- $log = $log."Multiple distance matrices : $multipwdfile\n";
-}
-
-
-$log = $log."Bootstraps : $bootstraps\n";
-
-if ( $use_pwd_based_methods == 1 ) {
- $log = $log."Prgrm to calculate pairwise dist. : TREE-PUZZLE (version: $PUZZLE_VERSION)\n";
-}
-
-
-if ( $use_fastme == 1 ) {
- $log = $log."Program to calculate tree : FastME (version: $FASTME_VERSION)\n";
- $log = $log."Method for intial tree in FastME : $fastme_init_tree_opt\n";
- $log = $log."Tree swapping (NNI) in FastME : balanced (default)\n";
-}
-if ( $use_phylip_nj == 1 ) {
- $log = $log."Program to calculate tree : PHYLIP NEIGHBOR NJ (version: $PHYLIP_VERSION)\n";
-}
-if ( $use_phylip_fitch_fm == 1 ) {
- $log = $log."Program to calculate tree : PHYLIP FITCH Fitch-Margoliash (version: $PHYLIP_VERSION)\n";
-}
-if ( $use_phylip_fitch_me == 1 ) {
- $log = $log."Program to calculate tree : PHYLIP FITCH Minimal Evolution (version: $PHYLIP_VERSION)\n";
-}
-if ( $use_bionj == 1 ) {
- $log = $log."Program to calculate tree : BIONJ (version: $BIONJ_VERSION)\n";
-}
-if ( $use_weighbor == 1 ) {
- $log = $log."Program to calculate tree : Weighbor [no invariable sites, b=14] (version: $WEIGHBOR_VERSION)\n";
-}
-if ( $use_raxml == 1 ) {
- $log = $log."Program to calculate tree : RAxML [$RAXML_MODEL_BASE] (uses its own bootstraps, if bootstrapped: -f $RAXML_ALGORITHM) (version: $RAXML_VERSION)\n";
-}
-if ( $use_phyml == 1 ) {
- $log = $log."Program to calculate tree : PHYML (MLE for gamma distr param and proportion of inv sites) (version: $PHYML_VERSION)\n";
- $log = $log."# of rel subst rate categories : $phyml_rel_substitution_rate_cat\n";
-}
-if ( $use_proml == 1 ) {
- $log = $log."Program to calculate tree : PHYLIP PROML (uses PAM unless JTT selected) (version: $PHYLIP_VERSION)\n";
-}
-if ( $use_protpars == 1 ) {
- $log = $log."Program to calculate tree : PHYLIP PROTPARS (with global rearrangements) (version: $PHYLIP_VERSION)\n";
-}
-if ( $use_phylip_fitch_fm == 1 || $use_phylip_fitch_me == 1 || $use_protpars == 1 || $use_proml ) {
- $log = $log."Number of jumbles (input order rand): $jumbles\n";
-
-}
-if ( $use_phylip_fitch_fm == 1 || $use_phylip_fitch_me == 1 || $use_proml ) {
- if ( $use_global_rearr == 1 ) {
- $log = $log."Global rearrangements : true\n";
- }
- else {
- $log = $log."Global rearrangements : false\n";
-
- }
-}
-
-if ( $bootstraps > 0 ) {
- $log = $log."Prgrm to calculate ML branch lenghts: TREE-PUZZLE (version: $PUZZLE_VERSION)\n";
-}
-
-$log = $log."Model : ";
-if ( $matrix == 0 ) {
- $log = $log."JTT (Jones et al. 1992)\n";
-}
-elsif ( $matrix == 1 ) {
- $log = $log."PAM (Dayhoff et al. 1978)\n";
-}
-elsif ( $matrix == 2 ) {
- $log = $log."BLOSUM 62 (Henikoff-Henikoff 92)\n";
-}
-elsif ( $matrix == 3 ) {
- $log = $log."mtREV24 (Adachi-Hasegawa 1996)\n";
-}
-elsif ( $matrix == 5 ) {
- $log = $log."VT (Mueller-Vingron 2000)\n";
-}
-elsif ( $matrix == 6 ) {
- $log = $log."WAG (Whelan-Goldman 2000)\n";
-}
-elsif ( $matrix == 7 ) {
- $log = $log."auto in TREE-PUZZLE\n";
-}
-elsif ( $matrix == 8 ) {
- $log = $log."DCMut (Kosial and Goldman, 2005) in PHYML and RAxML, VT in TREE-PUZZLE\n";
-}
-else {
- &dieWithUnexpectedError( "Unknown model: matrix=$matrix" );
-}
-if ( $use_raxml == 1 || $use_phyml == 1 ) {
- if ( $estimate_invar_sites == 1 ) {
- $log = $log."Estimate proportion of invariable sites in RAXML and/or PHYML: true\n";
- }
- else {
- $log = $log."Estimate proportion of invariable sites in RAXML and/or PHYML: false (proportion \"0.0\" is used in PHYML)\n";
- }
-}
-
-$log = $log."Model of rate heterogeneity (PUZZLE): ";
-if ( $rate_heterogeneity == 1 ) {
- $log = $log."8 Gamma distributed rates\n";
-}
-elsif ( $rate_heterogeneity == 2 ) {
- $log = $log."Two rates (1 invariable + 1 variable)\n";
-}
-elsif ( $rate_heterogeneity == 3 ) {
- $log = $log."Mixed (1 invariable + 8 Gamma rates)\n";
-}
-else {
- $log = $log."Uniform rate\n";
-}
-$log = $log."Seed for random number generators : $seed\n";
-if ( $exact_parameter_est == 1 ) {
- $log = $log."Exact parameter estimates in TREE-PUZZLE\n";
-}
-
-$log = $log."Start time/date : ".`date`;
-
-
-
-
-# That's where the mischief starts....
-# ------------------------------------
-
-$ii = 0;
-
-srand();
-my $range = 1000000;
-my $random_number = int( rand( $range ) );
-
-if ( $temp_dir eq "" ) {
- $temp_dir = $TEMP_DIR_DEFAULT;
-}
-
-$temp_dir = $temp_dir.$random_number.$ii;
-
-while ( -e $temp_dir ) {
- $ii++;
- $temp_dir = $temp_dir.$random_number.$ii;
-}
-
-mkdir( $temp_dir, 0700 )
-|| die "\n\n$0: Could not create <<$temp_dir>>: $!\n\n";
-
-unless ( ( -e $temp_dir ) && ( -d $temp_dir ) ) {
- die "\n\n$0: <<$temp_dir>> does not exist, or is not a directory.\n\n";
-}
-
-
-&cp( $infile, $temp_dir."/INFILE" );
-unless ( chmod ( 0600, $temp_dir."/INFILE" ) ) {
- warn "\n\n$0: Could not chmod. $!\n\n";
-}
-$infile = "INFILE";
-
-chdir ( $temp_dir )
-|| die "\n\n$0: Could not chdir to <<$temp_dir>>: $!\n\n";
-
-&cp( $infile, "infile" );
-@out = &getNumberOfSeqsAndAas( $infile );
-$number_of_seqs = $out[ 0 ];
-$number_of_aa = $out[ 1 ];
-
-my $SEQBOOT_OUTFILE = "seqboot_outfile";
-
-if ( $bootstraps > 1 && ( $use_pwd_based_methods == 1
- || $use_phyml == 1
- || $use_proml == 1
- || $use_protpars == 1 ) ) {
- &executeSeqboot( $seed, $bootstraps );
- &mv( "outfile", $SEQBOOT_OUTFILE );
- &rm( "infile" );
-}
-
-&cp( $infile, "align" );
-
-if ( $use_pwd_based_methods == 1 ) {
- # Calculating the pairwise distances (saved in file "infile"): "puzzle"
- if ( $bootstraps > 1 ) {
- &executePuzzleBootstrapped( $SEQBOOT_OUTFILE,
- $matrix,
- $number_of_seqs,
- $exact_parameter_est,
- $rate_heterogeneity );
-
- $pwdfile = $SEQBOOT_OUTFILE.".dist";
- }
- else {
- &executePuzzle( "infile",
- $matrix,
- $number_of_seqs,
- $exact_parameter_est,
- $rate_heterogeneity );
- $pwdfile = "infile.dist";
- }
-}
-
-&rm( "infile" );
-
-# Methods based on alignment
-# --------------------------
-my $OUTTREE_RAXML = "outtree_rax";
-my $OUTTREE_PHYML = "outtree_phyml";
-my $OUTTREE_PROML = "outtree_proml";
-my $OUTTREE_PROTPARS = "outtree_protpars";
-
-my $CONSENSUS_RAXML = "consensus_raxml";
-my $CONSENSUS_PHYML = "consensus_phyml";
-my $CONSENSUS_PROML = "consensus_proml";
-my $CONSENSUS_PROTPARS = "consensus_protpars";
-
-my $OUTTREES_ALL = "outtrees_all";
-my $all_count = 0;
-
-if ( $use_raxml == 1 ) {
-
- my $model = "---";
- if ( $matrix == 0 ) {
- $model = "JTT";
- }
- elsif ( $matrix == 1 ) {
- $model = "DAYHOFF";
- }
- elsif ( $matrix == 2 ) {
- $model = "BLOSUM62";
- }
- elsif ( $matrix == 3 ) {
- $model = "MTREV";
- }
- elsif ( $matrix == 5 ) {
- $model = "VT";
- }
- elsif ( $matrix == 6 ) {
- $model = "WAG";
- }
- elsif ( $matrix == 7 ) {
- $model = "VT";
- }
- elsif ( $matrix == 8 ) {
- $model = "DCMUT";
- }
- else {
- &dieWithUnexpectedError( "Unknown model: matrix=$matrix" );
- }
-
- print( "\n========== RAxML begin =========\n\n" );
- # Six arguments:
- # 1. DNA or Amino-Acids sequence filename (PHYLIP format)
- # 2. Model, eg. PROTGAMMAIVT
- # 3. Replicates (bootstrap)
- # 4. Seed for bootstrap
- # 5. Output suffix
- # 6. Algorithm (only for bootstrap, default otherwise)
- my $invar = "";
- if ( $estimate_invar_sites == 1 ) {
- $invar = "I";
- }
-
- # NOTE. RaxML does its own bootstrapping.
- &executeRaxml( "align", $RAXML_MODEL_BASE.$invar.$model."F", $bootstraps, $seed, "xxx", $RAXML_ALGORITHM );
- print( "\n========== RAxML end =========\n\n" );
-
- &rm( "RAxML_log.xxx" );
- &rm( "RAxML_parsimonyTree.xxx" );
- &mv( "RAxML_info.xxx", $outfile."_raxml_info" );
- if ( $bootstraps > 1 ) {
- &rm( "RAxML_bestTree.xxx" );
- &mv( "RAxML_bipartitions.xxx", $CONSENSUS_RAXML );
- &append( "RAxML_bootstrap.xxx", $OUTTREES_ALL );
- if ( $keep_multiple_trees == 1 ) {
- &mv( "RAxML_bootstrap.xxx", $multitreefile_raxml );
- }
- else {
- &rm( "RAxML_bootstrap.xxx" );
- }
- $all_count++;
- }
- else {
- &mv( "RAxML_result.xxx", $OUTTREE_RAXML );
- }
-}
-
-
-if ( $use_phyml == 1 ) {
-
- my $model = "---";
- if ( $matrix == 0 ) {
- $model = "JTT";
- }
- elsif ( $matrix == 1 ) {
- $model = "Dayhoff";
- }
- elsif ( $matrix == 2 ) {
- $model = "Blosum62";
- }
- elsif ( $matrix == 3 ) {
- $model = "MtREV";
- }
- elsif ( $matrix == 5 ) {
- $model = "VT";
- }
- elsif ( $matrix == 6 ) {
- $model = "WAG";
- }
- elsif ( $matrix == 7 ) {
- $model = "VT";
- }
- elsif ( $matrix == 8 ) {
- $model = "DCMut";
- }
- else {
- &dieWithUnexpectedError( "Unknown model: matrix=$matrix" );
- }
-
- my $input = "";
- if ( $bootstraps > 1 ) {
- $input = $SEQBOOT_OUTFILE;
- }
- else {
- $input = "align";
- }
- print( "\n========== PHYML begin =========\n\n" );
- # Six arguments:
- # 1. DNA or Amino-Acids sequence filename (PHYLIP format)
- # 2. number of data sets to analyse (ex:3)
- # 3. Model: JTT | MtREV | Dayhoff | WAG | VT | DCMut | Blosum62 (Amino-Acids)
- # 4. number of relative substitution rate categories (ex:4), positive integer
- # 5. starting tree filename (Newick format), your tree filename | BIONJ for a distance-based tree
- # 6. 1 to estimate proportion of invariable sites, otherwise, fixed proportion "0.0" is used
- # PHYML produces several results files :
- # <sequence file name>_phyml_lk.txt : likelihood value(s)
- # <sequence file name>_phyml_tree.txt : inferred tree(s)
- # <sequence file name>_phyml_stat.txt : detailed execution stats
- &executePhyml( $input, $bootstraps, $model, $phyml_rel_substitution_rate_cat, "BIONJ", $estimate_invar_sites );
- print( "\n========== PHYML end =========\n\n" );
-
- &rm( $input."_phyml_lk.txt" );
- &mv( $input."_phyml_tree.txt", $OUTTREE_PHYML );
- if ( -e $outfile."_phyml_stat" ) {
- &rm( $outfile."_phyml_stat" );
- }
- &mv( $input."_phyml_stat.txt", $outfile."_phyml_stat" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_PHYML, $OUTTREES_ALL );
- $all_count++;
- }
-
-}
-
-if ( $use_proml == 1 ) {
- my $input = "";
- if ( $bootstraps > 1 ) {
- $input = $SEQBOOT_OUTFILE;
- }
- else {
- $input = "align";
- }
- print( "\n========== PHYLIP PROML begin =========\n\n" );
- # Five arguments:
- # 1. name of alignment file (in correct format!)
- # 2. number of bootstraps
- # 3. jumbles: 0: do not jumble; >=1 number of jumbles
- # 4. seed for random number generator
- # 5. 1 for PAM instead of JTT
- my $use_pam = 1;
- if ( $matrix == 0 ) {
- $use_pam = 0;
- }
- &executeProml( $input, $bootstraps, $jumbles, $seed, $use_pam, $use_global_rearr );
- print( "\n========== PHYLIP PROML end =========\n\n" );
- &mv( "outtree", $OUTTREE_PROML );
- &rm( "outfile" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_PROML, $OUTTREES_ALL );
- $all_count++;
- }
-}
-
-
-if ( $use_protpars == 1 ) {
- my $input = "";
- if ( $bootstraps > 1 ) {
- $input = $SEQBOOT_OUTFILE;
- }
- else {
- $input = "align";
- }
- print( "\n========== PHYLIP PROTPARS begin =========\n\n" );
- &executeProtpars( $input, $bootstraps, $jumbles, $seed );
- print( "\n========== PHYLIP PROTPARS end =========\n\n" );
- &mv( "outtree", $OUTTREE_PROTPARS );
- &rm( $outfile."_protpars_outfile" );
- &mv( "outfile", $outfile."_protpars_outfile" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_PROTPARS, $OUTTREES_ALL );
- $all_count++;
- }
-}
-
-
-
-# Methods based on PWD
-# --------------------
-my $OUTTREE_FASTME = "outtree_fastme";
-my $OUTTREE_PHYLIP_NJ = "outtree_phylip_nj";
-my $OUTTREE_PHYLIP_FM = "outtree_phylip_fm";
-my $OUTTREE_PHYLIP_ME = "outtree_phylip_me";
-my $OUTTREE_BIONJ = "outtree_bionj";
-my $OUTTREE_WEIGHBOR = "outtree_weighbor";
-
-my $CONSENSUS_FASTME = "consensus_fastme";
-my $CONSENSUS_PHYLIP_NJ = "consensus_phylip_nj";
-my $CONSENSUS_PHYLIP_FM = "consensus_phylip_fm";
-my $CONSENSUS_PHYLIP_ME = "consensus_phylip_me";
-my $CONSENSUS_BIONJ = "consensus_bionj";
-my $CONSENSUS_WEIGHBOR = "consensus_weighbor";
-my $CONSENSUS_ALL = "consensus_all";
-
-
-if ( $use_fastme == 1 ) {
- print( "\n========== FASTME begin =========\n\n" );
- &executeFastme( $pwdfile, $bootstraps, $fastme_init_tree_opt );
- print( "\n========== FASTME end ===========\n\n" );
- &rm( "output.d" );
- &mv( "output.tre", $OUTTREE_FASTME );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_FASTME, $OUTTREES_ALL );
- $all_count++;
- }
-}
-if ( $use_phylip_nj ) {
- print( "\n========== PHYLIP NEIGHBOR begin =========\n\n" );
- &executeNeighbor( $pwdfile, $bootstraps, $seed, 0 );
- print( "\n========== PHYLIP NEIGHBOR end =========\n\n" );
- &mv( "outtree", $OUTTREE_PHYLIP_NJ );
- &rm( "outfile" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_PHYLIP_NJ, $OUTTREES_ALL );
- $all_count++;
- }
-}
-if ( $use_phylip_fitch_fm ) {
- print( "\n========== PHYLIP FITCH FM begin =========\n\n" );
- &executeFitch( $pwdfile, $bootstraps, $seed, $jumbles, 0, "FM", $use_global_rearr );
- print( "\n========== PHYLIP FITCH FM end =========\n\n" );
- &mv( "outtree", $OUTTREE_PHYLIP_FM );
- &rm( "outfile" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_PHYLIP_FM, $OUTTREES_ALL );
- $all_count++;
- }
-}
-if ( $use_phylip_fitch_me ) {
- print( "\n========== PHYLIP FITCH ME begin =========\n\n" );
- &executeFitch( $pwdfile, $bootstraps, $seed, $jumbles, 0, "ME", $use_global_rearr );
- print( "\n========== PHYLIP FITCH ME end =========\n\n" );
- &mv( "outtree", $OUTTREE_PHYLIP_ME );
- &rm( "outfile" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_PHYLIP_ME, $OUTTREES_ALL );
- $all_count++;
- }
-}
-if ( $use_bionj ) {
- print( "\n========== BIONJ begin =========\n\n" );
- &executeBionj( $pwdfile, $OUTTREE_BIONJ );
- print( "\n========== BIONJ end =========\n\n" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_BIONJ, $OUTTREES_ALL );
- $all_count++;
- }
-}
-if ( $use_weighbor ) {
- print( "\n========== WEIGHBOR begin =========\n\n" );
- &executeWeighbor( $number_of_aa, 14, $pwdfile, $OUTTREE_WEIGHBOR );
- print( "\n========== WEIGHBOR end =========\n\n" );
- if ( $bootstraps > 1 ) {
- &append( $OUTTREE_WEIGHBOR, $OUTTREES_ALL );
- $all_count++;
- }
-}
-
-
-
-if ( $bootstraps > 1 ) {
- # Consense:
- if ( $use_fastme == 1 ) {
- &consense( $OUTTREE_FASTME, $CONSENSUS_FASTME );
- }
- if ( $use_phylip_nj == 1 ) {
- &consense( $OUTTREE_PHYLIP_NJ, $CONSENSUS_PHYLIP_NJ );
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- &consense( $OUTTREE_PHYLIP_FM, $CONSENSUS_PHYLIP_FM );
- }
- if ( $use_phylip_fitch_me == 1 ) {
- &consense( $OUTTREE_PHYLIP_ME, $CONSENSUS_PHYLIP_ME );
- }
- if ( $use_bionj == 1 ) {
- &consense( $OUTTREE_BIONJ, $CONSENSUS_BIONJ );
- }
- if ( $use_weighbor == 1 ) {
- &consense( $OUTTREE_WEIGHBOR, $CONSENSUS_WEIGHBOR );
- }
- if ( $use_phyml == 1 ) {
- &consense( $OUTTREE_PHYML, $CONSENSUS_PHYML );
- }
- if ( $use_proml == 1 ) {
- &consense( $OUTTREE_PROML, $CONSENSUS_PROML );
- }
- if ( $use_protpars == 1 ) {
- &consense( $OUTTREE_PROTPARS, $CONSENSUS_PROTPARS );
- }
- if ( $all_count > 1 ) {
- &consense( $OUTTREES_ALL, $CONSENSUS_ALL );
- }
- else {
- &rm( $OUTTREES_ALL );
- }
-
- my $INTREE_FOR_PUZZLE = "intree"; #why so serious?
- &rm( $INTREE_FOR_PUZZLE );
- system( "touch", $INTREE_FOR_PUZZLE )
- && die("\n\n$0: could not \"touch $INTREE_FOR_PUZZLE\": $!\n\n");
-
- if ( $use_fastme == 1 ) {
- &append( $CONSENSUS_FASTME, $INTREE_FOR_PUZZLE );
- }
- if ( $use_phylip_nj == 1 ) {
- &append( $CONSENSUS_PHYLIP_NJ, $INTREE_FOR_PUZZLE );
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- &append( $CONSENSUS_PHYLIP_FM, $INTREE_FOR_PUZZLE );
- }
- if ( $use_phylip_fitch_me == 1 ) {
- &append( $CONSENSUS_PHYLIP_ME, $INTREE_FOR_PUZZLE );
- }
- if ( $use_bionj == 1 ) {
- &append( $CONSENSUS_BIONJ, $INTREE_FOR_PUZZLE );
- }
- if ( $use_weighbor == 1 ) {
- &append( $CONSENSUS_WEIGHBOR, $INTREE_FOR_PUZZLE );
- }
- if ( $use_raxml == 1 ) {
- # Needed, because TREE-PUZZLE adds internal labels for all subsequent trees
- # when evaluating given trees (this seems a strange behaviour).
- removeSupportValues( $CONSENSUS_RAXML, $CONSENSUS_RAXML."_support_removed" );
- &append( $CONSENSUS_RAXML."_support_removed", $INTREE_FOR_PUZZLE );
- &rm( $CONSENSUS_RAXML."_support_removed" );
- }
- if ( $use_phyml == 1 ) {
- &append( $CONSENSUS_PHYML, $INTREE_FOR_PUZZLE );
- }
- if ( $use_proml == 1 ) {
- &append( $CONSENSUS_PROML, $INTREE_FOR_PUZZLE );
- }
- if ( $use_protpars == 1 ) {
- &append( $CONSENSUS_PROTPARS, $INTREE_FOR_PUZZLE );
- }
- if ( $all_count > 1 ) {
- &append( $CONSENSUS_ALL, $INTREE_FOR_PUZZLE );
- }
-
-
- # Puzzle for ML branch lenghts:
- # The alignment is read from infile by default.
- # The tree is read from intree by default.
- &rm( "infile" );
- &mv( "align", "infile" ); # align = original alignment in phylip interleaved.
-
- &executePuzzleToCalculateBranchLenghts( $matrix,
- $exact_parameter_est,
- $rate_heterogeneity );
-
- my $OUTTREE_PUZZLE = "outtree_puzzle";
-
- &rm( $outfile."_puzzle_outfile" );
-
- &mv( "outfile", $outfile."_puzzle_outfile" );
- &mv( "outtree", $OUTTREE_PUZZLE );
- &rm( "outdist" );
- &rm( "intree" );
-
-
- # Transfer
- # --------
- my $counter = 0;
- if ( $use_fastme == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_FASTME, $fastme_outtree, $counter++ );
- &rm( $CONSENSUS_FASTME );
- }
- if ( $use_phylip_nj == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYLIP_NJ, $phylip_nj_outtree, $counter++ );
- &rm( $CONSENSUS_PHYLIP_NJ );
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYLIP_FM, $phylip_fm_outtree, $counter++ );
- &rm( $CONSENSUS_PHYLIP_FM );
- }
- if ( $use_phylip_fitch_me == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYLIP_ME, $phylip_me_outtree, $counter++ );
- &rm( $CONSENSUS_PHYLIP_ME );
- }
- if ( $use_bionj == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_BIONJ, $bionj_outtree, $counter++ );
- &rm( $CONSENSUS_BIONJ );
- }
- if ( $use_weighbor == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_WEIGHBOR, $weighbor_outtree, $counter++ );
- &rm( $CONSENSUS_WEIGHBOR );
- }
- if ( $use_raxml == 1 ) {
- &to_phyloxml( $CONSENSUS_RAXML, $raxml_outtree, 1, 1 );
- $counter++;
- }
- if ( $use_phyml == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PHYML, $phyml_outtree, $counter++ );
- &rm( $CONSENSUS_PHYML );
- }
- if ( $use_proml == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PROML, $proml_outtree, $counter++ );
- &rm( $CONSENSUS_PROML );
- }
- if ( $use_protpars == 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_PROTPARS, $protpars_outtree, $counter++ );
- &rm( $CONSENSUS_PROTPARS );
- }
- if ( $all_count > 1 ) {
- &executeSupportTransfer( $OUTTREE_PUZZLE, $CONSENSUS_ALL, $all_outtree, $counter++ );
- &rm( $CONSENSUS_ALL );
- }
-
- # Clean up
- # --------
- &rm( $OUTTREE_PUZZLE );
- &rm( $SEQBOOT_OUTFILE );
- if ( $keep_multiple_trees == 1 ) {
- if ( $use_fastme == 1 ) {
- &mv( $OUTTREE_FASTME, $multitreefile_fastme );
- }
- if ( $use_phylip_nj == 1 ) {
- &mv( $OUTTREE_PHYLIP_NJ, $multitreefile_phylip_nj );
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- &mv( $OUTTREE_PHYLIP_FM, $multitreefile_phylip_fm );
- }
- if ( $use_phylip_fitch_me == 1 ) {
- &mv( $OUTTREE_PHYLIP_ME, $multitreefile_phylip_me );
- }
- if ( $use_bionj == 1 ) {
- &mv( $OUTTREE_BIONJ, $multitreefile_bionj );
- }
- if ( $use_weighbor == 1 ) {
- &mv( $OUTTREE_WEIGHBOR, $multitreefile_weighbor );
- }
- if ( $use_phyml == 1 ) {
- &mv( $OUTTREE_PHYML, $multitreefile_phyml );
- }
- if ( $use_proml == 1 ) {
- &mv( $OUTTREE_PROML, $multitreefile_proml );
- }
- if ( $use_protpars == 1 ) {
- &mv( $OUTTREE_PROTPARS, $multitreefile_protpars );
- }
- &mv( $pwdfile, $multipwdfile );
- }
- else {
- if ( $use_fastme == 1 ) {
- &rm( $OUTTREE_FASTME );
- }
- if ( $use_phylip_nj == 1 ) {
- &rm( $OUTTREE_PHYLIP_NJ );
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- &rm( $OUTTREE_PHYLIP_FM );
- }
- if ( $use_phylip_fitch_me == 1 ) {
- &rm( $OUTTREE_PHYLIP_ME );
- }
- if ( $use_bionj == 1 ) {
- &rm( $OUTTREE_BIONJ );
- }
- if ( $use_weighbor == 1 ) {
- &rm( $OUTTREE_WEIGHBOR );
- }
- if ( $use_phyml == 1 ) {
- &rm( $OUTTREE_PHYML );
- }
- if ( $use_proml == 1 ) {
- &rm( $OUTTREE_PROML );
- }
- if ( $use_protpars == 1 ) {
- &rm( $OUTTREE_PROTPARS );
- }
- &rm( $pwdfile );
- }
- if ( $all_count > 1 ) {
- &rm( $OUTTREES_ALL );
- }
-} # if ( $bootstraps > 1 )
-else {
- &rm( "infile.dist" );
-
- &rm( "infile.puzzle" );
- if ( $use_fastme == 1 ) {
- &to_phyloxml( $OUTTREE_FASTME, $fastme_outtree, 0, 1 );
- }
- if ( $use_phylip_nj == 1 ) {
- &to_phyloxml( $OUTTREE_PHYLIP_NJ, $phylip_nj_outtree, 0, 1);
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- &to_phyloxml( $OUTTREE_PHYLIP_FM, $phylip_fm_outtree, 0, 1 );
- }
- if ( $use_phylip_fitch_me == 1 ) {
- &to_phyloxml( $OUTTREE_PHYLIP_ME, $phylip_me_outtree, 0, 1 );
- }
- if ( $use_bionj == 1 ) {
- &to_phyloxml( $OUTTREE_BIONJ, $bionj_outtree, 0, 1 );
- }
- if ( $use_weighbor == 1 ) {
- &to_phyloxml( $OUTTREE_WEIGHBOR, $weighbor_outtree, 0, 1 );
- }
- if ( $use_raxml == 1 ) {
- &to_phyloxml( $OUTTREE_RAXML, $raxml_outtree, 0, 1 );
- }
- if ( $use_phyml == 1 ) {
- &to_phyloxml( $OUTTREE_PHYML, $phyml_outtree, 0, 1 );
- }
- if ( $use_proml == 1 ) {
- &to_phyloxml( $OUTTREE_PROML, $proml_outtree, 0, 1 );
- }
- if ( $use_protpars == 1 ) {
- &to_phyloxml( $OUTTREE_PROTPARS, $protpars_outtree, 0, 1 );
- }
-} # if ( $bootstraps > 1 )
-
-&rm( $infile );
-&rm( "infile" );
-&rm( "align" );
-&rm( "align.reduced" );
-
-
-$log = $log."Finish time/date : ".`date`;
-
-if ( $bootstraps > 1 ) {
- $log = $log."Puzzle output file : ".$outfile."_puzzle_outfile\n";
-}
-$log = $log."Columns in alignment : $number_of_aa\n";
-$log = $log."Number of sequences in alignment : $number_of_seqs\n";
-if ( $all_count > 1 ) {
- $log = $log."Combined consensus : $all_outtree\n";
-}
-
-
-if ( $bootstraps > 1 ) {
- $log = $log."\n\n";
- $log = $log."Simple support value statistics (trees are numbered the same as for TREE PUZZLE output)\n";
- $log = $log."------------------------------- \n";
- $log = $log."\n";
-}
-
-open( OUT, ">$logfile" ) || die "\n$0: Cannot create file <<$logfile>>: $!\n";
-print OUT $log;
-close( OUT );
-
-if ( $bootstraps > 1 ) {
- # Simple support statistics
- # -------------------------
- my $SS_OUT = $temp_dir."/ss_out";
- my @phylos = ();
- my $ounter = 0;
- if ( $use_fastme == 1 ) {
- $phylos[ $ounter++ ] = $fastme_outtree;
- }
- if ( $use_phylip_nj == 1 ) {
- $phylos[ $ounter++ ] = $phylip_nj_outtree;
- }
- if ( $use_phylip_fitch_fm == 1 ) {
- $phylos[ $ounter++ ] = $phylip_fm_outtree;
- }
- if ( $use_phylip_fitch_me == 1 ) {
- $phylos[ $ounter++ ] = $phylip_me_outtree;
- }
- if ( $use_bionj == 1 ) {
- $phylos[ $ounter++ ] = $bionj_outtree;
- }
- if ( $use_weighbor == 1 ) {
- $phylos[ $ounter++ ] = $weighbor_outtree;
- }
- if ( $use_raxml == 1 ) {
- $phylos[ $ounter++ ] = $raxml_outtree;
- }
- if ( $use_phyml == 1 ) {
- $phylos[ $ounter++ ] = $phyml_outtree;
- }
- if ( $use_proml == 1 ) {
- $phylos[ $ounter++ ] = $proml_outtree;
- }
- if ( $use_protpars == 1 ) {
- $phylos[ $ounter++ ] = $protpars_outtree;
- }
- if ( $all_count > 1) {
- $phylos[ $ounter++ ] = $all_outtree;
- }
- &executeSupportStatistics( $SS_OUT, @phylos );
- &append( $SS_OUT, $logfile );
- &rm( $SS_OUT );
-
- # Append parts of puzzle output file
- # ----------------------------------
- if ( $all_count > 1 ) {
- &parsePuzzleOutfile( $outfile."_puzzle_outfile", $logfile );
- }
-}
-
-chdir( $current_dir )
-|| die "\n\n$0: Could not chdir to <<$current_dir>>: $!\n\n";
-
-rmdir( $temp_dir )
-|| print "\n\n$0: Warning: Could not remove <<$temp_dir>>: $!\n\n";
-
-print "\n\n\n$0 successfully comleted.\n\n";
-
-exit( 0 );
-
-
-
-# Methods:
-# --------
-
-
-# Six arguments:
-# 1. DNA or Amino-Acids sequence filename (PHYLIP format)
-# 2. Model, eg. PROTGAMMAIVT
-# 3. Replicates (bootstrap)
-# 4. Seed for bootstrap
-# 5. Output suffix
-# 6. Algorithm (only for bootstrap, default otherwise)
-# NOTE. RaxML does its own bootstrapping.
-sub executeRaxml {
- my $msa = $_[ 0 ];
- my $model = $_[ 1 ];
- my $replicates = $_[ 2 ];
- my $seed = $_[ 3 ];
- my $outfile_suffix = $_[ 4 ];
- my $algo = $_[ 5 ];
-
- &testForTextFilePresence( $msa );
- my $command = "$RAXML -m $model -s $msa -n $outfile_suffix";
-
- if ( $replicates > 1 ) {
- $command = $command . " -x $seed -N $replicates";
- if ( $algo ) {
- $command = $command . " -f $algo";
- }
- }
-
- print( "\n$command\n");
-
- system( $command )
- && &dieWithUnexpectedError( $command );
-
-}
-
-
-sub to_phyloxml {
- my $from = $_[ 0 ];
- my $to = $_[ 1 ];
- my $internal_names_are_boots = $_[ 2 ];
- my $extract_taxonomy = $_[ 3 ];
- &dieIfFileExists( $to );
- &dieIfFileNotExists( $from );
- my $command = "$NEWICK_TO_PHYLOXML -f=nn $from $to";
- if ( $internal_names_are_boots == 1 ) {
- $command = $command . " -i";
- }
- if ( $extract_taxonomy == 1 ) {
- $command = $command . " -xt";
- }
- system( $command )
- && die "$0: Could not execute \"$command \"";
- &rm( $from );
-}
-
-
-sub mv {
- my $from = $_[ 0 ];
- my $to = $_[ 1 ];
- &dieIfFileExists( $to );
- &dieIfFileNotExists( $from );
- system( "mv", $from, $to )
- && die "\n\n$0: could not move \"$from\" to \"$to\": $!\n\n";
-}
-
-sub cp {
- my $from = $_[ 0 ];
- my $to = $_[ 1 ];
- &dieIfFileExists( $to );
- &dieIfFileNotExists( $from );
-
- system( "cp", $from, $to )
- && die "\n\n$0: could not copy \"$from\" to \"$to\": $!\n\n";
-}
-
-sub rm {
- my $f = $_[ 0 ];
- unlink( $f );
-}
-
-sub consense {
- my $multi_in = $_[ 0 ];
- my $consense_out = $_[ 1 ];
- &executeConsense( $multi_in );
- &mv( "outtree", $consense_out );
- &rm( "outfile" );
-
-}
-
-
-
-# 1. file to be appended
-# 2. file to append to
-sub append {
- my $to_be_appended = $_[ 0 ];
- my $append_to = $_[ 1 ];
- &dieIfFileNotExists( $to_be_appended );
- system( "cat $to_be_appended >> $append_to" )
- && die "\n\n$0: could not execute \"cat $to_be_appended >> $append_to\": $!\n\n";
-
-}
-
-sub dieIfFileExists {
- my $file = $_[ 0 ];
- if ( -e $file ) {
- die "\n\n$0: \"$file\" already exists\n\n";
- }
-}
-
-sub dieIfFileNotExists {
- my $file = $_[ 0 ];
- unless ( ( -s $file ) && ( -f $file ) ) {
- die( "\n\n$0: \"$file\" does not exist or is empty" );
- }
-}
-
-
-
-
-# Two arguments:
-# 1. seed for random number generator
-# 2. number of bootstraps
-# Reads in "infile" by default.
-sub executeSeqboot {
-
- my $s = $_[ 0 ];
- my $bs = $_[ 1 ];
- my $verb = "";
-
- &testForTextFilePresence( $infile );
-
-
- $verb = "
-2";
-
- system( "$SEQBOOT << !
-r
-$bs$verb
-Y
-$s
-!" )
- && die "$0: Could not execute \"$SEQBOOT\"";
-
-}
-
-
-
-# One/two/three argument(s):
-# Reads in tree from "intree" by default. (Presence of "intree" automatically
-# switches into "User defined trees" mode.)
-# 1. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
-# 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
-# 2. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
-# 3. Model of rate heterogeneity:
-# 1 for "8 Gamma distributed rates"
-# 2 for "Two rates (1 invariable + 1 variable)"
-# 3 for "Mixed (1 invariable + 8 Gamma rates)"
-# otherwise: Uniform rate
-# Last modified: 09/08/03 (added 2nd and 3rd parameter)
-sub executePuzzleToCalculateBranchLenghts {
- my $matrix_option = $_[ 0 ];
- my $parameter_estimates_option = $_[ 1 ];
- my $rate_heterogeneity_option = $_[ 2 ];
- my $i = 0;
- my $mat = "";
- my $est = "";
- my $rate = "";
-
- unless ( ( -s "infile" ) && ( -f "infile" ) && ( -T "infile" ) ) {
- die "\n$0: executePuzzleToCalculateBranchLenghts: <<infile>> does not exist, is empty, or is not a plain textfile.\n";
- }
- unless ( ( -s "intree" ) && ( -f "intree" ) && ( -T "intree" ) ) {
- die "\n$0: executePuzzleToCalculateBranchLenghts: <<intree>> does not exist, is empty, or is not a plain textfile.\n";
- }
-
- $mat = setModelForPuzzle( $matrix_option );
- if ( $parameter_estimates_option ) {
- $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
- }
- if ( $rate_heterogeneity_option ) {
- $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
- }
-
-
- system( "$PUZZLE << !
-$mat$est$rate
-x
-y
-!" )
- && die "$0: Could not execute \"$PUZZLE\" (mat=$mat est=$est rate=$rate)";
-
-}
-
-# Two arguments:
-# 1. puzzle outfile
-# 2. file to append to
-sub parsePuzzleOutfile {
- my $puzzle_outfile = $_[ 0 ];
- my $file_to_append_to = $_[ 1 ];
- &testForTextFilePresence( $puzzle_outfile );
- open( OUT, ">>$file_to_append_to" ) || &dieWithUnexpectedError( "Cannot open \"$file_to_append_to\"" );
- open( IN, "$puzzle_outfile" ) || &dieWithUnexpectedError( "Cannot open file \"$puzzle_outfile\"" );
- my $return_line;
- my $read = 0;
- print OUT "\nTREE PUZZLE output\n";
- print OUT "------------------\n";
- while ( $return_line = <IN> ) {
- if ( $return_line =~/COMPARISON OF USER TREES/ ) {
- $read = 1;
- }
- elsif( $return_line =~/TIME STAMP/ ) {
- $read = 0;
- }
- elsif( $read ) {
- print OUT $return_line;
- }
- }
- close( IN );
- close( OUT );
-}
-
-# Three/four arguments:
-# 1. Name of file containing tree with correct branch lengths
-# 2. Name of file containing tree with correct bootstraps
-# 3. Outputfilename
-# 4. Index of tree with correct branch lengths, in case more than one in file
-# Last modified: 2007.11.27
-sub executeSupportTransfer {
- my $tree_with_bl = $_[ 0 ];
- my $tree_with_bs = $_[ 1 ];
- my $out = $_[ 2 ];
- my $index = $_[ 3 ];
-
- &testForTextFilePresence( $tree_with_bl );
- &testForTextFilePresence( $tree_with_bs );
- my $command = "$SUPPORT_TRANSFER $tree_with_bl $tree_with_bs $out $index";
- system( $command )
- && die "$0: Could not execute \"$command\"";
-}
-
-
-# Two or more arguments:
-# 1. outfile
-# 2. phylogeny 1 with support values
-# 3. phylogeny 2 with support values
-# 4. ...
-sub executeSupportStatistics {
- my $outfile = $_[ 0 ];
- &dieIfFileExists( $outfile );
- my $phylos = "";
- for( my $i = 1; $i < scalar(@_); ++$i ) {
- &testForTextFilePresence( $_[ $i ] );
- $phylos .= $_[ $i ]." ";
- }
- my $command = "$SUPPORT_STATISTICS -o=$outfile $phylos";
- system( "$command" )
- && die "$0: Could not execute \"$command\"";
-}
-
-
-sub getNumberOfSeqsAndAas {
- my $infile = $_[ 0 ];
- my $seqs = 0;
- my $aa = 0;
- open( IN, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
- while( <IN> ) {
- if ( $_ =~ /^\s*(\d+)\s+(\d+)\s*$/ ) {
- $seqs = $1;
- $aa = $2;
- }
- }
- close( IN );
-
- if ( $seqs == 0 || $aa == 0 ) {
- die( "\n$0: Could not get number of seqs and aa from: $infile" );
- }
- return $seqs, $aa;
-}
-
-
-
-sub removeSupportValues {
- my $infile = $_[ 0 ];
- my $outfile = $_[ 1 ];
- &testForTextFilePresence( $infile );
- open( OUT, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
- open( IN, "$infile" ) || &dieWithUnexpectedError( "Cannot open file \"$infile\"" );
- while ( my $line = <IN> ) {
- $line =~ s/\)\d+\.?\d*:/\):/g;
- print OUT "$line";
- }
- close( OUT );
- close( IN );
-}
-
-
-
-
-# Six arguments:
-# 1. name of alignment file (in correct format!)
-# 2. number of bootstraps
-# 3. jumbles: 0: do not jumble; >=1 number of jumbles
-# 4. seed for random number generator
-# 5. 1 for PAM instead of JTT
-# 6. 1 to use globale rearragements
-sub executeProml {
- my $align = $_[ 0 ];
- my $bs = $_[ 1 ];
- my $rand = $_[ 2 ];
- my $s = $_[ 3 ];
- my $use_pam = $_[ 4 ];
- my $use_global_rearr = $_[ 5 ];
- my $jumble = "";
- my $multi = "";
- my $pam = "";
-
- &testForTextFilePresence( $align );
-
- if ( $bs > 1 && $rand < 1 ) {
- $rand = 1;
- }
-
- if ( $rand >= 1 ) {
- $jumble = "
-J
-$s
-$rand";
- }
-
- if ( $bs > 1 ) {
- $multi = "
-M
-D
-$bs";
- }
-
-
- if ( $use_pam == 1 ) {
- $pam = "
-P
-P";
- }
-
- my $global = "";
- if ( $use_global_rearr == 1 ) {
- $global = "
-G";
- }
-
- system( "$PROML 2>&1 << !
-$align$jumble$multi$pam$global
-3
-Y
-!" )
- && &dieWithUnexpectedError( "Could not execute \"$PROML $align$jumble$multi$pam$global\"" );
- # 3: Do NOT print out tree
-
- return;
-
-} ## executeProml
-
-
-sub printUsage {
-
- print <<END;
-
-Copyright (C) 2002-2007 Christian M. Zmasek
-All rights reserved
-
-Author: Christian M. Zmasek
-phylosoft\@gmail.com
-http://www.phylosoft.org
-
- Requirements phylo_pl is part of the FORESTER collection of programs.
- ------------ Many of its global variables are set via forester.pm.
-
- Note. Use xt.pl (for Pfam alignments) or mt.pl (for other alignments)
- to run phylo_pl.pl on whole directories of alignments files.
-
- Usage
- -----
-
- phylo_pl.pl [-options] <input alignment in SELEX (Pfam), PHYLIP
- sequential format, or Clustal W output> <outputfile>
- [path/name for temporary directory to be created]
-
- Example:
- "% phylo_pl.pl -B100q\@1nbS9X IL5.aln IL5_tree"
-
- Options
- -------
-
- Bx : Number of bootstraps. B0: do not bootstrap. Default is 100 bootstrapps.
- The number of bootstrapps should be divisible by 10.
- J : Use JTT matrix (Jones et al. 1992) in TREE-PUZZLE and/or PHYML, RAXML, default: VT (Mueller-Vingron 2000).
- L : Use BLOSUM 62 matrix (Henikoff-Henikoff 92) in TREE-PUZZLE and/or PHYML, RAXML, default: VT.
- M : Use mtREV24 matrix (Adachi-Hasegawa 1996) in TREE-PUZZLE and/or PHYML, default: VT.
- W : Use WAG matrix (Whelan-Goldman 2000) in TREE-PUZZLE and/or PHYML, RAXML, default: VT.
- P : Use PAM matrix (Dayhoff et al. 1978) in TREE-PUZZLE and/or PHYML, RAXML, default: VT.
- D : Use DCMut matrix (Kosial and Goldman, 2005) in PHYML, RAXML, VT in TREE-PUZZLE.
- A : Let TREE-PUZZLE choose which matrix to use, default: VT
- E : Exact parameter estimates in TREE-PUZZLE, default: Approximate.
- Model of rate heterogeneity in TREE-PUZZLE (default: Uniform rate):
- g : 8 Gamma distributed rates
- t : Two rates (1 invariable + 1 variable)
- m : Mixed (1 invariable + 8 Gamma rates)
- q\@x: Use FastME, x: 1: GME
- 2: BME
- 3: NJ
- n : Use PHYLIP Neighbor (NJ).
- f : Use PHYLIP Fitch.
- e : Use PHYLIP Minimal Evolution.
- b : Use BIONJ.
- w : Use Weighbor.
- x : Use RAxML.
- y : Use PHYML.
- o : Use PHYLIP proml.
- p : Use PHYLIP protpars.
- rx : Number of relative substitution rate categories in PHYML (default is 4).
- jx : Number of jumbles (input order randomization) for PHYLIP FM, ME, PROTPARS, and PROML (default is 2) (random seed set with Sx).
- I : Estimate proportion of invariable sites in RAXML and/or PHYML (otherwise, proportion "0.0" is used in PHYML)
- G : to turn on global rearrangements in PHYLIP FM, ME, and PROML
- Sx : Seed for random number generator(s). Must be 4n+1. Default is 9.
- X : To keep multiple tree file (=trees from bootstrap resampled alignments) and
- pairwise distance matrix file (in case of bootstrap analysis).
-
-END
-
-} ## printUsage
git://source.jalview.org / jalview.git / commitdiff