lib/Bio/TreeIO/pag.pm

   1 #
   2 # BioPerl module for Bio::TreeIO::pag
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason-at-bioperl-dot-org>
   7 #
   8 # Copyright Jason Stajich
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::TreeIO::pag - Bio::TreeIO driver for Pagel format
  17
  18 =head1 SYNOPSIS
  19
  20   use Bio::TreeIO;
  21   my $in = Bio::TreeIO->new(-format => 'nexus',
  22                             -file   => 't/data/adh.mb_tree.nexus');
  23
  24   my $out = Bio::TreeIO->new(-format => 'pag');
  25   while( my $tree = $in->next_tree ) {
  26     $out->write_tree($tree);
  27   }
  28
  29 =head1 DESCRIPTION
  30
  31 Convert a Bio::TreeIO to Pagel format.
  32 More information here http://www.evolution.reading.ac.uk/index.html
  33
  34 =head1 FEEDBACK
  35
  36 =head2 Mailing Lists
  37
  38 User feedback is an integral part of the evolution of this and other
  39 Bioperl modules. Send your comments and suggestions preferably to
  40 the Bioperl mailing list.  Your participation is much appreciated.
  41
  42   bioperl-l@bioperl.org                  - General discussion
  43   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  44
  45 =head2 Support
  46
  47 Please direct usage questions or support issues to the mailing list:
  48
  49 I<bioperl-l@bioperl.org>
  50
  51 rather than to the module maintainer directly. Many experienced and
  52 reponsive experts will be able look at the problem and quickly
  53 address it. Please include a thorough description of the problem
  54 with code and data examples if at all possible.
  55
  56 =head2 Reporting Bugs
  57
  58 Report bugs to the Bioperl bug tracking system to help us keep track
  59 of the bugs and their resolution. Bug reports can be submitted via
  60 the web:
  61
  62   https://github.com/bioperl/bioperl-live/issues
  63
  64 =head1 AUTHOR - Jason Stajich
  65
  66 Email jason-at-bioperl-dot-org
  67
  68 =head1 APPENDIX
  69
  70 The rest of the documentation details each of the object methods.
  71 Internal methods are usually preceded with a _
  72
  73 =cut
  74
  75
  76 # Let the code begin...
  77
  78
  79 package Bio::TreeIO::pag;
  80 use strict;
  81
  82 our $TaxonNameLen = 10;
  83
  84 use base qw(Bio::TreeIO);
  85
  86 =head2 new
  87
  88  Title   : new
  89  Usage   : my $obj = Bio::TreeIO::pag->new();
  90  Function: Builds a new Bio::TreeIO::pag object
  91  Returns : an instance of Bio::TreeIO::pag
  92  Args    : -file/-fh for filename or filehandles
  93            -name_length for minimum name length (default = 10)
  94
  95 =cut
  96
  97 sub _initialize {
  98     my $self = shift;
  99     $self->SUPER::_initialize(@_);
 100     my ( $name_length ) = $self->_rearrange(
 101         [
 102             qw(NAME_LENGTH)
 103         ],
 104         @_
 105     );
 106     $self->name_length( defined $name_length ? $name_length : $TaxonNameLen );
 107 }
 108
 109 =head2 write_tree
 110
 111  Title   : write_tree
 112  Usage   :
 113  Function: Write a tree out in Pagel format
 114            Some options are only appropriate for bayesianmultistate and
 115            the simpler output is only proper for discrete
 116  Returns : none
 117  Args    : -no_outgroups => (number)
 118            -print_header => 0/1 (leave 0 for discrete, 1 for bayesianms)
 119            -special_node => special node - not sure what they wanted to do here
 120            -keep_outgroup => 0/1 (keep the outgroup node in the output)
 121            -outgroup_ancestor => Bio::Tree::Node (if we want to exclude or include the outgroup this is what we operate on)
 122            -tree_no       => a tree number label - only useful for BayesianMultistate
 123
 124
 125 =cut
 126
 127 sub write_tree {
 128     my ($self,$tree,@args) = @_;
 129     my ($keep_outgroup,
 130         $print_header,
 131         $no_outgroups,
 132         $special_node,
 133         $outgroup_ancestor,
 134         $tree_no) = (0,0,1);
 135     my $name_len = $self->name_length;
 136     if( @args ) {
 137         ($no_outgroups,
 138          $print_header,
 139          $special_node,
 140          $outgroup_ancestor,
 141          $tree_no,
 142          $keep_outgroup) = $self->_rearrange([qw(
 143                          NO_OUTGROUPS
 144                                                  PRINT_HEADER
 145                                                  SPECIAL_NODE
 146                                                  OUTGROUP_ANCESTOR
 147                                                  TREE_NO
 148                                                  KEEP_OUTGROUP
 149                          NAME_LENGTH)],@args);
 150     }
 151     my $newname_base = 1;
 152
 153     my $root = $tree->get_root_node;
 154     my $eps = 0.0001;
 155     my (%chars,%names);
 156     my @nodes = $tree->get_nodes;
 157     my $species_ct;
 158     my $traitct;
 159     for my $node ( @nodes ) {
 160         if ((defined $special_node) && ($node eq $special_node)) {
 161             my $no_of_tree_nodes = scalar(@nodes);
 162             my $node_name = sprintf("N%d",$no_of_tree_nodes+1);
 163             $names{$node->internal_id} = $node_name;
 164
 165         } elsif ($node->is_Leaf) {
 166             $species_ct++;
 167
 168             my $node_name = $node->id;
 169             if( length($node_name)> $name_len ) {
 170                 $self->warn( "Found a taxon name longer than $name_len letters, \n",
 171                              "name will be abbreviated.\n");
 172                 $node_name = substr($node_name, 0,$name_len);
 173             } else {
 174                 # $node_name = sprintf("%-".$TaxonNameLen."s",$node_name);
 175             }
 176             $names{$node->internal_id} = $node_name;
 177             my @tags = sort $node->get_all_tags;
 178             my @charstates = map { ($node->get_tag_values($_))[0] } @tags;
 179             $traitct = scalar @charstates unless defined $traitct;
 180             $chars{$node->internal_id} = [@charstates];
 181         } else {
 182             $names{$node->internal_id} = sprintf("N%d", $newname_base++);
 183         }
 184     }
 185
 186     # generate PAG representation
 187     if( $print_header ) {
 188         if ($keep_outgroup) {
 189             $self->_print(sprintf("%d %d\n",$species_ct,$traitct));
 190         } else {
 191             $self->_print( sprintf("%d %d\n",$species_ct-$no_outgroups,$traitct));
 192         }
 193     }
 194
 195     my @ancestors = ();
 196     if ($keep_outgroup) {
 197         push @ancestors, $root;
 198     } else {
 199         push @ancestors, ( $root, $outgroup_ancestor);
 200     }
 201     my @rest;
 202     foreach my $node (@nodes) {
 203         my $i = 0;
 204         foreach my $anc (@ancestors) {
 205             if ($anc && $node eq $anc) { $i = 1; last }
 206         }
 207         unless ($i > 0) {       # root not given in PAG
 208             my $current_name = $names{$node->internal_id};
 209             my $branch_length_to_output;
 210             if ($node->branch_length < $eps) {
 211                 my $msg_nodename = $current_name;
 212                 $msg_nodename =~ s/\s+$//;
 213                 warn( "TREE $tree_no, node \"$msg_nodename\": branch too ",
 214                       "short (", $node->branch_length, "): increasing length to ",
 215                       "$eps\n");
 216                 $branch_length_to_output = $eps;
 217             } else {
 218                 $branch_length_to_output = $node->branch_length;
 219             }
 220             my @line = ( $current_name,
 221                          $names{$node->ancestor->internal_id},
 222                          $branch_length_to_output);
 223
 224             if ($node->is_Leaf) {
 225                 push @line, @{$chars{$node->internal_id}};
 226                 $self->_print(join(',', @line),"\n");
 227             } else {
 228                 push @rest, \@line;
 229             }
 230         }
 231     }
 232     for ( @rest ) {
 233         $self->_print(join(',', @$_),"\n");
 234     }
 235 }
 236
 237 =head2 next_tree
 238
 239  Title   : next_tree
 240  Usage   :
 241  Function:
 242  Example :
 243  Returns :
 244  Args    :
 245
 246
 247 =cut
 248
 249 sub next_tree{
 250    my ($self,@args) = @_;
 251    $self->throw_not_implemented();
 252 }
 253
 254 =head2 name_length
 255
 256  Title   : name_length
 257  Usage   : $self->name_length(20);
 258  Function: set minimum taxon name length
 259  Returns : integer (length of name)
 260  Args    : integer
 261
 262 =cut
 263
 264 sub name_length {
 265     my ($self, $val) = @_;
 266     return $self->{'name_len'} = $val if $val;
 267     return $self->{'name_len'};
 268 }
 269
 270 1;