4 de julio de 2017

rooting and laddering Newick trees

Buenas,
esta entrada es para compartir un script Perl para enraizar árboles, muchos árboles, de manera automática, y ordenar los nodos por distancia. Pongo el título en inglés para los buscadores.
Rubén Sancho y yo estamos probando el software GRAMPA, que requiere una colección de árboles de genes precalculados en formato Newick, pero además los quiere enraizados y ordenados de manera ascendente. Como son más de mil árboles no lo queríamos hacer uno a uno con FigTree, por poner un ejemplo; queríamos automatizarlo y para ello aprovechamos los módulos Bio::TreeIO y Bio::Phylo. El primero es de Bioperl y ya lo tenía instalado, el segundo lo instalé con: $ sudo cpanm Bio::Phylo

El árbol de ejemplo, contenido en el fichero unrooted.ph, es:

(1_Sbic:0.047,2_Osat:0.068,(((((((((((3_B422:0.007,(((4_Barb:0.000,((5_Bret:0.003,(6_BsyE:0.000,7_BsyE:0.000):0.000):0.000,8_BsyC:0.000):0.000):0.000,9_Bpho:0.002):0.000,(10_BsyC:0.000,11_BsyC:0.000):0.001):0.000):0.008,((12_B422:0.005,13_Bpho:0.005):0.002,14_Barb:0.000):0.005):0.000,((((15_Bdis:0.000,16_Bdis:0.000):0.000,17_Bmex:0.033):0.013,18_Bhyb:0.001):0.014,(19_Bboi:0.021,((20_Bmex:0.017,21_Bsta:0.034):0.012,22_Bsta:0.000):0.003):0.020):0.012):0.002,23_Barb:0.000):0.000,24_Brup:0.012):0.000,25_BsyG:0.000):0.000,(26_Bpin:0.000,27_BsyG:0.000):0.000):0.011,28_Bsta:0.000):0.053,29_BsyG:0.000):0.036,(30_Barb:0.005,(31_Bpho:0.024,(32_BsyC:0.000,33_BsyE:0.000):0.027):0.029):0.000):0.005,34_Bboi:0.030):0.035);

Invocando el script lo enraizamos con el taxón elegido, el outgroup '_Sbic':

$ perl reroot_tree.pl unrooted.ph > rooted.ph

Obtenemos el fichero rooted.ph, que podemos visualizar con FigTree:

Árbol enraizado y con orden creciente de nodos.
Éste es el código de reroot_tree.pl :

 #!/usr/bin/perl -w  
   
 # Re-roots an input Newick tree with a user-defined outgroup and   
 # prints the resulting tree in ascending or descending node order.  
 # Based on https://github.com/phac-nml/snvphyl-tools/blob/master/rearrange_snv_matrix.pl  
   
 use strict;  
 use Bio::TreeIO;  
 use Bio::Phylo::IO;  
 use Bio::Phylo::Forest::Tree;  
    
 my $OUTGROUPSTRING = '_Sbic'; # change as needed  
 my $NODEORDER = 1; # 1:increasing, 0:decreasing  
 my ($outfound,$outnode,$sorted_newick) = (0);  
   
 die "# usage: $0 <tree.newick>\n" if(!$ARGV[0]);  
   
 # read input tree  
 my $input = new Bio::TreeIO(-file=>$ARGV[0],-format=>'newick');  
 my $intree = $input->next_tree();  
   
 # find outgroup taxon  
 for my $node ( $intree->get_nodes() ) {   
  if(defined($node->id()) && $node->id() =~ m/$OUTGROUPSTRING/) {  
   $outnode = $node;   
   $outfound = 1;  
   last;  
  }  
 }  
 if($outfound == 0) {  
  die "# cannot find outgroup $OUTGROUPSTRING in input tree $ARGV[0]\n";  
 }  
   
 # root in outgroup and sort in increasing order  
 $intree->reroot_at_midpoint($outnode);  
   
 # sort nodes in defined order and print  
 my $unsorted_tree = Bio::Phylo::IO->parse(  
  '-string' => $intree->as_text('newick'),  
  '-format' => 'newick'  
 )->first();  
   
 $unsorted_tree->ladderize($NODEORDER);  
   
 $sorted_newick = $unsorted_tree->to_newick();  
 $sorted_newick =~ s/'//g;  
   
 print $sorted_newick;  


Un saludo,
Bruno

No hay comentarios:

Publicar un comentario