Mercurial > repos > ucsb-phylogenetics > osiris_phylogenetics

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First commit of old repositories
author osiris_phylogenetics <ucsb_phylogenetics@lifesci.ucsb.edu>
date Tue, 11 Mar 2014 12:19:13 -0700
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5b9a38ec4a39
1 #!/usr/bin/perl -w
2
3 use strict;
4
5 #use FindBin;
6 #use lib "$FindBin::Bin/lib";
7 use Bio::TreeIO;
8 use Bio::Tree::Tree;
9
10 ###this script will find the phylogenetic distance between two species
11 #input is a tree, output filename, and table with pairwise distances
12 #usage:
13 #PD.pl <pairsTable> <treefile> <outfile> <yes|no>
14 # parse in newick/new hampshire format
15 my @species1;
16 my @species2;
17
18
19 my $half=$ARGV[3];
20 my $divtimebool;
21 if($half eq 'yes'){
22 $divtimebool=1;
23 }elsif($half eq 'no'){
24 $divtimebool=0;
25 }else{
26 die "Argument must contain yes or no for divergence times\n";
27 }
28 my $outfile = $ARGV[2];
29 open(OUT, ">$outfile") or die("Couldn't open output file $ARGV[2]\n");
30
31
32 my $pairsfile = $ARGV[0];
33 open(PAIRS, "$pairsfile") or die("Couldn't open input file $ARGV[0]\n");
34 while (<PAIRS>) {
35 chomp;
36 my $sp1;
37 my $sp2;
38 ($sp1, $sp2) = split("\t");
39 push(@species1, $sp1);
40 push(@species2, $sp2);
41 }
42
43 my $treefile = $ARGV[1];
44
45 for(my $i=0; $i < @species1; $i++){
46 print OUT $species1[$i]."\t".$species2[$i];
47 open(TREE, "$treefile") or die("Couldn't open output file $ARGV[1]\n");
48
49 my $treeio = new Bio::TreeIO('-format' => 'newick',
50 '-file' => $treefile);
51
52 while(my $tree = $treeio->next_tree){;
53 my $node1 = $tree->find_node(-id => $species1[$i]);
54 my $node2 = $tree->find_node(-id => $species2[$i]);
55 my $distances = $tree->distance(-nodes => [$node1,$node2]);
56
57 #ADD OPTION FOR DIVIDING BY 2 FOR DIVERGENCE TIMES
58 if($divtimebool==1){
59 $distances = $distances/2 ;
60 }
61 print OUT "\t".$distances;
62 }
63 print OUT "\n";
64 close(TREE);
65 }
66
67 close(PAIRS);
68 close(OUT);