Mercurial > repos > ucsb-phylogenetics > osiris_phylogenetics

comparison phylographics/makeRtrees.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
2
3 #This script generates an R script to print trees to a pdf file
4 #input is a table with treename<tab>newick tree
5 use strict;
6
7 my $filename = $ARGV[0];
8 my $outfile = $ARGV[1];
9 open FILE, $filename or die $!;
10 my $treetype = $ARGV[2];
11 my $extiplabels = $ARGV[3];
12 my $options;
13 my $labeltaxfile = $ARGV[4];
14 my $query = $ARGV[5];
15 my $midroot = $ARGV[6];
16 my %labelhash;
17 my $genecount=0;
18 my @genes;
19 my $markquery;
20 my $midpoint;
21
22 if($query eq 'yes'){
23 $markquery = 1;
24 }else{
25 $markquery = 0;
26 }
27 unless($labeltaxfile eq 'None'){
28 open LABELFILE, $labeltaxfile or die $!;
29 while (<LABELFILE>) {
30 chomp;
31 #get a line from the data file
32 my $currentinput = "$_";
33 if($currentinput =~ /\t/){
34 my @splitline = split(/\t/);
35 my $speciesname= $splitline[0];
36 $speciesname = "'".$speciesname."'";
37 my $treename = $splitline[1];
38 if(exists $labelhash{$treename}){
39 push @{ $labelhash{$treename} }, $speciesname;
40 }else{
41 push @{ $labelhash{$treename} }, $speciesname;
42 #$labelhash{$treename} = $speciesname;
43 $genecount ++;
44 push @genes, $treename;
45 }
46 }
47 }
48
49 }#end unless
50
51
52 if($extiplabels eq 'yes'){
53 $options = ", show.tip.label=FALSE";
54 }else{
55 $options = ", show.tip.label=TRUE";
56 }
57
58 print "require(ape);\n";
59 print "require(phangorn);\n";
60 print "pdf(file='$outfile');\n";
61
62 while (<FILE>) {
63 chomp;
64 #get a line from the data file
65 my $currentinput = "$_";
66 my @splitline = split(/\t/);
67 my $treename= $splitline[0];
68 my $tree = $splitline[1];
69 my $labelsvector;
70
71 #print the R commands to make tree graphics
72 print "raw_tree <- read.tree(text = '$tree');\n";
73 print "raw_tree\$edge.length[ is.na(raw_tree\$edge.length) ] <- 0 \n";
74 if($midroot eq 'yes'){
75 print "raw_tree <- midpoint(raw_tree)\n";
76 }
77 #Check if large tree, then make text size smaller
78 print "thetips <- raw_tree\$tip.label \n";
79 print "numtips <- length(thetips) \n";
80
81 #Make text smaller for trees with many tips
82 print "if(numtips>250){plot(raw_tree, cex=0.15, edge.width = 0.1, type='$treetype' $options)}else if(numtips>75){plot(raw_tree, cex=0.3, type='$treetype' $options)}else{plot(raw_tree, cex=0.6, type='$treetype' $options)};\n";
83 print "title('Tree File: $treename');\n";
84
85 #Add taxon labels, if optional file present and if labels exist for tree
86 if(exists $labelhash{$treename}){
87 $labelsvector = join ",", @{ $labelhash{$treename} };
88 $labelsvector = "tolabel <- c(".$labelsvector.")";
89 print "thetips <- raw_tree\$tip.label \n";
90 print $labelsvector."\n";
91 print "labels <- match(tolabel,thetips) \n";
92 print "tiplabels(tip=labels, pch=21, cex=1) \n";
93 }
94
95
96 #Add taxon labels if gene name contains QUERY - for readplacement
97 if($markquery == 1){
98 print "thetips <- raw_tree\$tip.label \n";
99 print "qlabels <- grep(\'QUERY\',thetips) \n";
100 print "tiplabels(tip=qlabels, pch=21, cex=1.1) \n";
101 print "l1labels <- grep(\'LANDMARK1\',thetips) \n";
102 print "tiplabels(tip=l1labels, pch=15, cex=.8, col='red') \n";
103 }
104 }
105
106 print "dev.off();\n";
107 close FILE;
108
109 #Testing hash arrays
110 #my %nums;
111 #my $test='odd';
112 #for my $n (4,5,6,10) {
113 # if ($n % 2) {
114 # push @{ $nums{$test} }, $n;
115 # } else {
116 # push @{ $nums{even} }, $n;
117 # }
118 #}
119 #
120 #print join ', ', @{ $nums{even} };
121 #print "\n\n";