1 package Text::Tradition::Stemma;
4 use Encode qw( decode_utf8 );
9 use Graph::Reader::Dot;
10 use IPC::Run qw/ run binary /;
11 use Text::Tradition::StemmaUtil qw/ phylip_pars_input /;
16 isa => 'Text::Tradition::Collation',
24 predicate => 'has_graph',
27 has distance_trees => (
29 isa => 'ArrayRef[Graph]',
30 writer => '_save_distance_trees',
31 predicate => 'has_distance_trees',
34 has distance_program => (
41 my( $self, $args ) = @_;
42 # If we have been handed a dotfile, initialize it into a graph.
43 if( exists $args->{'dot'} ) {
44 $self->graph_from_dot( $args->{'dot'} );
49 my( $self, $dotfh ) = @_;
51 binmode( $dotfh, ':utf8' );
52 my $reader = Graph::Reader::Dot->new();
53 my $graph = $reader->read_graph( $dotfh );
55 $self->graph( $graph );
56 # Go through the nodes and set any non-hypothetical node to extant.
57 foreach my $v ( $self->graph->vertices ) {
58 $self->graph->set_vertex_attribute( $v, 'class', 'extant' )
59 unless $self->graph->has_vertex_attribute( $v, 'class' );
62 warn "Failed to parse dot in $dotfh";
67 my( $self, $opts ) = @_;
69 # Get default and specified options
73 'hshape' => 'plaintext', # Shape for the hypothetical nodes
76 'fillcolor' => 'white',
77 'shape' => 'ellipse', # Shape for the extant nodes
80 'arrowhead' => 'open',
82 @graphopts{ keys %{$opts->{'graph'}} } = values %{$opts->{'graph'}}
84 @nodeopts{ keys %{$opts->{'node'}} } = values %{$opts->{'node'}}
86 @edgeopts{ keys %{$opts->{'edge'}} } = values %{$opts->{'edge'}}
90 push( @dotlines, 'digraph stemma {' );
91 ## Print out the global attributes
92 push( @dotlines, _make_dotline( 'graph', %graphopts ) ) if keys %graphopts;
93 push( @dotlines, _make_dotline( 'edge', %edgeopts ) ) if keys %edgeopts;
94 ## Delete our special attributes from the node set before continuing
95 my $hshape = delete $nodeopts{'hshape'};
96 my $htext = delete $nodeopts{'htext'};
97 push( @dotlines, _make_dotline( 'node', %nodeopts ) ) if keys %nodeopts;
99 # Add each of the nodes.
100 foreach my $n ( $self->graph->vertices ) {
101 if( $self->graph->get_vertex_attribute( $n, 'class' ) eq 'hypothetical' ) {
102 # Apply our display settings for hypothetical nodes.
103 push( @dotlines, _make_dotline( $n, 'shape' => $hshape, 'label' => $htext ) );
105 # Use the default display settings.
106 push( @dotlines, " $n;" );
109 # Add each of our edges.
110 foreach my $e ( $self->graph->edges ) {
111 my( $from, $to ) = @$e;
112 push( @dotlines, " $from -> $to;" );
114 push( @dotlines, '}' );
116 return join( "\n", @dotlines );
120 # Another version of dot output meant for graph editing, thus
125 push( @dotlines, 'digraph stemma {' );
126 my @real; # A cheap sort
127 foreach my $n ( sort $self->graph->vertices ) {
128 my $c = $self->graph->get_vertex_attribute( $n, 'class' );
129 $c = 'extant' unless $c;
130 if( $c eq 'extant' ) {
133 push( @dotlines, _make_dotline( $n, 'class' => $c ) );
136 # Now do the real ones
137 foreach my $n ( @real ) {
138 push( @dotlines, _make_dotline( $n, 'class' => 'extant' ) );
140 foreach my $e ( sort _by_vertex $self->graph->edges ) {
141 my( $from, $to ) = @$e;
142 push( @dotlines, " $from -> $to;" );
144 push( @dotlines, '}' );
145 return join( "\n", @dotlines );
149 my( $obj, %attr ) = @_;
151 foreach my $k ( keys %attr ) {
154 push( @pairs, "$k=\"$v\"" );
156 return sprintf( " %s [ %s ];", $obj, join( ', ', @pairs ) );
160 return $a->[0].$a->[1] cmp $b->[0].$b->[1];
163 # Render the stemma as SVG.
165 my( $self, $opts ) = @_;
166 my $dot = $self->as_dot( $opts );
167 my @cmd = qw/dot -Tsvg/;
169 my $dotfile = File::Temp->new();
171 # $dotfile->unlink_on_destroy(0);
172 binmode $dotfile, ':utf8';
174 push( @cmd, $dotfile->filename );
175 run( \@cmd, ">", binary(), \$svg );
176 $svg = decode_utf8( $svg );
182 my @wits = grep { $self->graph->get_vertex_attribute( $_, 'class' ) eq 'extant' }
183 $self->graph->vertices;
187 #### Methods for calculating phylogenetic trees ####
189 before 'distance_trees' => sub {
192 'program' => 'phylip_pars',
194 # TODO allow specification of method for calculating distance tree
195 if( !$self->has_distance_trees
196 || $args{'program'} ne $self->distance_program ) {
197 # We need to make a tree before we can return it.
198 my $dsub = 'run_' . $args{'program'};
199 my( $ok, $result ) = $self->$dsub();
201 # Save the resulting trees
202 my $trees = _parse_newick( $result );
203 $self->_save_distance_trees( $trees );
204 $self->distance_program( $args{'program'} );
206 warn "Failed to calculate distance trees: $result";
211 sub run_phylip_pars {
214 # Set up a temporary directory for all the default Phylip files.
215 my $phylip_dir = File::Temp->newdir();
216 # $phylip_dir->unlink_on_destroy(0);
217 # We need an infile, and we need a command input file.
218 open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
219 print MATRIX phylip_pars_input( $self->collation->make_alignment_table() );
222 open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
223 ## TODO any configuration parameters we want to set here
224 # U Search for best tree? Yes
225 # S Search option? More thorough search
226 # V Number of trees to save? 100
227 # J Randomize input order of species? No. Use input order
228 # O Outgroup root? No, use as outgroup species 1
229 # T Use Threshold parsimony? No, use ordinary parsimony
230 # W Sites weighted? No
231 # M Analyze multiple data sets? No
232 # I Input species interleaved? Yes
233 # 0 Terminal type (IBM PC, ANSI, none)? ANSI
234 # 1 Print out the data at start of run No
235 # 2 Print indications of progress of run Yes
236 # 3 Print out tree Yes
237 # 4 Print out steps in each site No
238 # 5 Print character at all nodes of tree No
239 # 6 Write out trees onto tree file? Yes
243 # And then we run the program.
244 my $program = File::Which::which( 'pars' );
245 unless( -x $program ) {
246 return( undef, "Phylip pars not found in path" );
250 # We need to run it in our temporary directory where we have created
251 # all the expected files.
252 local $CWD = $phylip_dir;
253 my @cmd = ( $program );
254 run \@cmd, '<', 'cmdfile', '>', '/dev/null';
256 # Now our output should be in 'outfile' and our tree in 'outtree',
257 # both in the temp directory.
260 if( -f "$phylip_dir/outtree" ) {
261 open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
265 return( 1, join( '', @outtree ) ) if @outtree;
268 if( -f "$phylip_dir/outfile" ) {
269 open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
273 push( @error, "Neither outtree nor output file was produced!" );
275 return( undef, join( '', @error ) );
281 # Parse the result into a tree
282 my $forest = Bio::Phylo::IO->parse(
286 # Turn the tree into a graph, starting with the root node
287 foreach my $tree ( @{$forest->get_entities} ) {
288 push( @trees, _graph_from_bio( $tree ) );
293 sub _graph_from_bio {
295 my $graph = Graph->new( 'undirected' => 1 );
296 # Give all the intermediate anonymous nodes a name.
298 foreach my $n ( @{$tree->get_entities} ) {
299 next if $n->get_name;
300 $n->set_name( $i++ );
302 my $root = $tree->get_root->get_name;
303 $graph->add_vertex( $root );
304 _add_tree_children( $graph, $root, $tree->get_root->get_children() );
308 sub _add_tree_children {
309 my( $graph, $parent, $tree_children ) = @_;
310 foreach my $c ( @$tree_children ) {
311 my $child = $c->get_name;
312 $graph->add_vertex( $child );
313 $graph->add_path( $parent, $child );
314 _add_tree_children( $graph, $child, $c->get_children() );
319 __PACKAGE__->meta->make_immutable;