1 package Text::Tradition::Stemma;
7 use Graph::Reader::Dot;
8 use IPC::Run qw/ run /;
10 use Text::Balanced qw/ extract_bracketed /;
14 isa => 'Text::Tradition::Collation',
18 has character_matrix => (
20 isa => 'ArrayRef[ArrayRef[Str]]',
21 writer => '_save_character_matrix',
22 predicate => 'has_character_matrix',
28 predicate => 'has_graph',
36 has distance_trees => (
38 isa => 'ArrayRef[Graph]',
39 writer => '_save_distance_trees',
40 predicate => 'has_distance_trees',
44 my( $self, $args ) = @_;
45 # If we have been handed a dotfile, initialize it into a graph.
46 if( exists $args->{'dot'} ) {
47 my $reader = Graph::Reader::Dot->new();
48 my $graph = $reader->read_graph( $args->{'dot'} );
50 ? $self->graph( $graph )
51 : warn "Failed to parse dot file " . $args->{'dot'};
54 # If we have a graph, calculate all the shortest paths between nodes,
55 # disregarding direction.
56 if( $self->has_graph ) {
58 if( $self->graph->is_directed ) {
59 # Make an undirected version.
60 $undirected = Graph->new( 'undirected' => 1 );
61 foreach my $v ( $self->graph->vertices ) {
62 $undirected->add_vertex( $v );
64 foreach my $e ( $self->graph->edges ) {
65 $undirected->add_edge( @$e );
68 $undirected = $self->graph;
70 $self->apsp( $undirected->APSP_Floyd_Warshall() );
74 before 'distance_trees' => sub {
77 # TODO allow specification of method for calculating distance tree
78 if( $args{'recalc'} || !$self->has_distance_trees ) {
79 # We need to make a tree before we can return it.
80 my( $ok, $result ) = $self->run_phylip_pars();
82 $self->_save_distance_trees( _parse_newick( $result ) );
84 warn "Failed to calculate distance tree: $result";
89 sub make_character_matrix {
91 unless( $self->collation->linear ) {
92 warn "Need a linear graph in order to make an alignment table";
95 my $table = $self->collation->make_alignment_table;
96 # Push the names of the witnesses to initialize the rows of the matrix.
97 my @matrix = map { [ $self->_normalize_ac( $_ ) ] } @{$table->[0]};
99 foreach my $token_index ( 1 .. $#{$table} ) {
100 # First implementation: make dumb alignment table, caring about
101 # nothing except which reading is in which position.
102 my @chars = convert_characters( $table->[$token_index] );
103 foreach my $idx ( 0 .. $#matrix ) {
104 push( @{$matrix[$idx]}, $chars[$idx] );
107 $self->_save_character_matrix( \@matrix );
111 my( $self, $witname ) = @_;
112 my $ac = $self->collation->ac_label;
113 if( $witname =~ /(.*)\Q$ac\E$/ ) {
114 $witname = $1 . '_ac';
116 return sprintf( "%-10s", $witname );
119 sub convert_characters {
121 # This is a simple algorithm that treats every reading as different.
122 # Eventually we will want to be able to specify how relationships
123 # affect the character matrix.
124 my %unique = ( '__UNDEF__' => 'X',
128 foreach my $word ( @$row ) {
129 if( $word && !exists $unique{$word} ) {
130 $unique{$word} = chr( 65 + $ctr );
134 if( scalar( keys %unique ) > 8 ) {
135 warn "Have more than 8 variants on this location; phylip will break";
137 my @chars = map { $_ ? $unique{$_} : $unique{'__UNDEF__' } } @$row;
141 sub phylip_pars_input {
143 $self->make_character_matrix unless $self->has_character_matrix;
145 my $rows = scalar @{$self->character_matrix};
146 my $columns = scalar @{$self->character_matrix->[0]} - 1;
147 $matrix .= "\t$rows\t$columns\n";
148 foreach my $row ( @{$self->character_matrix} ) {
149 $matrix .= join( '', @$row ) . "\n";
154 sub run_phylip_pars {
157 # Set up a temporary directory for all the default Phylip files.
158 my $phylip_dir = File::Temp->newdir();
159 # $phylip_dir->unlink_on_destroy(0);
160 # We need an infile, and we need a command input file.
161 open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
162 print MATRIX $self->phylip_pars_input();
165 open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
166 ## TODO any configuration parameters we want to set here
167 # U Search for best tree? Yes
168 # S Search option? More thorough search
169 # V Number of trees to save? 100
170 # J Randomize input order of species? No. Use input order
171 # O Outgroup root? No, use as outgroup species 1
172 # T Use Threshold parsimony? No, use ordinary parsimony
173 # W Sites weighted? No
174 # M Analyze multiple data sets? No
175 # I Input species interleaved? Yes
176 # 0 Terminal type (IBM PC, ANSI, none)? ANSI
177 # 1 Print out the data at start of run No
178 # 2 Print indications of progress of run Yes
179 # 3 Print out tree Yes
180 # 4 Print out steps in each site No
181 # 5 Print character at all nodes of tree No
182 # 6 Write out trees onto tree file? Yes
186 # And then we run the program.
188 my $PHYLIP_PATH = '/Users/tla/Projects/phylip-3.69/exe';
189 my $program = "pars";
190 if( $^O eq 'darwin' ) {
191 $program = "$PHYLIP_PATH/$program.app/Contents/MacOS/$program";
193 $program = "$PHYLIP_PATH/$program";
197 # We need to run it in our temporary directory where we have created
198 # all the expected files.
199 local $CWD = $phylip_dir;
200 my @cmd = ( $program );
201 run \@cmd, '<', 'cmdfile', '>', '/dev/null';
203 # Now our output should be in 'outfile' and our tree in 'outtree',
204 # both in the temp directory.
207 if( -f "$phylip_dir/outtree" ) {
208 open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
212 return( 1, join( '', @outtree ) ) if @outtree;
215 if( -f "$phylip_dir/outfile" ) {
216 open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
220 push( @error, "Neither outtree nor output file was produced!" );
222 return( undef, join( '', @error ) );
228 # Parse the result into a tree
229 my $forest = Bio::Phylo::IO->parse(
233 # Turn the tree into a graph, starting with the root node
234 foreach my $tree ( @{$forest->get_entities} ) {
235 push( @trees, _graph_from_bio( $tree ) );
240 sub _graph_from_bio {
242 my $graph = Graph->new( 'undirected' => 1 );
243 # Give all the intermediate anonymous nodes a name.
245 foreach my $n ( @{$tree->get_entities} ) {
246 next if $n->get_name;
247 $n->set_name( $i++ );
249 my $root = $tree->get_root->get_name;
250 $graph->add_vertex( $root );
251 _add_tree_children( $graph, $root, $tree->get_root->get_children() );
255 sub _add_tree_children {
256 my( $graph, $parent, $tree_children ) = @_;
257 foreach my $c ( @$tree_children ) {
258 my $child = $c->get_name;
259 $graph->add_vertex( $child );
260 $graph->add_path( $parent, $child );
261 _add_tree_children( $graph, $child, $c->get_children() );
266 __PACKAGE__->meta->make_immutable;