[Bio] / FigKernelScripts / make_proml_tree.pl Repository:
ViewVC logotype

View of /FigKernelScripts/make_proml_tree.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (download) (as text) (annotate)
Thu Feb 1 20:49:10 2007 UTC (12 years, 9 months ago) by golsen
Branch: MAIN
CVS Tags: mgrast_dev_08112011, rast_rel_2009_05_18, mgrast_dev_08022011, rast_rel_2014_0912, rast_rel_2008_06_18, myrast_rel40, rast_rel_2008_06_16, mgrast_dev_05262011, rast_rel_2008_12_18, mgrast_dev_04082011, rast_rel_2008_07_21, rast_rel_2010_0928, rast_2008_0924, mgrast_version_3_2, mgrast_dev_12152011, rast_rel_2008_04_23, mgrast_dev_06072011, rast_rel_2008_09_30, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, mgrast_rel_2008_0924, mgrast_rel_2008_1110_v2, rast_rel_2009_02_05, rast_rel_2011_0119, mgrast_rel_2008_0625, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, rast_rel_2008_10_09, mgrast_dev_04012011, rast_release_2008_09_29, mgrast_rel_2008_0806, mgrast_rel_2008_0923, mgrast_rel_2008_0919, rast_rel_2009_07_09, rast_rel_2010_0827, mgrast_rel_2008_1110, myrast_33, rast_rel_2011_0928, rast_rel_2008_09_29, mgrast_rel_2008_0917, rast_rel_2008_10_29, mgrast_dev_04052011, mgrast_dev_02222011, rast_rel_2009_03_26, mgrast_dev_10262011, rast_rel_2008_11_24, rast_rel_2008_08_07, HEAD
Changes since 1.5: +198 -33 lines
Replace make_proml_tree with entirely new version.  Move the former
version to make_proml_tree_0.

#! /usr/bin/perl
#
#   make_proml_tree_2:  A command line interface to the proml program
#
my $usage = <<"End_of_Usage";

Usage:

   make_proml_tree_2 < Alignment > Trees_with_likelihood_comment

options:
  -a Alpha         Alpha parameter for gamma distribution of rates
  -b GammaBins     Number of rate bins in gamma approximation (2 - 9)
  -c RatesCatFile  Rates and categories file (rates \\n categories)
  -d Directory     Directory for temp files. If it exists, files are retained.
  -g               Do global rearrangements
  -h HMMFile       User-defined rate HMM file [Not functional]
  -i FracInvar     Fraction of invariant sites
  -j JumbleSeed    Random addition order seed (odd)
  -l               User-defined branch lengths on user trees
  -m Model         Amino acid change model (JTT (D) | PMB | PAM)
  -n NJumbles      Number of random addition orders
  -p Persistance   Number of sites with correlated rates
  -s               Slower, more accurate calculation
  -t Tmp           Directory for directory of temporary files (D = /tmp)
  -u UserTreeFile  File of user trees
  -v CoefOfVar     Coeficient of variation for gamma distribution for rates
  -w Weights       Site weights file

End_of_Usage


use proml;
use gjoseqlib;
use gjonewicklib;

my %options = ( tree_format  => 'gjo' );

while ( @ARGV && $ARGV[0] =~ /^-/ )
{
    my $flag = shift @ARGV;
    my ( $file, $value, @trees );

    if    ( $flag =~ s/^-a// )
    {
        $value = $flag || shift @ARGV;
        $value > 0 or usage( "Bad value of alpha parameter: $value\n" );
        $options{ alpha } = $value;
    }
    elsif ( $flag =~ s/^-b// )
    {
        $value = $flag || shift @ARGV;
        $value > 0 && $value <= 9
            or usage( "Bad value for gamma bins parameter: $value" );
        $options{ gamma_bins } = $value;
    }
    elsif ( $flag =~ s/^-c// )
    {
        $file = $flag || shift @ARGV;
        -f $file
            and ( $value = read_categories_file( $file ) )
            or  usage( "Bad category file: $file" );
        $options{ categories } = $value;
    }
    elsif ( $flag =~ s/^-d// )
    {
        $value = $flag || shift @ARGV;
        $value
            or usage( "Bad value for directory for temp files: $value" );
        $options{ tmp_dir } = $value;
    }
    elsif ( $flag =~ s/^-g// )
    {
        $options{ global } = 1;
    }
    elsif ( $flag =~ s/^-h// )
    {
        $file = $flag || shift @ARGV;
        -f $file
            or  usage( "Bad HMM file: $file" );
    }
    elsif ( $flag =~ s/^-i// )
    {
        $value = $flag || shift @ARGV;
        $value >= 0 && $value < 1
            or usage( "Bad value for fraction of invariant sites: $value" );
        $options{ invar_frac } = $value;
    }
    elsif ( $flag =~ s/^-j// )
    {
        $value = $flag || shift @ARGV;
        $value > 0 && ( $value % 2 == 1 )
            or usage( "Bad value for jumble seed: $value" );
        $options{ jumble_seed } = $value;
    }
    elsif ( $flag =~ s/^-l// )
    {
        $options{ user_lengths } = 1;
    }
    elsif ( $flag =~ s/^-m// )
    {
        $value = $flag || shift @ARGV;
        $value
            or usage( "Bad value for evolution model: $value" );
        $options{ n_jumble } = $value;
    }
    elsif ( $flag =~ s/^-n// )
    {
        $value = $flag || shift @ARGV;
        $value > 0
            or usage( "Bad value for number of random addition orders: $value" );
        $options{ n_jumble } = $value;
    }
    elsif ( $flag =~ s/^-p// )
    {
        $value = $flag || shift @ARGV;
        $value > 0
            or usage( "Bad value for site rate correlation persistance length: $value" );
        $options{ persistance } = $value;
    }
    elsif ( $flag =~ s/^-s// )
    {
        $options{ slow } = 1;
    }
    elsif ( $flag =~ s/^-t// )
    {
        $value = $flag || shift @ARGV;
        -d $value
            or usage( "Bad value for tmp directory: $value" );
        $options{ tmp } = $value;
    }
    elsif ( $flag =~ s/^-u// )
    {
        $file = $flag || shift @ARGV;
        -f $file
            and ( @trees = gjonewicklib::read_newick_trees( $file ) )
            or  usage( "Bad tree file: $file" );
        $options{ user_trees } = \@trees;
    }
    elsif ( $flag =~ s/^-v// )
    {
        $value = $flag || shift @ARGV;
        $value >= 0
            or usage( "Bad value for coeficient of variation of gamma distribution: $value" );
        $options{ coef_of_var } = $value;
    }
    elsif ( $flag =~ s/^-w// )
    {
        my $file = $flag || shift @ARGV;
        -f $file
            and ( $value = read_weights_file( $file ) )
            or  usage( "Bad weights file: $file" );
        $options{ weights } = $value;
    }
    else
    {
        usage( "Bad flag: $flag" );
    }
}

my @align = gjoseqlib::read_fasta();

my @results = proml::proml( \@align, \%options );

foreach ( @results )
{
    my ( $tree, $likelihood ) = @$_;
    printf '[%12.4f]  ', $likelihood;
    gjonewicklib::writeNewickTree( $tree );
}

exit;


sub usage
{
    foreach ( @_ ) { print STDERR "$_\n" }
    print STDERR $usage;
    exit;
}


sub read_categories_file
{
    open( CATS, "<$_[0]" ) or return undef;
    $_ = <CATS>;
    chomp;
    s/[\t ,]+/ /g;
    my @rates = split;
    my $categories = join( '', map { chomp; s/\D+//g; $_ } <CATS> );
    close( CATS );
    [ \@rates, $categories ]
}


sub read_weights_file
{
    open( WGTS, "<$_[0]" ) or return undef;
    my $weights = join( '', map { chomp; s/\D+//g; $_ } <WGTS> );
    close( WGTS );
    $weights
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3