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 IPC::Run qw/ run binary /;
  14 use Text::Tradition::Error;
  15 use Text::Tradition::Stemma;
  16 @EXPORT_OK = qw/ character_input phylip_pars parse_newick newick_to_svg /;
  17
  18 =head1 NAME
  19
  20 Text::Tradition::StemmaUtil - standalone utilities for distance tree calculations
  21
  22 =head1 DESCRIPTION
  23
  24 This package contains a set of utilities for running phylogenetic analysis on
  25 text collations.
  26
  27 =head1 SUBROUTINES
  28
  29 =head2 character_input( $tradition, $opts )
  30
  31 Returns a character matrix string suitable for Phylip programs, which
  32 corresponds to the given alignment table.  See Text::Tradition::Collation
  33 for a description of the alignment table format. Options include:
  34
  35 =over
  36
  37 =item * exclude_layer - Exclude layered witnesses from the character input,
  38 using only the 'main' text of the witnesses in the tradition.
  39
  40 =item * collapse - A reference to an array of relationship names that should
  41 be treated as equivalent for the purposes of generating the character matrix.
  42
  43 =back
  44
  45 =cut
  46
  47 sub character_input {
  48     my ( $tradition, $opts ) = @_;
  49     my $table = $tradition->collation->alignment_table;
  50     if( $opts->{exclude_layer} ) {
  51         # Filter out all alignment table rows that do not correspond
  52         # to a named witness - these are the layered witnesses.
  53         my $newtable = { alignment => [] };
  54         foreach my $row ( $table->{alignment} ) {
  55                 if( $tradition->has_witness( $row->{witness} ) ) {
  56                         push( @{$newtable->{alignment}}, $row );
  57                 }
  58         }
  59         $table = $newtable;
  60     }
  61     my $character_matrix = _make_character_matrix( $table, $opts );
  62     my $input = '';
  63     my $rows = scalar @{$character_matrix};
  64     my $columns = scalar @{$character_matrix->[0]} - 1;
  65     $input .= "\t$rows\t$columns\n";
  66     foreach my $row ( @{$character_matrix} ) {
  67         $input .= join( '', @$row ) . "\n";
  68     }
  69     return $input;
  70 }
  71
  72 sub _make_character_matrix {
  73     my( $table, $opts ) = @_;
  74     # Push the names of the witnesses to initialize the rows of the matrix.
  75     my @matrix = map { [ _normalize_witname( $_->{'witness'} ) ] }
  76                                 @{$table->{'alignment'}};
  77     foreach my $token_index ( 0 .. $table->{'length'} - 1) {
  78         my @pos_tokens = map { $_->{'tokens'}->[$token_index] }
  79                                                         @{$table->{'alignment'}};
  80         my @pos_readings = map { $_ ? $_->{'t'} : $_ } @pos_tokens;
  81         my @chars = _convert_characters( \@pos_readings, $opts );
  82         foreach my $idx ( 0 .. $#matrix ) {
  83             push( @{$matrix[$idx]}, $chars[$idx] );
  84         }
  85     }
  86     return \@matrix;
  87 }
  88
  89 # Helper function to make the witness name something legal for pars
  90
  91 sub _normalize_witname {
  92     my( $witname ) = @_;
  93     $witname =~ s/\s+/ /g;
  94     $witname =~ s/[\[\]\(\)\:;,]//g;
  95     $witname = substr( $witname, 0, 10 );
  96     return sprintf( "%-10s", $witname );
  97 }
  98
  99 sub _convert_characters {
 100     my( $row, $opts ) = @_;
 101     # This is a simple algorithm that treats every reading as different.
 102     # Eventually we will want to be able to specify how relationships
 103     # affect the character matrix.
 104     my %unique = ( '__UNDEF__' => 'X',
 105                    '#LACUNA#'  => '?',
 106                  );
 107     my %equivalent;
 108     my %count;
 109     my $ctr = 0;
 110     foreach my $rdg ( @$row ) {
 111         next unless $rdg;
 112         next if $rdg->is_lacuna;
 113                 next if exists $unique{$rdg->text};
 114                 if( ref( $opts->{'collapse'} ) eq 'ARRAY' ) {
 115                         my @exclude_types = @{$opts->{'collapse'}};
 116                         my @set = $rdg->related_readings( sub { my $rel = shift;
 117                                 $rel->colocated && grep { $rel->type eq $_ } @exclude_types } );
 118                         push( @set, $rdg );
 119                         my $char = chr( 65 + $ctr++ );
 120                         map { $unique{$_->text} = $char } @set;
 121                         $count{$rdg->text} += scalar @set;
 122                 } else {
 123                         $unique{$rdg->text} = chr( 65 + $ctr++ );
 124                         $count{$rdg->text}++;
 125                 }
 126     }
 127     # Try to keep variants under 8 by lacunizing any singletons.
 128     if( scalar( keys %unique ) > 8 ) {
 129                 foreach my $word ( keys %count ) {
 130                         if( $count{$word} == 1 ) {
 131                                 $unique{$word} = '?';
 132                         }
 133                 }
 134     }
 135     my %u = reverse %unique;
 136     if( scalar( keys %u ) > 8 ) {
 137         warn "Have more than 8 variants on this location; phylip will break";
 138     }
 139     my @chars = map { $_ ? $unique{$_->text} : $unique{'__UNDEF__' } } @$row;
 140     return @chars;
 141 }
 142
 143 =head2 phylip_pars( $character_matrix )
 144
 145 Runs Phylip Pars on the given character matrix.  Returns results in Newick format.
 146
 147 =cut
 148
 149 sub phylip_pars {
 150         my( $charmatrix ) = @_;
 151     # Set up a temporary directory for all the default Phylip files.
 152     my $phylip_dir = File::Temp->newdir();
 153     # $phylip_dir->unlink_on_destroy(0);
 154     # We need an infile, and we need a command input file.
 155     open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
 156     print MATRIX $charmatrix;
 157     close MATRIX;
 158
 159     open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
 160     ## TODO any configuration parameters we want to set here
 161 #   U                 Search for best tree?  Yes
 162 #   S                        Search option?  More thorough search
 163 #   V              Number of trees to save?  100
 164 #   J     Randomize input order of species?  No. Use input order
 165 #   O                        Outgroup root?  No, use as outgroup species 1
 166 #   T              Use Threshold parsimony?  No, use ordinary parsimony
 167 #   W                       Sites weighted?  No
 168 #   M           Analyze multiple data sets?  No
 169 #   I            Input species interleaved?  Yes
 170 #   0   Terminal type (IBM PC, ANSI, none)?  ANSI
 171 #   1    Print out the data at start of run  No
 172 #   2  Print indications of progress of run  Yes
 173 #   3                        Print out tree  Yes
 174 #   4          Print out steps in each site  No
 175 #   5  Print character at all nodes of tree  No
 176 #   6       Write out trees onto tree file?  Yes
 177     print CMD "Y\n";
 178     close CMD;
 179
 180     # And then we run the program.
 181     my $program = File::Which::which( 'pars' );
 182     unless( -x $program ) {
 183                 throw( "Phylip pars not found in path" );
 184     }
 185
 186     {
 187         # We need to run it in our temporary directory where we have created
 188         # all the expected files.
 189         local $CWD = $phylip_dir;
 190         my @cmd = ( $program );
 191         run \@cmd, '<', 'cmdfile', '>', '/dev/null';
 192     }
 193     # Now our output should be in 'outfile' and our tree in 'outtree',
 194     # both in the temp directory.
 195
 196     my @outtree;
 197     if( -f "$phylip_dir/outtree" ) {
 198         open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
 199         @outtree = <TREE>;
 200         close TREE;
 201     }
 202     return join( '', @outtree ) if @outtree;
 203
 204         # If we got this far, we are about to throw an error.
 205     my @error;
 206     if( -f "$phylip_dir/outfile" ) {
 207         open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
 208         @error = <OUTPUT>;
 209         close OUTPUT;
 210     } else {
 211         push( @error, "Neither outtree nor output file was produced!" );
 212     }
 213     throw( join( '', @error ) );
 214 }
 215
 216 =head2 parse_newick( $newick_string )
 217
 218 Parses the given Newick tree(s) into one or more Stemma objects with
 219 undirected graphs.
 220
 221 =cut
 222
 223 sub parse_newick {
 224     my $newick = shift;
 225     my @stemmata;
 226     # Parse the result into a tree
 227     my $forest = Bio::Phylo::IO->parse(
 228         -format => 'newick',
 229         -string => $newick,
 230         );
 231     # Turn the tree into a graph, starting with the root node
 232     foreach my $tree ( @{$forest->get_entities} ) {
 233         my $stemma = Text::Tradition::Stemma->new(
 234                 graph => _graph_from_bio( $tree ),
 235                 is_undirected => 1 );
 236         push( @stemmata, $stemma );
 237     }
 238     return \@stemmata;
 239 }
 240
 241 sub _graph_from_bio {
 242     my $tree = shift;
 243     my $graph = Graph->new( 'undirected' => 1 );
 244     # Give all the intermediate anonymous nodes a name.
 245     my $i = 0;
 246     my $classes = {};
 247     foreach my $n ( @{$tree->get_terminals} ) {
 248         # The terminal nodes are our named witnesses.
 249                 $classes->{$n->get_name} = 'extant';
 250         }
 251         foreach my $n ( @{$tree->get_internals} ) {
 252         unless( defined $n->get_name && $n->get_name ne '' ) {
 253                 # Get an integer, make sure it's a unique name
 254                 while( exists $classes->{$i} ) {
 255                         $i++;
 256                 }
 257                 $n->set_name( $i++ );
 258         }
 259         $classes->{$n->get_name} = 'hypothetical';
 260     }
 261     _add_tree_children( $graph, $classes, undef, [ $tree->get_root ]);
 262     return $graph;
 263 }
 264
 265 sub _add_tree_children {
 266     my( $graph, $classes, $parent, $tree_children ) = @_;
 267     foreach my $c ( @$tree_children ) {
 268         my $child = $c->get_name;
 269         $graph->add_vertex( $child );
 270         $graph->set_vertex_attribute( $child, 'class', $classes->{$child} );
 271         $graph->add_path( $parent, $child ) if defined $parent;
 272         _add_tree_children( $graph, $classes, $child, $c->get_children() );
 273     }
 274 }
 275
 276 =head2 newick_to_svg( $newick_string )
 277
 278 Uses the FigTree utility (if installed) to transform the given Newick tree(s)
 279 into a graph visualization.
 280
 281 =cut
 282
 283 sub newick_to_svg {
 284         my $newick = shift;
 285     my $program = File::Which::which( 'figtree' );
 286     unless( -x $program ) {
 287                 throw( "FigTree commandline utility not found in path" );
 288     }
 289     my $svg;
 290     my $nfile = File::Temp->new();
 291     print $nfile $newick;
 292     close $nfile;
 293         my @cmd = ( $program, '-graphic', 'SVG', $nfile );
 294     run( \@cmd, ">", binary(), \$svg );
 295     return decode_utf8( $svg );
 296 }
 297
 298 sub throw {
 299         Text::Tradition::Error->throw(
 300                 'ident' => 'StemmaUtil error',
 301                 'message' => $_[0],
 302                 );
 303 }
 304
 305 1;
 306
 307 =head1 LICENSE
 308
 309 This package is free software and is provided "as is" without express
 310 or implied warranty.  You can redistribute it and/or modify it under
 311 the same terms as Perl itself.
 312
 313 =head1 AUTHOR
 314
 315 Tara L Andrews E<lt>aurum@cpan.orgE<gt>