[scpubgit/stemmatology.git] / lib / Text / Tradition / StemmaUtil.pm

package Text::Tradition::StemmaUtil;

use strict;
use warnings;
use Exporter 'import';
use vars qw/ @EXPORT_OK /;
use Bio::Phylo::IO;
use Encode qw( decode_utf8 );
use File::chdir;
use File::Temp;
use File::Which;
use Graph;
use Graph::Reader::Dot;
use IPC::Run qw/ run binary /;
use Text::Tradition::Error;
@EXPORT_OK = qw/ make_character_matrix character_input phylip_pars 
				 parse_newick newick_to_svg /;

=head1 NAME

Text::Tradition::StemmaUtil - standalone utilities for distance tree calculations

=head1 DESCRIPTION

This package contains a set of utilities for running phylogenetic analysis on
text collations.

=head1 SUBROUTINES

=cut

sub _make_character_matrix {
    my( $table ) = @_;
    # Push the names of the witnesses to initialize the rows of the matrix.
    my @matrix = map { [ _normalize_witname( $_->{'witness'} ) ] } 
    				@{$table->{'alignment'}};
    foreach my $token_index ( 0 .. $table->{'length'} - 1) {
        # First implementation: make dumb alignment table, caring about
        # nothing except which reading is in which position.
        my @pos_readings = map { $_->{'tokens'}->[$token_index] }
        						@{$table->{'alignment'}};
        my @pos_text = map { $_ ? $_->{'t'} : $_ } @pos_readings;
        my @chars = _convert_characters( \@pos_text );
        foreach my $idx ( 0 .. $#matrix ) {
            push( @{$matrix[$idx]}, $chars[$idx] );
        }
    }
    return \@matrix;
} 

# Helper function to make the witness name something legal for pars

sub _normalize_witname {
    my( $witname ) = @_;
    $witname =~ s/\s+/ /g;
    $witname =~ s/[\[\]\(\)\:;,]//g;
    $witname = substr( $witname, 0, 10 );
    return sprintf( "%-10s", $witname );
}

sub _convert_characters {
    my $row = shift;
    # This is a simple algorithm that treats every reading as different.
    # Eventually we will want to be able to specify how relationships
    # affect the character matrix.
    my %unique = ( '__UNDEF__' => 'X',
                   '#LACUNA#'  => '?',
                 );
    my %count;
    my $ctr = 0;
    foreach my $word ( @$row ) {
        if( $word && !exists $unique{$word} ) {
            $unique{$word} = chr( 65 + $ctr );
            $ctr++;
        }
        $count{$word}++ if $word;
    }
    # Try to keep variants under 8 by lacunizing any singletons.
    if( scalar( keys %unique ) > 8 ) {
		foreach my $word ( keys %count ) {
			if( $count{$word} == 1 ) {
				$unique{$word} = '?';
			}
		}
    }
    my %u = reverse %unique;
    if( scalar( keys %u ) > 8 ) {
        warn "Have more than 8 variants on this location; phylip will break";
    }
    my @chars = map { $_ ? $unique{$_} : $unique{'__UNDEF__' } } @$row;
    return @chars;
}

=head2 character_input( $alignment_table )

Returns a character matrix string suitable for Phylip programs, which 
corresponds to the given alignment table.  See Text::Tradition::Collation 
for a description of the alignment table format.

=cut

sub character_input {
    my $table = shift;
    my $character_matrix = _make_character_matrix( $table );
    my $input = '';
    my $rows = scalar @{$character_matrix};
    my $columns = scalar @{$character_matrix->[0]} - 1;
    $input .= "\t$rows\t$columns\n";
    foreach my $row ( @{$character_matrix} ) {
        $input .= join( '', @$row ) . "\n";
    }
    return $input;
}

=head2 phylip_pars( $character_matrix )

Runs Phylip Pars on the given character matrix.  Returns results in Newick format.

=cut

sub phylip_pars {
	my( $charmatrix ) = @_;
    # Set up a temporary directory for all the default Phylip files.
    my $phylip_dir = File::Temp->newdir();
    # $phylip_dir->unlink_on_destroy(0);
    # We need an infile, and we need a command input file.
    open( MATRIX, ">$phylip_dir/infile" ) or die "Could not write $phylip_dir/infile";
    print MATRIX $charmatrix;
    close MATRIX;

    open( CMD, ">$phylip_dir/cmdfile" ) or die "Could not write $phylip_dir/cmdfile";
    ## TODO any configuration parameters we want to set here
#   U                 Search for best tree?  Yes
#   S                        Search option?  More thorough search
#   V              Number of trees to save?  100
#   J     Randomize input order of species?  No. Use input order
#   O                        Outgroup root?  No, use as outgroup species 1
#   T              Use Threshold parsimony?  No, use ordinary parsimony
#   W                       Sites weighted?  No
#   M           Analyze multiple data sets?  No
#   I            Input species interleaved?  Yes
#   0   Terminal type (IBM PC, ANSI, none)?  ANSI
#   1    Print out the data at start of run  No
#   2  Print indications of progress of run  Yes
#   3                        Print out tree  Yes
#   4          Print out steps in each site  No
#   5  Print character at all nodes of tree  No
#   6       Write out trees onto tree file?  Yes
    print CMD "Y\n";
    close CMD;

    # And then we run the program.
    my $program = File::Which::which( 'pars' );
    unless( -x $program ) {
		throw( "Phylip pars not found in path" );
    }

    {
        # We need to run it in our temporary directory where we have created
        # all the expected files.
        local $CWD = $phylip_dir;
        my @cmd = ( $program );
        run \@cmd, '<', 'cmdfile', '>', '/dev/null';
    }
    # Now our output should be in 'outfile' and our tree in 'outtree',
    # both in the temp directory.

    my @outtree;
    if( -f "$phylip_dir/outtree" ) {
        open( TREE, "$phylip_dir/outtree" ) or die "Could not open outtree for read";
        @outtree = <TREE>;
        close TREE;
    }
    return join( '', @outtree ) if @outtree;

	# If we got this far, we are about to throw an error.
    my @error;
    if( -f "$phylip_dir/outfile" ) {
        open( OUTPUT, "$phylip_dir/outfile" ) or die "Could not open output for read";
        @error = <OUTPUT>;
        close OUTPUT;
    } else {
        push( @error, "Neither outtree nor output file was produced!" );
    }
    throw( join( '', @error ) );
}

=head2 parse_newick( $newick_string )

Parses the given Newick tree(s) into one or more undirected Graph objects.

=cut

sub parse_newick {
    my $newick = shift;
    my @trees;
    # Parse the result into a tree
    my $forest = Bio::Phylo::IO->parse( 
        -format => 'newick',
        -string => $newick,
        );
    # Turn the tree into a graph, starting with the root node
    foreach my $tree ( @{$forest->get_entities} ) {
        push( @trees, _graph_from_bio( $tree ) );
    }
    return \@trees;
}

=head2 newick_to_svg( $newick_string )

Uses the FigTree utility (if installed) to transform the given Newick tree(s)
into a graph visualization.

=cut

sub newick_to_svg {
	my $newick = shift;
    my $program = File::Which::which( 'figtree' );
    unless( -x $program ) {
		throw( "FigTree commandline utility not found in path" );
    }
    my $svg;
    my $nfile = File::Temp->new();
    print $nfile $newick;
    close $nfile;
	my @cmd = ( $program, '-graphic', 'SVG', $nfile );
    run( \@cmd, ">", binary(), \$svg );
    return decode_utf8( $svg );
}

sub _graph_from_bio {
    my $tree = shift;
    my $graph = Graph->new( 'undirected' => 1 );
    # Give all the intermediate anonymous nodes a name.
    my $i = 0;
    foreach my $n ( @{$tree->get_entities} ) {
        next if $n->get_name;
        $n->set_name( $i++ );
    }
    my $root = $tree->get_root->get_name;
    $graph->add_vertex( $root );
    _add_tree_children( $graph, $root, $tree->get_root->get_children() );
    return $graph;
}

sub _add_tree_children {
    my( $graph, $parent, $tree_children ) = @_;
    foreach my $c ( @$tree_children ) {
        my $child = $c->get_name;
        $graph->add_vertex( $child );
        $graph->add_path( $parent, $child );
        _add_tree_children( $graph, $child, $c->get_children() );
    }
}

sub throw {
	Text::Tradition::Error->throw( 
		'ident' => 'StemmaUtil error',
		'message' => $_[0],
		);
}

1;

=head1 LICENSE

This package is free software and is provided "as is" without express
or implied warranty.  You can redistribute it and/or modify it under
the same terms as Perl itself.

=head1 AUTHOR

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