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

Annotation of /FigKernelScripts/place_properties_on_tree.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (view) (download) (as text)

1 : overbeek 1.5 ########################################################################
2 : overbeek 1.3 #
3 :     # This program takes as input a file describing the presence/absence of
4 :     # families in each of the genomes (the tips of the tree). It produces
5 :     # the families.on.tree file, which extends the presence/absence data to
6 :     # nonleaf nodes. The properties in that file are just
7 :     #
8 :     # [Property,Node,Values]
9 :     #
10 :     # where Values contains
11 :     #
12 : overbeek 1.4 # 0 - does not include the family
13 : overbeek 1.3 # 1 - does include the family
14 :     # 0,1 - uncertain
15 :     #
16 : overbeek 1.1 use Carp;
17 :     use strict;
18 :     use Data::Dumper;
19 :     use Getopt::Long;
20 :     use SeedEnv;
21 :     use tree_utilities;
22 :     my $usage = "usage: place_properties_on_tree -t tree -p property-table -e extended-property-table\n";
23 :     my $treeF;
24 :     my $propF;
25 :     my $extended_propF;
26 :    
27 :     my $rc = GetOptions('t=s' => \$treeF,
28 :     'p=s' => \$propF,
29 :     'e=s' => \$extended_propF);
30 :    
31 :     if ((! $rc) ||
32 :     (! -s $treeF) || (! -s $propF) || (! open(EX,">$extended_propF")))
33 :     {
34 :     print STDERR $usage; exit ;
35 :     }
36 :    
37 :     my $tree;
38 :     ($tree = &parse_newick_tree(join("",`cat $treeF`)))
39 :     || die "could not parse the tree";
40 : overbeek 1.5 my $tips = &tips_of_tree($tree);
41 : overbeek 1.1 open(PROP,"<$propF") || die "could not open $propF";
42 :     my %prop_occ;
43 : overbeek 1.5 my %genomes;
44 : overbeek 1.1 while (defined($_ = <PROP>))
45 :     {
46 :     chop;
47 :     my($prop,$tip,$val) = split(/\t/,$_);
48 : overbeek 1.5 $genomes{$tip} = 1;
49 : overbeek 1.1 $prop_occ{$prop}->{$tip} = $val;
50 :     }
51 :     close(PROP);
52 : overbeek 1.5 my $tree = &subtree($tree,\%genomes);
53 : overbeek 1.1 my @props = keys(%prop_occ);
54 :     foreach my $prop (@props)
55 :     {
56 :     my $occ_on_nodes = &place_prop_on_nodes($tree,$prop_occ{$prop});
57 :     foreach my $node (keys(%$occ_on_nodes))
58 :     {
59 :     my $vals = $occ_on_nodes->{$node};
60 :     my $set = join(",",@$vals);
61 :     print EX join("\t",($prop,$node,$set)),"\n";
62 :     }
63 :     }
64 :    
65 :     close(EX);
66 :    
67 :     sub place_prop_on_nodes {
68 :     my($tree,$occ_on_tips) = @_;
69 :    
70 :     my $occH = {};
71 :     foreach my $tip (keys(%$occ_on_tips))
72 :     {
73 : overbeek 1.6 my $val = defined($occ_on_tips->{$tip}) ? $occ_on_tips->{$tip} : 'unknown';
74 : overbeek 1.5 $occH->{$tip} = [$val];
75 : overbeek 1.1 }
76 :     &pass1($tree,$occH);
77 :     &pass2($tree,$occH);
78 :     return $occH;
79 :     }
80 :    
81 :     sub pass1 {
82 :     my($node,$occH) = @_;
83 :    
84 :     if (@{$node->[2]} > 1) # if not child
85 :     {
86 :     my $i;
87 :     my @current_values;
88 :     for ($i=1; ($i < @{$node->[2]}); $i++)
89 :     {
90 :     &pass1($node->[2]->[$i],$occH);
91 :     push(@current_values,$occH->{$node->[2]->[$i]->[0]});
92 :     }
93 :     $occH->{$node->[0]} = &iu(\@current_values);
94 :     }
95 :     }
96 :    
97 :     sub pass2 {
98 :     my($node,$occH) = @_;
99 :    
100 :     if (@{$node->[2]} > 1) # if not child
101 :     {
102 :     my $i;
103 :     for ($i=1; ($i < @{$node->[2]}); $i++)
104 :     {
105 : overbeek 1.2 my $overlap = &i([$occH->{$node->[0]},$occH->{$node->[2]->[$i]->[0]}]);
106 : overbeek 1.1 if (@$overlap > 0)
107 :     {
108 :     $occH->{$node->[2]->[$i]->[0]} = $overlap;
109 :     }
110 :     &pass2($node->[2]->[$i],$occH);
111 :     }
112 :     }
113 :     }
114 :     sub iu {
115 :     my($values) = @_;
116 :    
117 :     my $x = &i($values);
118 :     if (@$x > 0)
119 :     {
120 :     return $x;
121 :     }
122 :     return &u($values);
123 :     }
124 :    
125 :     sub i {
126 :     my($values) = @_;
127 :     my %counts;
128 :     # print &Dumper($values);
129 :     foreach my $v (@$values)
130 :     {
131 :     if (! ref($v)) { confess('BAD') }
132 :     foreach my $x (@$v)
133 :     {
134 :     $counts{$x}++;
135 :     }
136 :     }
137 : overbeek 1.2 my @tmp = sort grep { $counts{$_} == @$values } keys(%counts);
138 : overbeek 1.1 return \@tmp;
139 :     }
140 :    
141 :     sub u {
142 :     my($values) = @_;
143 :    
144 :     my %counts;
145 :     foreach my $v (@$values)
146 :     {
147 :     foreach my $x (@$v)
148 :     {
149 :     $counts{$x}++;
150 :     }
151 :     }
152 : overbeek 1.2 return [sort keys(%counts)];
153 : overbeek 1.1 }
154 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3