1 package Text::Tradition::Stemma;
4 use Encode qw( decode_utf8 );
9 use Graph::Reader::Dot;
10 use IPC::Run qw/ run binary /;
15 isa => 'Text::Tradition::Collation',
23 predicate => 'has_graph',
26 has distance_trees => (
28 isa => 'ArrayRef[Graph]',
29 writer => '_save_distance_trees',
30 predicate => 'has_distance_trees',
34 my( $self, $args ) = @_;
35 # If we have been handed a dotfile, initialize it into a graph.
36 if( exists $args->{'dot'} ) {
37 $self->graph_from_dot( $args->{'dot'} );
42 my( $self, $dotfh ) = @_;
44 binmode( $dotfh, ':utf8' );
45 my $reader = Graph::Reader::Dot->new();
46 my $graph = $reader->read_graph( $dotfh );
48 $self->graph( $graph );
49 # Go through the nodes and set any non-hypothetical node to extant.
50 foreach my $v ( $self->graph->vertices ) {
51 $self->graph->set_vertex_attribute( $v, 'class', 'extant' )
52 unless $self->graph->has_vertex_attribute( $v, 'class' );
55 warn "Failed to parse dot in $dotfh";
60 my( $self, $opts ) = @_;
62 # Get default and specified options
66 'hshape' => 'plaintext', # Shape for the hypothetical nodes
69 'fillcolor' => 'white',
70 'shape' => 'ellipse', # Shape for the extant nodes
73 'arrowhead' => 'open',
75 @graphopts{ keys %{$opts->{'graph'}} } = values %{$opts->{'graph'}}
77 @nodeopts{ keys %{$opts->{'node'}} } = values %{$opts->{'node'}}
79 @edgeopts{ keys %{$opts->{'edge'}} } = values %{$opts->{'edge'}}
83 push( @dotlines, 'digraph stemma {' );
84 ## Print out the global attributes
85 push( @dotlines, _make_dotline( 'graph', %graphopts ) ) if keys %graphopts;
86 push( @dotlines, _make_dotline( 'edge', %edgeopts ) ) if keys %edgeopts;
87 ## Delete our special attributes from the node set before continuing
88 my $hshape = delete $nodeopts{'hshape'};
89 my $htext = delete $nodeopts{'htext'};
90 push( @dotlines, _make_dotline( 'node', %nodeopts ) ) if keys %nodeopts;
92 # Add each of the nodes.
93 foreach my $n ( $self->graph->vertices ) {
94 if( $self->graph->get_vertex_attribute( $n, 'class' ) eq 'hypothetical' ) {
95 # Apply our display settings for hypothetical nodes.
96 push( @dotlines, _make_dotline( $n, 'shape' => $hshape, 'label' => $htext ) );
98 # Use the default display settings.
99 push( @dotlines, " $n;" );
102 # Add each of our edges.
103 foreach my $e ( $self->graph->edges ) {
104 my( $from, $to ) = @$e;
105 push( @dotlines, " $from -> $to;" );
107 push( @dotlines, '}' );
109 return join( "\n", @dotlines );
113 # Another version of dot output meant for graph editing, thus
118 push( @dotlines, 'digraph stemma {' );
119 my @real; # A cheap sort
120 foreach my $n ( sort $self->graph->vertices ) {
121 my $c = $self->graph->get_vertex_attribute( $n, 'class' );
122 $c = 'extant' unless $c;
123 if( $c eq 'extant' ) {
126 push( @dotlines, _make_dotline( $n, 'class' => $c ) );
129 # Now do the real ones
130 foreach my $n ( @real ) {
131 push( @dotlines, _make_dotline( $n, 'class' => 'extant' ) );
133 foreach my $e ( sort _by_vertex $self->graph->edges ) {
134 my( $from, $to ) = @$e;
135 push( @dotlines, " $from -> $to;" );
137 push( @dotlines, '}' );
138 return join( "\n", @dotlines );
142 my( $obj, %attr ) = @_;
144 foreach my $k ( keys %attr ) {
147 push( @pairs, "$k=\"$v\"" );
149 return sprintf( " %s [ %s ];", $obj, join( ', ', @pairs ) );
153 return $a->[0].$a->[1] cmp $b->[0].$b->[1];
156 # Render the stemma as SVG.
158 my( $self, $opts ) = @_;
159 my $dot = $self->as_dot( $opts );
160 my @cmd = qw/dot -Tsvg/;
162 my $dotfile = File::Temp->new();
164 # $dotfile->unlink_on_destroy(0);
165 binmode $dotfile, ':utf8';
167 push( @cmd, $dotfile->filename );
168 run( \@cmd, ">", binary(), \$svg );
169 $svg = decode_utf8( $svg );
175 my @wits = grep { $self->graph->get_vertex_attribute( $_, 'class' ) eq 'extant' }
176 $self->graph->vertices;
180 #### Methods for calculating phylogenetic trees ####
182 before 'distance_trees' => sub {
185 # TODO allow specification of method for calculating distance tree
186 if( $args{'recalc'} || !$self->has_distance_trees ) {
187 # We need to make a tree before we can return it.
188 my( $ok, $result ) = $self->run_phylip_pars();
190 # Save the resulting trees
191 my $trees = _parse_newick( $result );
192 $self->_save_distance_trees( $trees );
194 warn "Failed to calculate distance trees: $result";
199 sub make_character_matrix {
201 unless( $self->collation->linear ) {
202 warn "Need a linear graph in order to make an alignment table";
205 my $table = $self->collation->make_alignment_table;
206 # Push the names of the witnesses to initialize the rows of the matrix.
207 my @matrix = map { [ $self->_normalize_ac( $_ ) ] } @{$table->[0]};
208 foreach my $token_index ( 1 .. $#{$table} ) {
209 # First implementation: make dumb alignment table, caring about
210 # nothing except which reading is in which position.
211 my @chars = convert_characters( $table->[$token_index] );
212 foreach my $idx ( 0 .. $#matrix ) {
213 push( @{$matrix[$idx]}, $chars[$idx] );
220 my( $self, $witname ) = @_;
221 my $ac = $self->collation->ac_label;
222 if( $witname =~ /(.*)\Q$ac\E$/ ) {
223 $witname = $1 . '_ac';
225 return sprintf( "%-10s", $witname );
228 sub convert_characters {
230 # This is a simple algorithm that treats every reading as different.
231 # Eventually we will want to be able to specify how relationships
232 # affect the character matrix.
233 my %unique = ( '__UNDEF__' => 'X',
238 foreach my $word ( @$row ) {
239 if( $word && !exists $unique{$word} ) {
240 $unique{$word} = chr( 65 + $ctr );
243 $count{$word}++ if $word;
245 # Try to keep variants under 8 by lacunizing any singletons.
246 if( scalar( keys %unique ) > 8 ) {
247 foreach my $word ( keys %count ) {
248 if( $count{$word} == 1 ) {
249 $unique{$word} = '?';
253 my %u = reverse %unique;
254 if( scalar( keys %u ) > 8 ) {
255 warn "Have more than 8 variants on this location; phylip will break";
257 my @chars = map { $_ ? $unique{$_} : $unique{'__UNDEF__' } } @$row;
261 sub phylip_pars_input {
263 my $character_matrix = $self->make_character_matrix;
265 my $rows = scalar @{$character_matrix};
266 my $columns = scalar @{$character_matrix->[0]} - 1;
267 $input .= "\t$rows\t$columns\n";
268 foreach my $row ( @{$character_matrix} ) {
269 $input .= join( '', @$row ) . "\n";
274 sub run_phylip_pars {
277 # Set up a temporary directory for all the default Phylip files.
278 my $phylip_dir = File::Temp->newdir();
279 # $phylip_dir->unlink_on_destroy(0);
280 # We need an infile, and we need a command input file.
281 open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
282 print MATRIX $self->phylip_pars_input();
285 open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
286 ## TODO any configuration parameters we want to set here
287 # U Search for best tree? Yes
288 # S Search option? More thorough search
289 # V Number of trees to save? 100
290 # J Randomize input order of species? No. Use input order
291 # O Outgroup root? No, use as outgroup species 1
292 # T Use Threshold parsimony? No, use ordinary parsimony
293 # W Sites weighted? No
294 # M Analyze multiple data sets? No
295 # I Input species interleaved? Yes
296 # 0 Terminal type (IBM PC, ANSI, none)? ANSI
297 # 1 Print out the data at start of run No
298 # 2 Print indications of progress of run Yes
299 # 3 Print out tree Yes
300 # 4 Print out steps in each site No
301 # 5 Print character at all nodes of tree No
302 # 6 Write out trees onto tree file? Yes
306 # And then we run the program.
307 my $program = File::Which::which( 'pars' );
308 unless( -x $program ) {
309 return( undef, "Phylip pars not found in path" );
313 # We need to run it in our temporary directory where we have created
314 # all the expected files.
315 local $CWD = $phylip_dir;
316 my @cmd = ( $program );
317 run \@cmd, '<', 'cmdfile', '>', '/dev/null';
319 # Now our output should be in 'outfile' and our tree in 'outtree',
320 # both in the temp directory.
323 if( -f "$phylip_dir/outtree" ) {
324 open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
328 return( 1, join( '', @outtree ) ) if @outtree;
331 if( -f "$phylip_dir/outfile" ) {
332 open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
336 push( @error, "Neither outtree nor output file was produced!" );
338 return( undef, join( '', @error ) );
344 # Parse the result into a tree
345 my $forest = Bio::Phylo::IO->parse(
349 # Turn the tree into a graph, starting with the root node
350 foreach my $tree ( @{$forest->get_entities} ) {
351 push( @trees, _graph_from_bio( $tree ) );
356 sub _graph_from_bio {
358 my $graph = Graph->new( 'undirected' => 1 );
359 # Give all the intermediate anonymous nodes a name.
361 foreach my $n ( @{$tree->get_entities} ) {
362 next if $n->get_name;
363 $n->set_name( $i++ );
365 my $root = $tree->get_root->get_name;
366 $graph->add_vertex( $root );
367 _add_tree_children( $graph, $root, $tree->get_root->get_children() );
371 sub _add_tree_children {
372 my( $graph, $parent, $tree_children ) = @_;
373 foreach my $c ( @$tree_children ) {
374 my $child = $c->get_name;
375 $graph->add_vertex( $child );
376 $graph->add_path( $parent, $child );
377 _add_tree_children( $graph, $child, $c->get_children() );
382 __PACKAGE__->meta->make_immutable;