analysis/lib/Text/Tradition/StemmaUtil.pm

   1 package Text::Tradition::StemmaUtil;
   2
   3 use strict;
   4 use warnings;
   5 use Exporter 'import';
   6 use vars qw/ @EXPORT_OK /;
   7 use Bio::Phylo::IO;
   8 use Encode qw( decode_utf8 );
   9 use File::chdir;
  10 use File::Temp;
  11 use File::Which;
  12 use Graph;
  13 use Graph::Reader::Dot;
  14 use IPC::Run qw/ run binary /;
  15 use Text::Tradition::Error;
  16 @EXPORT_OK = qw/ read_graph display_graph editable_graph
  17         character_input phylip_pars parse_newick newick_to_svg /;
  18
  19 =head1 NAME
  20
  21 Text::Tradition::StemmaUtil - standalone utilities for stemma graph display and
  22 distance tree calculations
  23
  24 =head1 DESCRIPTION
  25
  26 This package contains a set of utilities for displaying arbitrary stemmata and
  27 running phylogenetic analysis on text collations.
  28
  29 =head1 SUBROUTINES
  30
  31 =head2 read_graph( $dotstr) {
  32
  33 Parses the graph specification on the filehandle in $dotstr and returns a Graph
  34 object. This subroutine works around some deficiencies in Graph::Reader::Dot.
  35
  36 =cut
  37
  38 sub read_graph {
  39         my $dotstr = shift;
  40         # Graph::Reader::Dot does not handle bare non-ASCII Unicode word characters.
  41         # We get around this by wrapping all words in double quotes, as long as they
  42         # aren't already wrapped, and as long as they aren't the initial '(di)?graph'.
  43         # Also need to deal correctly with the graph title.
  44         if( $dotstr =~ /^\s*((di)?graph)\s+(.*?)\s*\{(.*)$/s ) {
  45                 my( $decl, $ident, $rest ) = ( $1, $3, $4 );
  46                 unless( substr( $ident, 0, 1 ) eq '"' ) {
  47                         $ident = '"'.$ident.'"';
  48                 }
  49                 $rest =~ s/(?<!")\b(\w+)\b(?!")/"$1"/g;
  50                 $dotstr = "$decl $ident { $rest";
  51         }
  52
  53         # Now open a filehandle onto the string and pass it to Graph::Reader::Dot.
  54         my $dotfh;
  55         open $dotfh, '<', \$dotstr;
  56         binmode $dotfh, ':utf8';
  57         my $reader = Graph::Reader::Dot->new();
  58         # Redirect STDOUT in order to trap any error messages - syntax errors
  59         # are evidently not fatal.
  60         my $graph;
  61         my $reader_out;
  62         my $reader_err;
  63         {
  64                 local(*STDOUT);
  65                 open( STDOUT, ">", \$reader_out );
  66                 local(*STDERR);
  67                 open( STDERR, ">", \$reader_err );
  68                 $graph = $reader->read_graph( $dotfh );
  69                 close STDOUT;
  70                 close STDERR;
  71         }
  72         if( $reader_out && $reader_out =~ /error/s ) {
  73                 throw( "Error trying to parse dot: $reader_out" );
  74         } elsif( !$graph ) {
  75                 throw( "Failed to create graph from dot" );
  76         }
  77         # Wrench the graph identifier out of the graph
  78         ## HORRIBLE HACK but there is no API access to the graph identifier!
  79         $graph->set_graph_attribute( 'name', $graph->[4]->{'name'} );
  80
  81         # Correct for implicit graph -> digraph quirk of reader
  82         if( $reader_err && $reader_err =~ /graph will be treated as digraph/ ) {
  83                 my $udgraph = $graph->undirected_copy;
  84                 $udgraph->set_graph_attributes( $graph->get_graph_attributes );
  85                 foreach my $v ( $graph->vertices ) {
  86                         $udgraph->set_vertex_attributes( $v, $graph->get_vertex_attributes( $v ) );
  87                 }
  88                 $graph = $udgraph;
  89         }
  90
  91         return $graph;
  92 }
  93
  94 =head2 display_graph( $graph, $opts )
  95
  96 Returns a dot specification intended for display, according to the logical
  97 attributes of the witnesses.
  98
  99 =cut
 100
 101 sub display_graph {
 102     my( $graph, $opts ) = @_;
 103
 104     # Get default and specified options
 105     my %graphopts = (
 106         # 'ratio' => 1,
 107         'bgcolor' => 'transparent',
 108     );
 109     my %nodeopts = (
 110                 'fontsize' => 11,
 111                 'style' => 'filled',
 112                 'fillcolor' => 'white',
 113                 'color' => 'white',
 114                 'shape' => 'ellipse',   # Shape for the extant nodes
 115         );
 116         my %edgeopts = (
 117                 'arrowhead' => 'none',
 118         );
 119         @graphopts{ keys %{$opts->{'graph'}} } = values %{$opts->{'graph'}}
 120                 if $opts->{'graph'};
 121         @nodeopts{ keys %{$opts->{'node'}} } = values %{$opts->{'node'}}
 122                 if $opts->{'node'};
 123         @edgeopts{ keys %{$opts->{'edge'}} } = values %{$opts->{'edge'}}
 124                 if $opts->{'edge'};
 125
 126         my $gdecl = $graph->is_directed ? 'digraph' : 'graph';
 127         my $gname = $opts->{'name'} ? '"' . $opts->{'name'} . '"'
 128                 : 'stemma';
 129         my @dotlines;
 130         push( @dotlines, "$gdecl $gname {" );
 131         ## Print out the global attributes
 132         push( @dotlines, _make_dotline( 'graph', %graphopts ) ) if keys %graphopts;
 133         push( @dotlines, _make_dotline( 'edge', %edgeopts ) ) if keys %edgeopts;
 134         push( @dotlines, _make_dotline( 'node', %nodeopts ) ) if keys %nodeopts;
 135
 136         # Add each of the nodes.
 137     foreach my $n ( $graph->vertices ) {
 138         my %vattr = ( 'id' => $n );  # Set the SVG element ID to the sigil itself
 139         if( $graph->has_vertex_attribute( $n, 'label' ) ) {
 140                 $vattr{'label'} = $graph->get_vertex_attribute( $n, 'label' );
 141         }
 142                 push( @dotlines, _make_dotline( $n, %vattr ) );
 143     }
 144     # Add each of our edges.
 145     foreach my $e ( $graph->edges ) {
 146         my( $from, $to ) = map { _dotquote( $_ ) } @$e;
 147         my $connector = $graph->is_directed ? '->' : '--';
 148         push( @dotlines, "  $from $connector $to;" );
 149     }
 150     push( @dotlines, '}' );
 151
 152     return join( "\n", @dotlines );
 153 }
 154
 155
 156 =head2 editable_graph( $graph, $opts )
 157
 158 Returns a dot specification of a stemma graph with logical witness features,
 159 intended for editing the stemma definition.
 160
 161 =cut
 162
 163 sub editable_graph {
 164         my( $graph, $opts ) = @_;
 165
 166         # Create the graph
 167         my $join = ( $opts && exists $opts->{'linesep'} ) ? $opts->{'linesep'} : "\n";
 168         my $gdecl = $graph->is_undirected ? 'graph' : 'digraph';
 169         my $gname = exists $opts->{'name'} ? '"' . $opts->{'name'} . '"'
 170                 : 'stemma';
 171         my @dotlines;
 172         push( @dotlines, "$gdecl $gname {" );
 173         my @real; # A cheap sort
 174     foreach my $n ( sort $graph->vertices ) {
 175         my $c = $graph->get_vertex_attribute( $n, 'class' );
 176         $c = 'extant' unless $c;
 177         if( $c eq 'extant' ) {
 178                 push( @real, $n );
 179         } else {
 180                         push( @dotlines, _make_dotline( $n, 'class' => $c ) );
 181                 }
 182     }
 183         # Now do the real ones
 184         foreach my $n ( @real ) {
 185                 push( @dotlines, _make_dotline( $n, 'class' => 'extant' ) );
 186         }
 187         foreach my $e ( sort _by_vertex $graph->edges ) {
 188                 my( $from, $to ) = map { _dotquote( $_ ) } @$e;
 189                 my $conn = $graph->is_undirected ? '--' : '->';
 190                 push( @dotlines, "  $from $conn $to;" );
 191         }
 192     push( @dotlines, '}' );
 193     return join( $join, @dotlines );
 194 }
 195
 196 sub _make_dotline {
 197         my( $obj, %attr ) = @_;
 198         my @pairs;
 199         foreach my $k ( keys %attr ) {
 200                 my $v = _dotquote( $attr{$k} );
 201                 push( @pairs, "$k=$v" );
 202         }
 203         return sprintf( "  %s [ %s ];", _dotquote( $obj ), join( ', ', @pairs ) );
 204 }
 205
 206 sub _dotquote {
 207         my( $str ) = @_;
 208         return $str if $str =~ /^[A-Za-z0-9_]+$/;
 209         $str =~ s/\"/\\\"/g;
 210         $str = '"' . $str . '"';
 211         return $str;
 212 }
 213
 214 sub _by_vertex {
 215         return $a->[0].$a->[1] cmp $b->[0].$b->[1];
 216 }
 217
 218 =head2 character_input( $tradition, $opts )
 219
 220 Returns a character matrix string suitable for Phylip programs, which
 221 corresponds to the given alignment table.  See Text::Tradition::Collation
 222 for a description of the alignment table format. Options include:
 223
 224 =over
 225
 226 =item * exclude_layer - Exclude layered witnesses from the character input,
 227 using only the 'main' text of the witnesses in the tradition.
 228
 229 =item * collapse - A reference to an array of relationship names that should
 230 be treated as equivalent for the purposes of generating the character matrix.
 231
 232 =back
 233
 234 =cut
 235
 236 sub character_input {
 237     my ( $tradition, $opts ) = @_;
 238     my $table = $tradition->collation->alignment_table;
 239     if( $opts->{exclude_layer} ) {
 240         # Filter out all alignment table rows that do not correspond
 241         # to a named witness - these are the layered witnesses.
 242         my $newtable = { alignment => [], length => $table->{length} };
 243         foreach my $row ( @{$table->{alignment}} ) {
 244                 if( $tradition->has_witness( $row->{witness} ) ) {
 245                         push( @{$newtable->{alignment}}, $row );
 246                 }
 247         }
 248         $table = $newtable;
 249     }
 250     my $character_matrix = _make_character_matrix( $table, $opts );
 251     my $input = '';
 252     my $rows = scalar @{$character_matrix};
 253     my $columns = scalar @{$character_matrix->[0]} - 1;
 254     $input .= "\t$rows\t$columns\n";
 255     foreach my $row ( @{$character_matrix} ) {
 256         $input .= join( '', @$row ) . "\n";
 257     }
 258     return $input;
 259 }
 260
 261 sub _make_character_matrix {
 262     my( $table, $opts ) = @_;
 263     # Push the names of the witnesses to initialize the rows of the matrix.
 264     my @matrix = map { [ _normalize_witname( $_->{'witness'} ) ] }
 265                                 @{$table->{'alignment'}};
 266     foreach my $token_index ( 0 .. $table->{'length'} - 1) {
 267         my @pos_tokens = map { $_->{'tokens'}->[$token_index] }
 268                                                         @{$table->{'alignment'}};
 269         my @pos_readings = map { $_ ? $_->{'t'} : $_ } @pos_tokens;
 270         my @chars = _convert_characters( \@pos_readings, $opts );
 271         foreach my $idx ( 0 .. $#matrix ) {
 272             push( @{$matrix[$idx]}, $chars[$idx] );
 273         }
 274     }
 275     return \@matrix;
 276 }
 277
 278 # Helper function to make the witness name something legal for pars
 279
 280 sub _normalize_witname {
 281     my( $witname ) = @_;
 282     $witname =~ s/\s+/ /g;
 283     $witname =~ s/[\[\]\(\)\:;,]//g;
 284     $witname = substr( $witname, 0, 10 );
 285     return sprintf( "%-10s", $witname );
 286 }
 287
 288 sub _convert_characters {
 289     my( $row, $opts ) = @_;
 290     # This is a simple algorithm that treats every reading as different.
 291     # Eventually we will want to be able to specify how relationships
 292     # affect the character matrix.
 293     my %unique = ( '__UNDEF__' => 'X',
 294                    '#LACUNA#'  => '?',
 295                  );
 296     my %equivalent;
 297     my %count;
 298     my $ctr = 0;
 299     foreach my $rdg ( @$row ) {
 300         next unless $rdg;
 301         next if $rdg->is_lacuna;
 302                 next if exists $unique{$rdg->text};
 303                 if( ref( $opts->{'collapse'} ) eq 'ARRAY' ) {
 304                         my @exclude_types = @{$opts->{'collapse'}};
 305                         my @set = $rdg->related_readings( sub { my $rel = shift;
 306                                 $rel->colocated && grep { $rel->type eq $_ } @exclude_types } );
 307                         push( @set, $rdg );
 308                         my $char = chr( 65 + $ctr++ );
 309                         map { $unique{$_->text} = $char } @set;
 310                         $count{$rdg->text} += scalar @set;
 311                 } else {
 312                         $unique{$rdg->text} = chr( 65 + $ctr++ );
 313                         $count{$rdg->text}++;
 314                 }
 315     }
 316     # Try to keep variants under 8 by lacunizing any singletons.
 317     if( scalar( keys %unique ) > 8 ) {
 318                 foreach my $word ( keys %count ) {
 319                         if( $count{$word} == 1 ) {
 320                                 $unique{$word} = '?';
 321                         }
 322                 }
 323     }
 324     my %u = reverse %unique;
 325     if( scalar( keys %u ) > 8 ) {
 326         warn "Have more than 8 variants on this location; phylip will break";
 327     }
 328     my @chars = map { $_ ? $unique{$_->text} : $unique{'__UNDEF__' } } @$row;
 329     return @chars;
 330 }
 331
 332 =head2 phylip_pars( $character_matrix )
 333
 334 Runs Phylip Pars on the given character matrix.  Returns results in Newick format.
 335
 336 =cut
 337
 338 sub phylip_pars {
 339         my( $charmatrix ) = @_;
 340     # Set up a temporary directory for all the default Phylip files.
 341     my $phylip_dir = File::Temp->newdir();
 342     # $phylip_dir->unlink_on_destroy(0);
 343     # We need an infile, and we need a command input file.
 344     open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
 345     print MATRIX $charmatrix;
 346     close MATRIX;
 347
 348     open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
 349     ## TODO any configuration parameters we want to set here
 350 #   U                 Search for best tree?  Yes
 351 #   S                        Search option?  More thorough search
 352 #   V              Number of trees to save?  100
 353 #   J     Randomize input order of species?  No. Use input order
 354 #   O                        Outgroup root?  No, use as outgroup species 1
 355 #   T              Use Threshold parsimony?  No, use ordinary parsimony
 356 #   W                       Sites weighted?  No
 357 #   M           Analyze multiple data sets?  No
 358 #   I            Input species interleaved?  Yes
 359 #   0   Terminal type (IBM PC, ANSI, none)?  ANSI
 360 #   1    Print out the data at start of run  No
 361 #   2  Print indications of progress of run  Yes
 362 #   3                        Print out tree  Yes
 363 #   4          Print out steps in each site  No
 364 #   5  Print character at all nodes of tree  No
 365 #   6       Write out trees onto tree file?  Yes
 366     print CMD "Y\n";
 367     close CMD;
 368
 369     # And then we run the program.
 370     my $program = File::Which::which( 'pars' );
 371     unless( $program && -x $program ) {
 372                 throw( "Phylip pars not found in path" );
 373     }
 374
 375     {
 376         # We need to run it in our temporary directory where we have created
 377         # all the expected files.
 378         local $CWD = $phylip_dir;
 379         my @cmd = ( $program );
 380         run \@cmd, '<', 'cmdfile', '>', '/dev/null';
 381     }
 382     # Now our output should be in 'outfile' and our tree in 'outtree',
 383     # both in the temp directory.
 384
 385     my @outtree;
 386     if( -f "$phylip_dir/outtree" ) {
 387         open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
 388         @outtree = <TREE>;
 389         close TREE;
 390     }
 391     return join( '', @outtree ) if @outtree;
 392
 393         # If we got this far, we are about to throw an error.
 394     my @error;
 395     if( -f "$phylip_dir/outfile" ) {
 396         open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
 397         @error = <OUTPUT>;
 398         close OUTPUT;
 399     } else {
 400         push( @error, "Neither outtree nor output file was produced!" );
 401     }
 402     throw( join( '', @error ) );
 403 }
 404
 405 =head2 parse_newick( $newick_string )
 406
 407 Parses the given Newick tree(s) into one or more Stemma objects with
 408 undirected graphs.
 409
 410 =cut
 411
 412 sub parse_newick {
 413     my $newick = shift;
 414     # Parse the result into a set of trees and return them.
 415     my $forest = Bio::Phylo::IO->parse(
 416         -format => 'newick',
 417         -string => $newick,
 418         );
 419     return map { _graph_from_bio( $_ ) } @{$forest->get_entities};
 420 }
 421
 422 sub _graph_from_bio {
 423     my $tree = shift;
 424     my $graph = Graph->new( 'undirected' => 1 );
 425     # Give all the intermediate anonymous nodes a name.
 426     my $i = 0;
 427     my $classes = {};
 428     foreach my $n ( @{$tree->get_terminals} ) {
 429         # The terminal nodes are our named witnesses.
 430                 $classes->{$n->get_name} = 'extant';
 431         }
 432         foreach my $n ( @{$tree->get_internals} ) {
 433         unless( defined $n->get_name && $n->get_name ne '' ) {
 434                 # Get an integer, make sure it's a unique name
 435                 while( exists $classes->{$i} ) {
 436                         $i++;
 437                 }
 438                 $n->set_name( $i++ );
 439         }
 440         $classes->{$n->get_name} = 'hypothetical';
 441     }
 442     _add_tree_children( $graph, $classes, undef, [ $tree->get_root ]);
 443     return $graph;
 444 }
 445
 446 sub _add_tree_children {
 447     my( $graph, $classes, $parent, $tree_children ) = @_;
 448     foreach my $c ( @$tree_children ) {
 449         my $child = $c->get_name;
 450         $graph->add_vertex( $child );
 451         $graph->set_vertex_attribute( $child, 'class', $classes->{$child} );
 452         $graph->add_path( $parent, $child ) if defined $parent;
 453         _add_tree_children( $graph, $classes, $child, $c->get_children() );
 454     }
 455 }
 456
 457 =head2 newick_to_svg( $newick_string )
 458
 459 Uses the FigTree utility (if installed) to transform the given Newick tree(s)
 460 into a graph visualization.
 461
 462 =cut
 463
 464 sub newick_to_svg {
 465         my $newick = shift;
 466     my $program = File::Which::which( 'figtree' );
 467     unless( -x $program ) {
 468                 throw( "FigTree commandline utility not found in path" );
 469     }
 470     my $svg;
 471     my $nfile = File::Temp->new();
 472     print $nfile $newick;
 473     close $nfile;
 474         my @cmd = ( $program, '-graphic', 'SVG', $nfile );
 475     run( \@cmd, ">", binary(), \$svg );
 476     return decode_utf8( $svg );
 477 }
 478
 479 sub throw {
 480         Text::Tradition::Error->throw(
 481                 'ident' => 'StemmaUtil error',
 482                 'message' => $_[0],
 483                 );
 484 }
 485
 486 1;
 487
 488 =head1 LICENSE
 489
 490 This package is free software and is provided "as is" without express
 491 or implied warranty.  You can redistribute it and/or modify it under
 492 the same terms as Perl itself.
 493
 494 =head1 AUTHOR
 495
 496 Tara L Andrews E<lt>aurum@cpan.orgE<gt>