diff phylostatistics/PD.pl @ 0:5b9a38ec4a39 draft default tip

First commit of old repositories

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/phylostatistics/PD.pl	Tue Mar 11 12:19:13 2014 -0700
@@ -0,0 +1,68 @@
+#!/usr/bin/perl -w
+
+use strict;
+
+#use FindBin;
+#use lib "$FindBin::Bin/lib";
+use Bio::TreeIO;
+use Bio::Tree::Tree;
+
+###this script will find the phylogenetic distance between two species
+#input is a tree, output filename, and table with pairwise distances
+#usage:
+#PD.pl <pairsTable> <treefile> <outfile> <yes|no>
+# parse in newick/new hampshire format
+my @species1;
+my @species2;
+
+
+my $half=$ARGV[3];
+my $divtimebool;
+if($half eq 'yes'){
+	$divtimebool=1;
+}elsif($half eq 'no'){
+	$divtimebool=0;
+}else{
+	die "Argument must contain yes or no for divergence times\n";
+}
+my $outfile = $ARGV[2];
+open(OUT, ">$outfile") or die("Couldn't open output file $ARGV[2]\n");
+
+
+my $pairsfile = $ARGV[0];
+open(PAIRS, "$pairsfile") or die("Couldn't open input file $ARGV[0]\n");
+while (<PAIRS>) {
+        chomp;
+        my $sp1;
+        my $sp2;
+        ($sp1, $sp2) = split("\t");
+        push(@species1, $sp1);
+        push(@species2, $sp2);
+}
+
+my $treefile = $ARGV[1];
+
+for(my $i=0; $i < @species1; $i++){
+        print OUT $species1[$i]."\t".$species2[$i];
+        open(TREE, "$treefile") or die("Couldn't open output file $ARGV[1]\n");
+
+        my $treeio = new Bio::TreeIO('-format' => 'newick',
+                                   '-file'   => $treefile);
+
+        while(my $tree = $treeio->next_tree){;
+                my $node1 = $tree->find_node(-id => $species1[$i]);
+                my $node2 = $tree->find_node(-id => $species2[$i]);
+                my $distances = $tree->distance(-nodes => [$node1,$node2]);
+
+                #ADD OPTION FOR DIVIDING BY 2 FOR DIVERGENCE TIMES
+		if($divtimebool==1){
+                	$distances = $distances/2 ;
+		}
+                print OUT "\t".$distances;
+        }
+print OUT "\n";
+close(TREE);
+}
+
+close(PAIRS);
+close(OUT);