1 package Text::Tradition::Stemma;
4 use Encode qw( decode_utf8 );
9 use Graph::Reader::Dot;
10 use IPC::Run qw/ run binary /;
12 use Text::Balanced qw/ extract_bracketed /;
16 isa => 'Text::Tradition::Collation',
20 has character_matrix => (
22 isa => 'ArrayRef[ArrayRef[Str]]',
23 writer => '_save_character_matrix',
24 predicate => 'has_character_matrix',
30 predicate => 'has_graph',
38 has distance_trees => (
40 isa => 'ArrayRef[Graph]',
41 writer => '_save_distance_trees',
42 predicate => 'has_distance_trees',
46 my( $self, $args ) = @_;
47 # If we have been handed a dotfile, initialize it into a graph.
48 if( exists $args->{'dot'} ) {
49 # Open the file, assume UTF-8
50 open( my $dot, $args->{'dot'} ) or warn "Failed to read dot file";
51 # TODO don't bother if we haven't opened
52 binmode $dot, ":utf8";
53 my $reader = Graph::Reader::Dot->new();
54 my $graph = $reader->read_graph( $dot );
56 ? $self->graph( $graph )
57 : warn "Failed to parse dot file " . $args->{'dot'};
60 # If we have a graph, calculate all the shortest paths between nodes,
61 # disregarding direction.
62 if( $self->has_graph ) {
64 if( $self->graph->is_directed ) {
65 # Make an undirected version.
66 $undirected = Graph->new( 'undirected' => 1 );
67 foreach my $v ( $self->graph->vertices ) {
68 $undirected->add_vertex( $v );
70 foreach my $e ( $self->graph->edges ) {
71 $undirected->add_edge( @$e );
74 $undirected = $self->graph;
76 $self->apsp( $undirected->APSP_Floyd_Warshall() );
80 # Render the stemma as SVG.
83 # TODO add options for display, someday
84 my $dgraph = Graph::Convert->as_graph_easy( $self->graph );
85 # Set some class display attributes for 'hypothetical' and 'extant' nodes
86 $dgraph->set_attribute( 'flow', 'south' );
87 foreach my $n ( $dgraph->nodes ) {
88 if( $n->attribute( 'class' ) eq 'hypothetical' ) {
89 $n->set_attribute( 'shape', 'point' );
90 $n->set_attribute( 'pointshape', 'diamond' );
92 $n->set_attribute( 'shape', 'ellipse' );
96 # Render to svg via graphviz
97 my @cmd = qw/dot -Tsvg/;
99 my $dotfile = File::Temp->new();
101 # $dotfile->unlink_on_destroy(0);
102 binmode $dotfile, ':utf8';
103 print $dotfile $dgraph->as_graphviz();
104 push( @cmd, $dotfile->filename );
105 run( \@cmd, ">", binary(), \$svg );
106 $svg = decode_utf8( $svg );
110 #### Methods for calculating phylogenetic trees ####
112 before 'distance_trees' => sub {
115 # TODO allow specification of method for calculating distance tree
116 if( $args{'recalc'} || !$self->has_distance_trees ) {
117 # We need to make a tree before we can return it.
118 my( $ok, $result ) = $self->run_phylip_pars();
120 $self->_save_distance_trees( _parse_newick( $result ) );
122 warn "Failed to calculate distance tree: $result";
127 sub make_character_matrix {
129 unless( $self->collation->linear ) {
130 warn "Need a linear graph in order to make an alignment table";
133 my $table = $self->collation->make_alignment_table;
134 # Push the names of the witnesses to initialize the rows of the matrix.
135 my @matrix = map { [ $self->_normalize_ac( $_ ) ] } @{$table->[0]};
137 foreach my $token_index ( 1 .. $#{$table} ) {
138 # First implementation: make dumb alignment table, caring about
139 # nothing except which reading is in which position.
140 my @chars = convert_characters( $table->[$token_index] );
141 foreach my $idx ( 0 .. $#matrix ) {
142 push( @{$matrix[$idx]}, $chars[$idx] );
145 $self->_save_character_matrix( \@matrix );
149 my( $self, $witname ) = @_;
150 my $ac = $self->collation->ac_label;
151 if( $witname =~ /(.*)\Q$ac\E$/ ) {
152 $witname = $1 . '_ac';
154 return sprintf( "%-10s", $witname );
157 sub convert_characters {
159 # This is a simple algorithm that treats every reading as different.
160 # Eventually we will want to be able to specify how relationships
161 # affect the character matrix.
162 my %unique = ( '__UNDEF__' => 'X',
166 foreach my $word ( @$row ) {
167 if( $word && !exists $unique{$word} ) {
168 $unique{$word} = chr( 65 + $ctr );
172 if( scalar( keys %unique ) > 8 ) {
173 warn "Have more than 8 variants on this location; phylip will break";
175 my @chars = map { $_ ? $unique{$_} : $unique{'__UNDEF__' } } @$row;
179 sub phylip_pars_input {
181 $self->make_character_matrix unless $self->has_character_matrix;
183 my $rows = scalar @{$self->character_matrix};
184 my $columns = scalar @{$self->character_matrix->[0]} - 1;
185 $matrix .= "\t$rows\t$columns\n";
186 foreach my $row ( @{$self->character_matrix} ) {
187 $matrix .= join( '', @$row ) . "\n";
192 sub run_phylip_pars {
195 # Set up a temporary directory for all the default Phylip files.
196 my $phylip_dir = File::Temp->newdir();
197 # $phylip_dir->unlink_on_destroy(0);
198 # We need an infile, and we need a command input file.
199 open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
200 print MATRIX $self->phylip_pars_input();
203 open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
204 ## TODO any configuration parameters we want to set here
205 # U Search for best tree? Yes
206 # S Search option? More thorough search
207 # V Number of trees to save? 100
208 # J Randomize input order of species? No. Use input order
209 # O Outgroup root? No, use as outgroup species 1
210 # T Use Threshold parsimony? No, use ordinary parsimony
211 # W Sites weighted? No
212 # M Analyze multiple data sets? No
213 # I Input species interleaved? Yes
214 # 0 Terminal type (IBM PC, ANSI, none)? ANSI
215 # 1 Print out the data at start of run No
216 # 2 Print indications of progress of run Yes
217 # 3 Print out tree Yes
218 # 4 Print out steps in each site No
219 # 5 Print character at all nodes of tree No
220 # 6 Write out trees onto tree file? Yes
224 # And then we run the program.
226 my $PHYLIP_PATH = '/Users/tla/Projects/phylip-3.69/exe';
227 my $program = "pars";
228 if( $^O eq 'darwin' ) {
229 $program = "$PHYLIP_PATH/$program.app/Contents/MacOS/$program";
231 $program = "$PHYLIP_PATH/$program";
235 # We need to run it in our temporary directory where we have created
236 # all the expected files.
237 local $CWD = $phylip_dir;
238 my @cmd = ( $program );
239 run \@cmd, '<', 'cmdfile', '>', '/dev/null';
241 # Now our output should be in 'outfile' and our tree in 'outtree',
242 # both in the temp directory.
245 if( -f "$phylip_dir/outtree" ) {
246 open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
250 return( 1, join( '', @outtree ) ) if @outtree;
253 if( -f "$phylip_dir/outfile" ) {
254 open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
258 push( @error, "Neither outtree nor output file was produced!" );
260 return( undef, join( '', @error ) );
266 # Parse the result into a tree
267 my $forest = Bio::Phylo::IO->parse(
271 # Turn the tree into a graph, starting with the root node
272 foreach my $tree ( @{$forest->get_entities} ) {
273 push( @trees, _graph_from_bio( $tree ) );
278 sub _graph_from_bio {
280 my $graph = Graph->new( 'undirected' => 1 );
281 # Give all the intermediate anonymous nodes a name.
283 foreach my $n ( @{$tree->get_entities} ) {
284 next if $n->get_name;
285 $n->set_name( $i++ );
287 my $root = $tree->get_root->get_name;
288 $graph->add_vertex( $root );
289 _add_tree_children( $graph, $root, $tree->get_root->get_children() );
293 sub _add_tree_children {
294 my( $graph, $parent, $tree_children ) = @_;
295 foreach my $c ( @$tree_children ) {
296 my $child = $c->get_name;
297 $graph->add_vertex( $child );
298 $graph->add_path( $parent, $child );
299 _add_tree_children( $graph, $child, $c->get_children() );
304 __PACKAGE__->meta->make_immutable;