#!/usr/bin/perl -w my $usage = "\nUsage: $0 [-h] nexusFile\n". " -h: help\n". "Takes one nexus file, run model-test, and print out a paup-command \n". "file, which incorporate the model selected based on AIC.\n"; our ($opt_h); use File::Basename; use Getopt::Std; getopts('h') || die "$usage\n"; die "$usage\n" if (defined($opt_h)); die "ERROR: Please give only 1 file\n" if (@ARGV > 1); my $modeltestBin = "modeltest3.7"; my $mtBlock = "/usr/share/doc/modeltest-3.7/modelblockPAUPb10"; my $criterion = 'AIC'; # other choice is 'hLRT' #@ARGV = ('-') unless @ARGV; # take STDIN when no arg my $nex = shift; $base = basename($nex); if ($nex =~ /\.nex$/) { $base =~ s/\.nex$//; } elsif ($nex =~ /\.nx$/) { $base =~ s/\.nx$//; } system("mkdir -p modeltest.$base"); # run modeltest system("cd modeltest.$base; echo \"execute ../$nex; execute $mtBlock; quit WarnTSave=no;\" | paup 1>&2; $modeltestBin < model.scores > modeltest.out"); print "#NEXUS\nBEGIN PAUP;\nexecute $nex;\nset criterion=likelihood;\n"; my $encounter = ($criterion eq 'hLRT') ? 1 : 2; my ($cnt1, $cnt2) = (0,0); open MT, ") { if( /^\s*Model selected:\s*(.+)\s*/) { $cnt1++; if ($cnt1 == $encounter) { print "[! With $criterion, model selected: $1 ]\n" } } if (/^Lset/) { $cnt2++; if ($cnt2 == $encounter) { print "\n$_\n\n"; } } } close(MT); my $template = <<"EOF"; [ exhaustive search alltrees ; ] [ outgroup ''; ] [ heuristic search ] hsearch addseq=random nreps=10 rearrlimit=40000 limitperrep=yes; [ bootstrap nreps=50 treefile=boot.tre search=heuristic /addseq=random rearrlimit=3000 nreps=2 limitperrep=yes; ] [ save consensus tree contree all/strict=no majrule=yes treefile=consensus.tre; ] savetrees file=tree.tre root=y brlens=y format=NEXUS; savetrees file=tree.phy root=y brlens=y format=phylip; [ saveassum file=paup.assm; ] [ write brlen file ] [ log start file=brlens; describetrees /brlens=yes; log stop; ] End; EOF print "$template\n"; exit;