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
  81 use strict;
  82
  83 our $TaxonNameLen = 10;
  84
  85 use base qw(Bio::TreeIO);
  86
  87 =head2 new
  88
  89  Title   : new
  90  Usage   : my $obj = Bio::TreeIO::pag->new();
  91  Function: Builds a new Bio::TreeIO::pag object
  92  Returns : an instance of Bio::TreeIO::pag
  93  Args    : -file/-fh for filename or filehandles
  94            -name_length for minimum name length (default = 10)
  95
  96 =cut
  97
  98 sub _initialize {
  99     my $self = shift;
 100     $self->SUPER::_initialize(@_);
 101     my ( $name_length ) = $self->_rearrange(
 102         [
 103             qw(NAME_LENGTH)
 104         ],
 105         @_
 106     );
 107     $self->name_length( defined $name_length ? $name_length : $TaxonNameLen );
 108 }
 109
 110 =head2 write_tree
 111
 112  Title   : write_tree
 113  Usage   :
 114  Function: Write a tree out in Pagel format
 115            Some options are only appropriate for bayesianmultistate and
 116            the simpler output is only proper for discrete
 117  Returns : none
 118  Args    : -no_outgroups => (number)
 119            -print_header => 0/1 (leave 0 for discrete, 1 for bayesianms)
 120            -special_node => special node - not sure what they wanted to do here
 121            -keep_outgroup => 0/1 (keep the outgroup node in the output)
 122            -outgroup_ancestor => Bio::Tree::Node (if we want to exclude or include the outgroup this is what we operate on)
 123            -tree_no       => a tree number label - only useful for BayesianMultistate
 124
 125
 126 =cut
 127
 128 sub write_tree {
 129     my ($self,$tree,@args) = @_;
 130     my ($keep_outgroup,
 131         $print_header,
 132         $no_outgroups,
 133         $special_node,
 134         $outgroup_ancestor,
 135         $tree_no) = (0,0,1);
 136     my $name_len = $self->name_length;
 137     if( @args ) {
 138         ($no_outgroups,
 139          $print_header,
 140          $special_node,
 141          $outgroup_ancestor,
 142          $tree_no,
 143          $keep_outgroup) = $self->_rearrange([qw(
 144                          NO_OUTGROUPS
 145                                                  PRINT_HEADER
 146                                                  SPECIAL_NODE
 147                                                  OUTGROUP_ANCESTOR
 148                                                  TREE_NO
 149                                                  KEEP_OUTGROUP
 150                          NAME_LENGTH)],@args);
 151     }
 152     my $newname_base = 1;
 153
 154     my $root = $tree->get_root_node;
 155     my $eps = 0.0001;
 156     my (%chars,%names);
 157     my @nodes = $tree->get_nodes;
 158     my $species_ct;
 159     my $traitct;
 160     for my $node ( @nodes ) {
 161         if ((defined $special_node) && ($node eq $special_node)) {
 162             my $no_of_tree_nodes = scalar(@nodes);
 163             my $node_name = sprintf("N%d",$no_of_tree_nodes+1);
 164             $names{$node->internal_id} = $node_name;
 165
 166         } elsif ($node->is_Leaf) {
 167             $species_ct++;
 168
 169             my $node_name = $node->id;
 170             if( length($node_name)> $name_len ) {
 171                 $self->warn( "Found a taxon name longer than $name_len letters, \n",
 172                              "name will be abbreviated.\n");
 173                 $node_name = substr($node_name, 0,$name_len);
 174             } else {
 175                 # $node_name = sprintf("%-".$TaxonNameLen."s",$node_name);
 176             }
 177             $names{$node->internal_id} = $node_name;
 178             my @tags = sort $node->get_all_tags;
 179             my @charstates = map { ($node->get_tag_values($_))[0] } @tags;
 180             $traitct = scalar @charstates unless defined $traitct;
 181             $chars{$node->internal_id} = [@charstates];
 182         } else {
 183             $names{$node->internal_id} = sprintf("N%d", $newname_base++);
 184         }
 185     }
 186
 187     # generate PAG representation
 188     if( $print_header ) {
 189         if ($keep_outgroup) {
 190             $self->_print(sprintf("%d %d\n",$species_ct,$traitct));
 191         } else {
 192             $self->_print( sprintf("%d %d\n",$species_ct-$no_outgroups,$traitct));
 193         }
 194     }
 195
 196     my @ancestors = ();
 197     if ($keep_outgroup) {
 198         push @ancestors, $root;
 199     } else {
 200         push @ancestors, ( $root, $outgroup_ancestor);
 201     }
 202     my @rest;
 203     foreach my $node (@nodes) {
 204         my $i = 0;
 205         foreach my $anc (@ancestors) {
 206             if ($anc && $node eq $anc) { $i = 1; last }
 207         }
 208         unless ($i > 0) {       # root not given in PAG
 209             my $current_name = $names{$node->internal_id};
 210             my $branch_length_to_output;
 211             if ($node->branch_length < $eps) {
 212                 my $msg_nodename = $current_name;
 213                 $msg_nodename =~ s/\s+$//;
 214                 warn( "TREE $tree_no, node \"$msg_nodename\": branch too ",
 215                       "short (", $node->branch_length, "): increasing length to ",
 216                       "$eps\n");
 217                 $branch_length_to_output = $eps;
 218             } else {
 219                 $branch_length_to_output = $node->branch_length;
 220             }
 221             my @line = ( $current_name,
 222                          $names{$node->ancestor->internal_id},
 223                          $branch_length_to_output);
 224
 225             if ($node->is_Leaf) {
 226                 push @line, @{$chars{$node->internal_id}};
 227                 $self->_print(join(',', @line),"\n");
 228             } else {
 229                 push @rest, \@line;
 230             }
 231         }
 232     }
 233     for ( @rest ) {
 234         $self->_print(join(',', @$_),"\n");
 235     }
 236 }
 237
 238 =head2 next_tree
 239
 240  Title   : next_tree
 241  Usage   :
 242  Function:
 243  Example :
 244  Returns :
 245  Args    :
 246
 247
 248 =cut
 249
 250 sub next_tree{
 251    my ($self,@args) = @_;
 252    $self->throw_not_implemented();
 253 }
 254
 255 =head2 name_length
 256
 257  Title   : name_length
 258  Usage   : $self->name_length(20);
 259  Function: set minimum taxon name length
 260  Returns : integer (length of name)
 261  Args    : integer
 262
 263 =cut
 264
 265 sub name_length {
 266     my ($self, $val) = @_;
 267     return $self->{'name_len'} = $val if $val;
 268     return $self->{'name_len'};
 269 }
 270
 271 1;