Introducción al uso de módulos en Perl con aplicaciones en bioinformática

Presentación

Este material didáctico corresponde al último tema del Curso Fundamental de Posgrado ofrecido en el marco de los Programas de Posgrado en Ciencias Bioquímicas y Ciencias Biomédicas de la Universidad Nacional Autómoma de México (UNAM), semestre 2020-1, e impartido en el Centro de Ciencias Genómicas (CCG), Cuernavaca, Morelos, México.

Módulos del CPAN, requeridos para este tema

• \( Statistics::Descriptive \)
• \( Statistics::Distributions \)
• \( Bio::Perl \)
• \( cpanm \)

En Ubuntu, los puedes instalar usando:

sudo apt install perldoc cpanminus bioperl libstatistics-descriptive-perl libstatistics-distributions-perl

Estos módulos están en el directorio \( lib \) del repositorio GitHub.

Algunos de los módulos que vamos a desarrollar en este tema, serán usados para correr tuberías de análisis filogenéticos. Para ello necesitamos contar con los siguientes binarios instalados en nuestro sistema:

1. Software para alineamientos múltiples

sudo apt install muscle clustalo 

1. Software para estima filogenética

sudo apt install fasttree phyml

Estos binarios (Linux, 64bit) están en el directorio \( bin \) del repositorio GitHub.

¿Qué es un paquete?

Perl usa paquetes para particionar o dividir en compartimentos el espacio global de identificadores, es decir, de los nombres de variables, subrutinas, descriptores de archivos y formatos.

Cuando ejecutamos un script de Perl, todas las variables y nombres de subrutinas o descriptores de archivos (file handles) están definidos en el paquete principal, llamado package main.

Cada identificador global tiene de hecho dos partes: su nombre de paquete y su identificador. Estas dos “piezas” están separadas la una de la otra por un doble doble-punto ::

Así por ejemplo, la variable $stats::x es la variable global $x, que reside en el paquete \( stats \).

$main::x es la variable $x definida en el paquete \( main \). Ambas variables son completamente independientes.

¿Qué es un módulo?

De manera muy simple, podríamos definir a un módulo como la unidad de re-uso y distribución de código en Perl. Más formalmente, se trata de un paquete definido en un archivo que lleva el mismo nombre del paquete y tiene la terminación \( .pm \), es decir, algo como NombrePaquete.pm

Los módulos contienen funciones relacionadas y están diseñados para ser usados por programas que llaman a dichos módulos, o incluso otros módulos que llaman a las funciones contenidas en diversos módulos externos.

Cada módulo tiene una interfaz pública, que permite el acceso a las funciones y variables definidas en ellos. Llamamos clases a los módulos con interfaz dirigida a objetos. Como tales, son paquetes

Programa paquetes1.pl - uso de paquetes para particionar el espacio de símbolos

Este miniprograma ilustra el concepto de paquetes como divisores del espacio de nombres para evitar colisiones entre símbolos.

#!/usr/bin/env perl

# programa: paquetes1.pl
# El concepto de paquetes como divisores del espacio de nombres para evitar colisiones entre símbolos

# estamos en el paquete principal "main"
$secuencia = 'MVLLIATG';

# declaramos un nuevo paquete y declaramos una variable en dicho paquete
package Seq1;
$secuencia = 'AAAAAAAA';

# declaramos un segundo paquete y declaramos una variable en dicho paquete
package Seq2;
$secuencia = 'CCCCCCCC';

# regresamos al paquete main e imprimimos los valores de las 3 variables escalares secuencia
package main;

# imprimimos las secuencias, usando las variables completamente calificadas
print "\n\tmain es: $secuencia"
     ."\n\tSeq1 es: $Seq1::secuencia"
     ."\n\tSeq2 es: $Seq2::secuencia\n\n";

Programa paquetes1.pl - uso de paquetes para particionar el espacio de símbolos

• De nuevo la llamada a \( print() \)

print "\n\tmain es: $secuencia"
     ."\n\tSeq1 es: $Seq1::secuencia"
     ."\n\tSeq2 es: $Seq2::secuencia\n\n";

• Veamos la salida del programa

    main es: MVLLIATG
    Seq1 es: AAAAAAAA
    Seq2 es: CCCCCCCC

Programa paquetes2.pl

Veamos ahora cómo importar al script \( paquetes2.pl \) símbolos (variables) guardados en los archivos extermos Seq1.pm y Seq2.pm, cada uno de los cuales representa un módulo distinto (paquete guardado en un archivo).

• Echemos primero un vistazo a los archivos Seq1.pm y Seq2.pm.

# contents of Seq1.pm
package Seq1;
$secuencia = 'AAAAAAAA';

1;
-----------------------------------------------
# contents of Seq2.pm
package Seq2;
$secuencia = 'CCCCCCCC';

1;

Como ven, cada archivo inicia con la definición de un paquete con el mismo nombre que el archivo.

Además, los archivos \( *.pm \) terminan con una línea que contiene '1;' Es decir, terminan con un valor verdadero.

Programa paquetes2.pl - cont.

• Veamos ahora cómo se llaman los paquetes externos desde el script \( paquetes2.pl \).

La clave está en:

1. la llamada a la directiva use lib '/path/to/lib';
2. usar use Seq1; para llamar al paquete Seq1.pm (igual para Seq2.pm). ¡Nótese que no se usa la terminación '.*pm' en la llamada!

Estos paquetes estan guardados en los archivos Seq1.pm y Seq2.pm, y por tanto representan modulos (aunque todavia muy primitivos y poco funcionales, …)

Programa paquetes2.pl - llamada a un módulo - cont.

Ustedes tendrán que editar use lib '/path/to/lib'; para que apunte al directorio en el que han guardado los archivos *.pm.

#!/usr/bin/env perl
use strict;
use warnings;

my $host = `hostname`;

# programa: paquetes2.pl
# El concepto de paquetes como divisores del espacio de nombres para evitar colisiones entre símbolos.
# Importanmos símbolos desde dos paquetes externos usando la directiva 'use paquete;'.
# Estos paquetes estan guardados en los archivos Seq1.pm y Seq2.pm, y por tanto representan módulos
# (aunque todavia muy primitivos y poco funcionales, ...)

if ($host eq "alisio"){ 
   use lib qw(/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code);
}
else{
   # >>> EDITEN ESTA LINEA PARA APUNTAR AL PATH CORRECTO 
   #       EN EL QUE SE ENCUENTRAN LOS ARCHIVOS '*pm'    <<<
   use lib qw(/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code);
}

use Seq1;
use Seq2;

my $secuencia = 'MVLLIATG';

# imprimimos los valores de las 3 variables escalares secuencia
print "\n\tmain es: $secuencia"
     ."\n\tSeq1 es: $Seq1::secuencia"
     ."\n\tSeq2 es: $Seq2::secuencia\n\n";

Programa paquetes2.pl - llamada a un módulo - cont.

• de nuevo la llamada a \( print() \)

print "\n\tmain es: $secuencia"
     ."\n\tSeq1 es: $Seq1::secuencia"
     ."\n\tSeq2 es: $Seq2::secuencia\n\n";

• veamos la salida del progragma:

    main es: MVLLIATG
    Seq1 es: AAAAAAAA
    Seq2 es: CCCCCCCC

Diferencias entre las directivas \( require \) y \( use \)

El módulo es la unidad de reciclado de código en Perl. Para acceder a su funcionalidad, el módulo debe “exportar” funciones y variables, mientras que el script que quiere hacer uso de ella debe importar dichos símbolos.

Las directivas use y require cargan el módulo en tu programa, pero tienen un efecto o semántica diferentes.

• \( require \) carga los módulos durnte el tiempo de ejecución
• \( use \) carga los módulos durante el compilado del programa, mucho antes de que se ejecute

¿Qué significa ésto?

1. Un progrma que llama a un módulo con \( use \), pero éste no está disponible en el sistema, implica que el programa no va a empezar a correr, ya que no compila.
2. Además, \( use \) realiza in \( import \) implícito del paquete incluído en el módulo.

Por ello, una vez importado, el programa que llama al módulo no require ya usar el path completo para acceder al archivo *pm.

Asímismo, un símbolo importado con \( use \), ya no require del uso del nombre completamente cualificado del paquete.

Es decir, podemos usar las variables y subrutinas importadas como si fueran parte del paquete del script que las importa.

En esta sesión vamos a ir escribiendo progresivamente el módulo \( PhyloTools \) y el script \( run\_phyloTools.pl \).

Con ellos podremos hacer un análisis filogenético completo, incluyendo el alineamiento múltiple de secuencias, la estima filogenética bajo máxima verosimilitud, usando diferentes algoritmos y modelos, e incluso evaluar el ajuste de diversos modelos a los datos, usando tanto pruebas de razones de versomilitud (LRTs) como pruebas basadas en el contenido de información de Akaike (AIC).

Iremos construyendo dicho módulo y el scipt que lo llama paso a paso, explicando en el camino los componentes de un módulo.

Para empezar, descarguen y exploren los programas run_muscle.pl, \( run\_aligners.pl \), así como los módulos asociados \( PhyloTools1.pm \) y \( PhyloTools2.pm \)

La principal diferencia entre estos módulos radica en que el primero de ellos no exporta las subrutinas, mientras que el segundo sí lo hace.

Por ello, en \( run\_muscle.pl \) necesitamos llamar a la función de \( PhyloTools1.pm \) usando la declaración completa \( PhyloTools1::run\_muscle \).

En \( run\_aligners.pl \) podemos simplemente llamar a la subrutina \( run\_muscle \), ya que es exportada por \( PhyloTools2.pm \) gracias al uso del módulo Exporter, que viene con Perl.

Veamos los detalles:

1. Exploremos \( PhyloTools1.pm \) y corramos \( run\_muscle.pl \).

# Modulo de ejemplo del curso "Introduccion a la programacion en Perl para bioinformatica"
package PhyloTools1;
use strict;
use warnings;

use File::Basename;
use Carp;

my $muscle_bin = '/usr/bin/muscle';

my @binaries = ($muscle_bin);

foreach my $bin (@binaries)
{
  check_is_installed($bin);
}

#---------------------------------- SUBROUTINES ---------------------------------------

sub check_is_installed
{
   my @binaries = @_;
   foreach my $bin (@binaries)
   {
       if( ! -x "$bin")
       {
          croak "\n#ERROR! $bin is not installed. Please install or check path setting\n\tExiting now!\n";
       }
   }
}

\( PhyloTools2.pm \) y \( run\_aligners.pl \)

Como comentábamos arriba, \( PhyloTools1.pm \) no exporta las subrutinas, mientras que \( PhyloTools2.pm \) sí lo hace, gracias al uso del módulo core \( Exporter \), tal y como se ve en la cabecera del archivo \( PhyloTools2.pm \), mostrada seguidamente:

# Pablo Vinuesa, http://www.ccg.unam.mx/~vinuesa
# 2018 CCG/UNAM, Mexico, 
# Modulo de ejemplo del curso "Introduccion a la programacion en Perl para bioinformatica"
# el codigo y tutorial estan disponibles en: https://github.com/vinuesa/curso_perl4bioinfo

package PhyloTools2;
use strict;
use warnings;

our $VERSION = '0.2'; # v0.2 17 Nov. 2018
our(@ISA, @EXPORT);  # our($VAR1, @ARR1 ...) are package global symbols!
use Exporter;

@ISA    = qw(Exporter);                  # activamos el proceso de exportacion
@EXPORT = qw( run_muscle run_clustalo);  # exportamos automaticamente las subrutinas indicadas

Un módulo puede hacer disponibles variables y funciones a otros módulos o scripts mediante un proceso denominado exportación, un proceso que se inicia llamando a \( import \) al usar la directiva \( use\ lib \)

El módulo core \( Exporter \) implementa un método \( import \) estándar para módulos, que permite exportar símbolos de un módulo. \( Exporter \) utiliza las variables globales \( @EXPORT\_OK \) y \( @EXPORT \), que contienen una lista de símbolos para exportar, cuando sean requeridos.

Basados en la documentación de \( Exporter \), que podemos consultar con el comando \( perldoc\ Exporter \), vemos que:

package YourModule;
require Exporter;    # load the Exporter module from @INC
@ISA = qw(Exporter); # use YourModule implicitly calls YourModule->import()
@EXPORT_OK = qw(munge frobnicate);  # symbols to export on request by calling script

y en el script que llama al módulo escribimos:

use YourModule qw(frobnicate);  # import listed symbols; 
                             # implicitly calls YourModule->import()
frobnicate($left, $right)       # calls YourModule::frobnicate()

Veamos la cabecera de nuestro script \( run\_aligners.pl \), que contiene las directivas:

• use lib '/home/vinuesa/Cursos/Perl4bioInfo/my_code';
• use PhyloTools2;

#!/usr/bin/env perl

use strict;
use warnings;

use File::Basename;
use Cwd;
use Carp;

my $cwd = cwd();
my $hostname = '';

if ($hostname eq "alisio")
{ 
   # BEGIN { unshift @INC, '/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code' ;} 
   use lib qw(/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code) ;
}
else
{
   # >>> EDIT THIS LINE TO SET THE CORRECT PATH TO THE *pm FILES <<<
   # BEGIN { unshift @INC, '/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code' ;}
   use lib '/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code' ;
}

use PhyloTools2; # use PhyloTools2 qw(run_clustalo run_muscle); # <<< if using @EXPORT_OK()

my $VERSION = 0.2; # v0.2_16Oct2018
my $progname = basename($0);

print_help() unless @ARGV;

\( PhyloTools2.pm \) y \( run\_aligners.pl \)

Gracias a ello, en run_aligners.pl la llamada a la función se hace directamente:

...

if( $alignment_algorithm eq "muscle")
{
   my $mus_aln = run_muscle($fasta_file, 1); # <== llamada directas a la sub
}
else
{
  my $cluo_aln =  run_clustalo($fasta_file); # <== llamada directas a la sub
}

...

Corramos el script sin argumentos para ver la ayuda:

host=$(hostname)
if [ $host == "alisio" ]; then
   script_dir=/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code
else
   script_dir=$HOME/Dropbox/Cursos/perl4bioinfo/my_code
fi

$script_dir/run_aligners.pl

  USAGE:
    run_aligners.pl version 0.3 requires 2 arguments:

    - Required
        1. A fasta file name 

    - Optional
        2. an alignment algorithm name <clustalo|muscle>

  AIM:
    will run clustalo or muscle to align the input fasta file

  CWD:
    /home/vinuesa/Cursos/perl4bioinfo    

Ahora con argumentos:

host=$(hostname)
if [ $host == "alisio" ]; then
   script_dir=/home/vinuesa/Dropbox/Cursos/perl4bioinfo/my_code
else
   script_dir=$HOME/Dropbox/Cursos/perl4bioinfo/my_code
fi

$script_dir/run_aligners.pl GDP_12_prokEuc.faa

# &run_clustalo is running /usr/local/bin/clustalo -i GDP_12_prokEuc.faa -o GDP_12_prokEuc_cluoAln.faa --force
# &run_clustalo returned aligned outfile GDP_12_prokEuc_cluoAln.faa

Diferencias entre \( @EXPORT \) y \( @EXPORT\_OK \)

Es importante resaltar que la interfaz de usuario del módulo \( miModulo.pm \) puede:

1. Exportar variables por defecto con @EXPORT(&myFunc1 $VAR1);

   • que implica que si desde el script hacemos \( use\ miModulo \), podemos simplemente usar \( myFunc() \)

2. Exportar solo variables requeridas por el script que las llama con @EXPORT_OK(&myFunc1 $VAR1);

   • que implica que desde el script necesitamos importar explícitamente los símbolso que quermos con: use miModulo qw(&myFunc1 $VAR1);. Entonces podemos simplemente usar \( myFunc() \)

Diferencias entre declaraciones léxica \( my \) y global \( our \) de variables bajo el pragma \( use\ strict \)

1. la declaración léxica de una variable con my $variable la localiza a un “scope” o bloque
2. la declaración global de una variable con our $variable la hace global en el paquete y en el script que la llama

El método \( import \) heredado del paquete \( Exporter \) implícitamente llamado cuando el script principal usa \( use\ myModule \), busca variables globales en \( myModule \).

Por tanto, para poder exportar variables desde un módulo, éstas tienen que ser globales.

Dado que usamos el pragma \( use\ strict; \), necesitamos declarar a las variables globales con \( our \).

En resumen:

las variables globales a exportar de un módulo, que declaramos con \( our \), típicamente son:

• $VERSION
• \( @EXPORT \)
• \( @EXPORT\_OK \)

En conclusión:

es importante recordar que:

• Un paquete es simplemente una manera de agrupar variables y funciones relacionadas, pero sin conferir privacidad. Cualquier símbolo en un paquete es, por definición, global y accesible desde cualquier lugar. Los paquetes sólo agrupan, no esconden!.
• Al definir un módulo, es decir, un paquete en un archivo *pm, imponemos un “scope” o bloque mediado por el archivo, haciendo las variables de un módulo privadas de él.

Machote minimalista para escribir módulos que exportan funciones

package PackageName;
use strict;
use warnings;

our(@ISA, @EXPORT, $VERSION);
use Exporter;

$VERSION = '0.01';
@ISA     = qw(Exporter);
@EXPORT  = qw(sub1 sub2 …);

sub sub1
{
  # tu código
}

sub sub2
{
  # tu código
}

1; # recuerda que el módulo debe terminar retornando un valor verdadero!

Un módulo no está completo si no tiene documentación.

En general, toda la documentación de Perl se escribe en formato POD (Plain Old Documentation).

Se trata de un lenguaje de marcado que permite incluir bloques de documentación en un script o módulo de perl y su lectura desde la línea de comandos usando el programa \( perldoc \).

• cuando el intérprete \( perl \) encuentra un signo = al inicio de una linea, tal como =head1
• cambia al modo de procesamiento de POD, hasta que encuentra una línea que inicia con =cut

Para compilar el código, perl ignora estos bloques (como muchos programadores …)

1. Descarga \( PhyloTools3.pm \) y analiza su contenido, en particular los bloques \( POD \)
2. Corre \( perldoc\ PhyloTool3s.pm \)
3. Ejecuta \( podchecker\ PhyloTools3.pm \)

Aprende a usar \( POD \) estudiando \( perldoc\ perlpod \) (y \( perldoc\ perlpodspec \))

El arreglo \( @INC \)

¿Dónde busca Perl los módulos instalados en el sistema?

El compilar perl, se pasa una lista de directorios, que es sistema-específica, en la que se van a guardar los módulos estándar que se distribuyen con Perl, así como módulo instalados con poderes de superusuario.

La lista de esos directorios la exporta Perl al arreglo \( @INC \)

Para ver los detalles de configuración de la compilación del intérprete \( perl \), usa la siguiente línea, que incluye, al final de la copiosa salida, los directorios que se incluyen en el arreglo \( @INC \).

perl -V | tail -10; # imprime sólo las últimas líneas de la copiosa salida

  @INC:
    /etc/perl
    /usr/local/lib/x86_64-linux-gnu/perl/5.26.1
    /usr/local/share/perl/5.26.1
    /usr/lib/x86_64-linux-gnu/perl5/5.26
    /usr/share/perl5
    /usr/lib/x86_64-linux-gnu/perl/5.26
    /usr/share/perl/5.26
    /usr/local/lib/site_perl
    /usr/lib/x86_64-linux-gnu/perl-base

Cualquier módulo instalado en dichos directorios será accesible desde cualquier script que se corra en el sistema al invocar la directiva \( use\ module \)

Listar módulos de la distribución estándar

Perl viene con una gran cantidad de módulos pre-instalados y listos para ser usados.

Por ejemplo, la mayor parte de los > 66MB de la distribución v.5.14.2 corresponde a módulos, concretamente 652!

Puedes usar el módulo \( Module::CoreList \) para ver o contar los módulos que vienen con cada distribución.

¿Sabes qué versión de Perl estás corriendo?

perl -v; # imprime la versión de Perl

This is perl 5, version 26, subversion 1 (v5.26.1) built for x86_64-linux-gnu-thread-multi
(with 67 registered patches, see perl -V for more detail)

Copyright 1987-2017, Larry Wall

Perl may be copied only under the terms of either the Artistic License or the
GNU General Public License, which may be found in the Perl 5 source kit.

Complete documentation for Perl, including FAQ lists, should be found on
this system using "man perl" or "perldoc perl".  If you have access to the
Internet, point your browser at http://www.perl.org/, the Perl Home Page.

¿Sabes qué versión del kernel Linux estás corriendo y qué distribución?

uname -a; # imprime la versión del kernel de Linux
echo '-----------------------------------------------------------------------'
lsb_release -a # imprime detalles de la distribución de Linux

Linux alisio 4.15.0-70-generic #79-Ubuntu SMP Tue Nov 12 10:36:11 UTC 2019 x86_64 x86_64 x86_64 GNU/Linux
-----------------------------------------------------------------------
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 18.04.3 LTS
Release:    18.04
Codename:   bionic

Ejecuta este one-liner, para ver el número de \( módulos core \) en tu distirbución de Perl.

perl -MModule::CoreList -e 'print join "\n", Module::CoreList->find_modules(qr/.*/), "\n"' | wc -l

893

Es decir, Ubuntu 18.04.1 LTS corre Perl v5.26.1, con 893 módulos core.

Además de estos módulos de la distribución estándar, los programadores de Perl pueden usar cualquiera de los ~200,000 módulos existentes en CPAN (Comprehensive Perl Archive Network).

En conjunto, estos módulos permiten al programador de Perl escribir con gran facilidad programas para todo tipo de aplicaciones, desde programación Web y de bases de datos a estadística, graficado y bioinformática, por mencionar sólo algunas de las áreas relevantes para este curso.

Veremos cómo buscar módulos en CPAN y su instalación un poco más adelante. Ahora nos enfocaremos a explorar algunos módulos de la distribución estándar, aprendiendo a manejar \( perldoc \) y a usar módulos con interfaz funcional y dirigida a objetos.

Uso de \( perldoc \) para explorar la documentación de módulos de Perl

\( perldoc \) es un programa muy útil para explorar la copiosa documentación que viene con Perl y cada uno de sus módulos. Veamos algunos ejemplos:

Explora y aprende Perl con \( perldoc \), ejemplos de uso

• Explora la documentación de Perl: \( perldoc\ perl \)
• Ver el índice o tabla de contenidos de Perl: \( perldoc\ perltoc \)
• Ver la tabla de contenidos de perlfaq: \( perldoc\ perlfaq \)
• Ver documentación sobre sintaxis de Perl: \( perldoc\ perlsyn \)
• Ver documentación de variables reservadas: \( perldoc\ perlvar \)
• Ver la documentación de un módulo particular: \( perldoc\ nombreModulo \);
• Ver el código de un módulo: \( perldoc\ –m\ nombreModulo \)
• Ver la documentación de una función integrada: \( perldoc\ –f\ split \)
• Localiza dónde se instaló un módulo de Perl: \( perldoc\ –l nombreModulo \)
• Interroga las FAQs (Frequently Asqued Questions) : \( perldoc\ –q\ sort \)

Directorios y archivos

• \( Cwd \) Obtén el nombre de la ruta de acceso absoluta del directorio de trabajo actual
• \( File::Compare \) Compara el contenido de dos archivos
• \( File::Find \) Atraviesa un sistema de archivos para buscar un archivo (como find de UNIX)
• \( File::Glob \) Unix-like file globbing
• \( File::Spec \) Usa operaciones con archivos de manera portable

Procesado de cadenas, parseo de textos y listas

• \( Text::CSV\_XS \) Parseo de datos en formato CSV (separado por comas)
• \( Text::Wrap \) Formatea párrafos a un ancho determinado
• \( List::Util \) - A selection of general-utility list subroutines

Errores y avisos

• \( Carp \) provee las funciones \( carp \), \( croak \), \( confess \) y \( cluck \) para morir o avisar de errores desde la perspectiva del que invoca las llamadas, dando la traza de la pila

Interfaz de usuario: parseo de opciones y parámetros

• \( Getopt::Long \) Procesado de opciones de línea de comando en formato extendido –help|h
• \( Getopt::Std \) Procesado de switches de línea de comando de un solo carácter –f -fxz

Web, HTML, XML, bases de datos

• \( Catalyst \) Es el marco/framework de programación Web más popular de Perl
• \( CGI \) y \( CGI::Simple \) Acceso a y procesado de formas CGI y generación de HTML
• \( Net::FTP \) Conexión e interacción con sitios FTP
• \( LWP \) y \( LWP::Simple \) Las librerías del World-Wide Web para Perl
• \( WWW::Mechanize \) para extraer información de páginas Web
• \( HTML::Parser \) parseo de documentos HTML
• \( XML::Twig \) parseo de XML
• \( XML::Compile \) usa esquemas XML para traducir entre XML y Perl
• \( DBI \) el la interfaz abstrcta para bases de datos de Perl

y muchos, muchos más …

\( File::Basename \)

#!/usr/bin/env perl

use File::Basename;

$fullname = '/home/user/bin/my_perl_script.pl';
@suffixlist = qw(pl pm);

($name,$path,$suffix) = fileparse($fullname,@suffixlist);
$name = fileparse($fullname,@suffixlist);
$basename = basename($fullname,@suffixlist);
$dirname = dirname($fullname);

print "$name | $basename | $dirname | $suffix\n"

my_perl_script. | my_perl_script. | /home/user/bin | pl

• Exploren y corran ahora el script \( FileBasename.pl \)
• Exploren la documentación del módulo usando \( perldoc\ File::Basename \)

\( File::Basename \) - cont.

El módulo \( File::Basename \) detecta sobre qué tipo de sistema está corriendo, parseando correctamente los delimitadores de directorio, bien sean tipo UNIX ‘/’ o tipo Windows ‘\’, permitiendo escribir código portable, que corre bien en cualquier tipo de sistema operativo.

\( Getopt::Std \) - Process single-character switches with switch clustering

A medida que nuestro código se hace más complejo y necesitamos pasarle más parámetros, hacerlo como argumentos posicionales en la línea de comandos se hace cada vez más problemático.

Con el módulo base \( Getopt::Std \) podemos pasar fácilmente switches y parámetros como se muestra la sección de SYNOPSIS de la ayuda del módulo, la cual imprimimos con '\( perldoc\ Getopt::Std \)':

use Getopt::Std;

# -o & -i are boolean flags, -f takes an argument. Values in %opts
getopts('oif:', \%opts); 

Ejemplo de interfaz \( Getopt::Std \) para el script \( LRT\_calculator.pl \)

Veamos como ejemplo el script \( LRT\_calculator.pl \).

Si lo llamamos sin argumentos, imprime la ayuda:

. ./set_script_dir_by_hostname
$script_dir/LRT_calculator.pl

     LRT_calculator.pl vers. 0.1 usage [options]:
      -h  this help message
      -l  lnL0    (required)
      -L  lnL1    (required)
      -d  df      (required)
      -a  alpha   (optional, default:0.05)

  AIM:  
    LRT_calculator.pl uses as input the -lnL scores of a pair of nested H1 and H0 models  
       to calculate the corresponding LRT statistic [LRT = 2(lnL H1 - lnL H0)],
        and compute its significance at a user-defined alpha (critical value).

  USAGE EXAMPLE: 
    LRT_calculator.pl -l -19955.541 -L -19953.1245 -d 3
    LRT_calculator.pl -l -19955.541 -L -19950.1245 -d 3

Ejemplo de interfaz \( Getopt::Std \) para el script \( LRT\_calculator.pl \) - cont.

Veamos ahora las líneas clave del script que hacen uso de \( Getop::Std \)

my $a = 0.05; # set default alpha
my (%opts,$lnL0,$lnL1,$df);

getopts('l:L:d:a:h', \%opts);

if(($opts{'h'})||(scalar(keys(%opts))==0)) 
{ 
   print_help();
}

if(defined($opts{'h'})){ print_help(); }

if(defined($opts{'l'})){ $lnL0 = $opts{'l'}; }
else{ die "ERROR: no lnL score for the null hypothesis (lnL0) provided\n\n"; }

if(defined($opts{'L'})){ $lnL1 = $opts{'L'}; }
else{ die "ERROR: no lnL score for the alternate hypothesis (lnL1) provided\n\n"; }

if(defined($opts{'d'})){ $df = $opts{'d'}; }
else{ die "ERROR: no degrees of freedom (df) provided\n\n"; }

if(defined($opts{'a'})){ $a = $opts{'a'}; }

Ejemplo de interfaz \( Getopt::Std \) para el script \( LRT\_calculator.pl \) - cont.

Y para finalizar, corramos el script:

. ./set_script_dir_by_hostname
$script_dir/LRT_calculator.pl -l -19955.541 -L -19950.1245 -d 4

 LRT_calculator.pl v.0.1 running with the following parameters:

    lnL0=-19955.541, lnL1=-19950.1245, df=4, q=0.05


    LRT=10.8329999999987 | Chi-squared-crit-val (df=4; q=0.05) = 9.4877 | p-val=0.028506

>>> Significant ChiSq test! pchisq(0.028506) < alpha(0.05); LRT(10.8329999999987) > crit.val(9.4877) ==> reject the null model! 

Ejemplo de interfaz \( Getopt::Std \) para el script \( LRT\_calculator.pl \) - cont.

Noten que el órden de los parámetros no importa

. ./set_script_dir_by_hostname
$script_dir/LRT_calculator.pl -d 5 -L -19950.1245 -l -19955.541 -a 0.01

 LRT_calculator.pl v.0.1 running with the following parameters:

    lnL0=-19955.541, lnL1=-19950.1245, df=5, q=0.01


    LRT=10.8329999999987 | Chi-squared-crit-val (df=5; q=0.01) = 15.086 | p-val=0.054794

>>>Non-significant ChiSq test! pchisq(0.054794) > alpha(0.01); LRT(10.8329999999987) <= crit.val(15.086) ==> cannot reject the null model!

Vimos arriba que la distribución estándar de Perl trae muchos módulos.

Pero existen muchos más: CPAN es el repositorio público de módulos de Perl.

Lo accesan en https://www.cpan.org/.

Welcome to CPAN

The Comprehensive Perl Archive Network (CPAN) currently has 184,975 Perl modules
in 40,530 distributions, written by 13,870 authors, mirrored on 255 servers. 

...

¡No dejen de visitralo! Es un repositorio GIGANTESCO!

Tiene mucha información y código: 184,975 módulos!!!

Todos los lenguajes modernos tienen sus repositorios. Pero CPAN posiblemente sea el más grande.

Uno de los lemas de Perl es el de “no reinventar la rueda”.

Un programador avanzado de Perl, tiene un conocimiento amplio de CPAN, haciendo uso extensivo de los recursos allí disponibles.

Será realmente difícil que no encuentres un módulo en CPAN para ayudarte en resolver tus proyectos de programación.

Para buscar módulos en CPAN, lo mejor es usar uno de los buscadores especializados como Search-CPAN! Visita http://search.cpan.org/.

En este buscador pueden buscar por nombre de módulo o autor.

Una lista de “top-CPAN-authors” pueden verla aquí: http://thegestalt.org/simon/perl/wholecpan.html

Por ejemplo, hagan una búsqueda con el término Statistics para ver qué módulos encuentran.

Algunos módulos de estadística y matemáticas en CPAN

• \( Statistics::* \) Una gran gama de módulos
• \( Statistics::Descriptive \)- Module of basic descriptive statistical functions
• \( Statistics::Distributions \) - Perl module for calculating critical values and upper probabilities of common
• \( Statistics::Distribution::Generator \) - A way to compose complicated probability functions
• \( Statistics::R \) - Perl interface with the R statistical program
• \( Math::* \) - Una gran gama de módulos
• \( Math::BigInt \) - Arbitrary size integer/float math package

LRT_calculator.pl

Si bien ya habíamos corrido \( LRT\_calculator.pl \), no habiamos visto las líneas de llamada al módulo.

• Veamos primero la subrutina:

use Statistics::Distributions;

sub LRT_significance_calculator
{
   my ( $lnL0, $lnL1, $df, $a ) = @_;
   
   my $LRT = 2*( $lnL1 - $lnL0 );
   my $qchisq=Statistics::Distributions::chisqrdistr($df,$a);  
   my $pchisq=Statistics::Distributions::chisqrprob($df,$LRT);

   print STDOUT "\n\tLRT=$LRT | Chi-squared-crit-val (df=$df; q=$a) = $qchisq | p-val=$pchisq\n\n";
   
   return($LRT,$pchisq,$qchisq);
}

• y la llamada desde \( LRT\_calculator.pl \):

my ($LRT,$pchisq,$qchisq)=LRT_significance_calculator($lnL0, $lnL1, $df, $a);

if ($pchisq <= $a){ print ">>> Significant ChiSq test! pchisq($pchisq) < alpha($a); LRT($LRT) > crit.val($qchisq) ==> reject the null model! \n\n"}
else{ print ">>>Non-significant ChiSq test! pchisq($pchisq) > alpha($a); LRT($LRT) <= crit.val($qchisq) ==> cannot reject the null model!\n\n"}

Sencillo ¿veradad? Revisen arriba, en la sección “*2.4.2 Getopt::Std ejemplos de llamada al script \( LRT\_calculator.pl \).

Veamos ahora cómo insatalar módulos de CPAN. Tenemos varias opciones para ello. Sólo revisaremos algunas de ellas, en base a nuestros permisos y necesidades. Veamos algunos escenarios. Como ejemplo, veamos opciones para instalar el módulo \( Statistics::Distributions \) usado por el script \( LRT\_calculator.pl \) arriba mencionado.

Como adminstrador (root) queremos hacer instalaciones para todos los usuarios

Aquí tenemos nuevamente varias opiciones.

Usando el manejador de paquetes de Linux (\( apt \) en sistemas Debian)

1. busca cómo se llama el paquete con \( apt-cache\ search \)

apt-cache search Statistics::Distributions

1. instala como root

sudo apt install libstatistics-descriptive-perl

Si bien esta alternativa posiblemente no instalará la última versión de algunos módulos, sí permitirá que el sistema operativo automáticamente actualize el módulo, cuando una n nueva versión aparezca en el repositorio de Debian o Ubuntu, facilitando mucho el mantenimiento de múltiples módulos en el sistema. Es la vía que recomiendo.

Como usuario, queremos hacer una instalación en direcotrio bajo nuestro $HOME

Cuando trabajamos en un servidor, generalmente no tenemos permisos de administrador, por lo que no podemos usar la opción indicada arriba. En ese escenario, podemos:

1. Contactar al administrador para que nos instale el módulo
2. Generar un directorio bajo nuestro $HOME, donde instalamos los módulos.

Instalación manual con \( perl\ Makefile.PL \)

Veamos este segundo escenario. Vamos a instalar \( Statistics::Distributions \) que podemos encontrar aquí: https://metacpan.org/pod/Statistics::Distributions#AVAILABILITY

mkdir -p $HOME/CPAN/downloads
cd $HOME/CPAN/downloads

wget –c https://cpan.metacpan.org/authors/id/M/MI/MIKEK/Statistics-Distributions-1.02.tar.gz
tar –xvzf Statistics-Distributions-1.02.tar.gz
cd Statistics-Distributions-1.02

# Si no tienes poderes de administrador, instala el módulo en un directorio al que tengas acceso, por ejemplo: 
perl Makefile.PL INSTALL_BASE = $HOME/CPAN/
make
make test
make install

Como usuario, queremos hacer una instalación en direcotrio bajo nuestro $HOME

Instalación manual con \( perl\ Build.PL \)

Algunas distribuciones usan un archivo \( Build.PL \) para instalar módulos. El procedimiento es equivalente:

perl Buil.PL --install_base = $HOME/CPAN/
perl Build
Perl Build test
Perl Build install

Uso de paquete \( cpanm \) (CPANminus) para instalación de módulos

Una opción muy conveniente (mi favorita) para manejar módulos es mediante el módulo \( cpanm \) CPANminus, que puedes descargar desde aquí: https://cpan.metacpan.org/authors/id/M/MI/MIYAGAWA/App-cpanminus-1.7044.tar.gz. No es parte aún de la distribución estándar, pero es un módulo muy popular y útil.

En sistemas Linux/Debian puedes instalarlo en tu máquina con:

sudo apt install cpanm

Una vez instalado, explora la documentación con \( perldoc\ cpanm \)

Muestro la sección de SYNOPSIS

NAME
    cpanm - get, unpack build and install modules from CPAN

SYNOPSIS
      cpanm Test::More                                 # install Test::More
      cpanm MIYAGAWA/Plack-0.99_05.tar.gz              # full distribution path
      cpanm http://example.org/LDS/CGI.pm-3.20.tar.gz  # install from URL
      cpanm ~/dists/MyCompany-Enterprise-1.00.tar.gz   # install from a local file
      cpanm --interactive Task::Kensho                 # Configure interactively
      cpanm .                                          # install from local directory
      cpanm --installdeps .                            # install all the deps for the
                                                       # current directory
      cpanm -L extlib Plack                            # install Plack and all non-core 
                                                       # deps into extlib
      cpanm --mirror http://cpan.cpantesters.org/ DBI  # use the fast-syncing mirror
      cpanm --from https://cpan.metacpan.org/ Plack    # use only the HTTPS mirror

Uso de paquete \( cpanm \) (CPANminus) para instalación de módulos

Instalemos ahora el módulo \( Statistics::Descriptive \) usando \( cpanm \)

cpanm -i Statistics::Descriptive
cpanm --local-lib=~/perl5 local::lib && eval $(perl -I ~/perl5/lib/perl5/ -Mlocal::lib)

Si no corremos como superusuario (\( su \) o \( sudo \)) \( cpanm \) instala módulos en ~/perl5.

\( PERL5LIB \)

Si instalamos módulos como superusuario (\( su \) o \( sudo \)), éstos se instalan en uno de los directorios del arreglo especial \( @INC \), que contiene todas las rutas a directorios establecidos durante la compilación de \( perl \) para albergar módulos.

Podemos ver el contenido de @INC con el siguiente 1liner de Perl:

perl –e ‘print join “\n”, @INC, "\n"’

ó tecleando \( perl -V \), como ya vimos

Si no tenemos derechos de administrador y hemos instalado el módulo en un directorio en nuestro $HOME como ~/lib ó ~/CPAN/, este directorio no estará incluido en el arreglo \( @INC \).

Para que un script pueda encontrarlo, necesitamos indicar la ruta de acceso mediante una de las siguientes estrategias:

1. con la directiva \( use\ lib\ ‘/path/to/module/directory/’;\ use\ lib\ Path2Modulo; \)
2. mediante la declaración y exportación de la variable de ambiente \( PERL5LIB \)

# asumiendo que usas el shell Bash, añade esta línea a tu archivo ~/.bashrc o ~/.bash_profile
export PERL5LIB=$HOME/CPAN

Usando esta última estrategia, se garantiza que Perl siempre va a encontrar módulos instalados desde CPAN en un directorio personal. Es por tanto más portable y conveniente que la opción 1 (\( use\ lib \)).

Podemos listar y explorar los módulos instalados por el usuario (independientes de la distribución estándar) usando cualquiera de las sigientes opciones.

• \( instmodsh \) - un shell para explorar módulos instalados

perldoc instmodsh

Pueden explorar el código del módulo con un editor de texto como \( vi \), el editor por excelencia en Linux

vi $(which instmodsh)

• usando el módulo \( ExtUtils::Installed \)

Veamos su documentación:

perldoc ExtUtils::Installed

• Pueden usar el siguiente \( perl\ 1-liner \) llamando al módulo \( ExtUtils::Installed \) para ver la lista de módulos externos instalados

perl -MExtUtils::Installed -e 'my ($inst) = ExtUtils::Installed->new( skip_cwd => 1 ); my (@modules) = $inst->modules(); print join "\n",  @modules, "\n";'

• usando \( perldoc\ perllocal \)

  perldoc perllocal
  

• One-liners para encontrar todos los módulos instalados en el sistema

  
    find `perl -e 'print "@INC"' ` -name '*.pm' -print
    find `perl -e 'print "@INC"' ` -name '*.pm' -print | wc -l
    find `perl -e 'print "@INC"' ` -name '*.pm' -print | grep Statist

    # Perl 1liner para encontrar todos los módulos instalados en el sistema
    perl -MFile::Find=find -MFile::Spec::Functions -Tlwe \
        'find { wanted => sub{ print canonpath $_ if /\.pm\z/ }, no_chdir => 1 }, @INC'

Existen al menos dos paradigmas básicos de programación:

1. Programación por procedimientos (la que hemos visto/usado hasta ahora) http://en.wikipedia.org/wiki/Procedural_programming

2. Programación dirigida a objetos: las variables y código son encapsuladas en agrupaciones lógicas conocidas como clases. Las clases encapsulan descripciones generales de procesos u objetos de interés mediante la definición de interacciones entre objetos y sus métodos. http://en.wikipedia.org/wiki/Object-oriented_programming

Definiciones simplificadas de conceptos y términos de programación dirigida a objetos

Clases

Son módulos con una interfaz dirigida a objetos, y como tales, son paquetes

Objetos

Son un conjunto de datos con acciones y atributos asociados. Representan instancias particulares de una clase. Las acciones son funciones que un objeto puede hacer con los datos, las cuales vienen codificadas en los métodos del objeto.

Métodos

Son análogos a funciones, ya que tienen un nombre y reciben parámetros para operar sobre datos. Pero no son lo mismo ya que hay que llamarlos con una sintaxis particular que contiene el nombre de una clase (métodos de clases) o de un objeto (métodos de objetos).

Uso del módulo Statistics::Descriptive::Discrete con interfaz dirigida a objetos

• veamos el script statisticsDescriptiveDiscrete.pl

#!/usr/bin/env perl
use strict;
use warnings;
use File::Basename;
use Statistics::Descriptive::Discrete;

my $progname = basename($0); # statisticsDescriptiveDiscrete.pl
my $VERSION  = '0.01';
my $stats    = Statistics::Descriptive::Discrete->new();

print "# Please enter at least three integers between 0 and 100: e.g. 1 4 4 5 3 6 4\n";
chomp( my $data_string = <> );
my @digits = split /\s+/, $data_string;

die "# ERROR: need to enter at least three integers between 0 and 100: e.g. 1 4 4 5 3 6 4\n"
    if @digits < 3;

# make sure only integers were provided
my @non_digits = grep  { !/\d+/ } @digits;
die "# ERROR: need to enter at least three integers between 0 and 100: e.g. 1 4 4 5 3 6 4\n"
    if @non_digits ;

$stats->add_data(@digits); # añade los dígitos al objeto stats

print "count = ",              $stats->count(),              "\n";
print "uniq  = ",              $stats->uniq(),               "\n";
print "sum = ",                $stats->sum(),                "\n";
print "min = ",                $stats->min(),                "\n";
print "max = ",                $stats->max(),                "\n";
print "mean = ",               $stats->mean(),               "\n";
print "standard_deviation = ", $stats->standard_deviation(), "\n";
print "variance = ",           $stats->variance(),           "\n";
print "sample_range = ",       $stats->sample_range(),       "\n";
print "mode = ",               $stats->mode(),               "\n";
print "median = ",             $stats->median(),             "\n";

Uso del módulo Statistics::Descriptive::Discrete con interfaz dirigida a objetos

• la salida del script statisticsDescriptiveDiscrete.pl

# Please enter at least three integers between 0 and 100: e.g. 1 4 4 5 3 6 4
1 4 4 5 3 6 4     
count = 7
uniq  = 5
sum = 27
min = 1
max = 6
mean = 3.85714285714286
standard_deviation = 1.57359158493889
variance = 2.47619047619048
sample_range = 5
mode = 4
median = 4

Sintaxis básica para el uso de clases (módulos con interfaz dirigida a objetos):

1. Importación de la clase

   #use Class::Name
   use Statistics::Descriptive::Discrete;
   

2. Generación de una nueva instancia de la clase, es decir, generación de un objeto de la clase mediante un método de clase llamado constructor:

   # my $object = Class::Name->new(parameter1, parameter2 …); # con params.
   my $stats = Statistics::Descriptive::Discrete->new(); # sin parámetros
   

3. Uso de métodos de un objeto para realizar las acciones deseadas con los datos

   # $object->methodName(parameter1, parameter2);
   # $object->methodName(); # if no parameter required
   $stats->add_data(@digits);
   print "count = ", $stats->count(), "\n"; # llamada al método count del objeto $stats
   …
   

El projecto BioPerl de código fuente abierto es fruto de la colaboración entre bioinformáticos, biólogos y científicos del área de la computación con el objetivo de desarrollar una plataforma de código para resolver todo tipo de problemas de bioinformática.

El proyecto inició en 1995 y tuvo su lanzamiento oficial el 11 de Junio de 2002, por lo que claramente es el primer y mayor proyecto de código abierto para bioinformática actuamente existente.

BioPerl es un proyecto activo, financiado por la Open Bioinformatics Foundation, basado en módulos de Perl con el fin de facilitar la administración y manipulación de información relacionada con ciencias de la vida, en particular secuecias genómicas, análisis de genes y proteínas, e interacción con diversos tipos de bases de datos.

En términos prácticos, BioPerl es un gran repositorio de módulos de Perl para bioinformática con interfaz dirigida a objetos, llamados clases, que cuenta con una excelente documentación y tutoriales.

¿Porqué una interfaz dirigida a objetos?

Resulta que diseñar y abstraer el código en términos de clases y objetos es la técnica de programación más flexible, y en última instancia simple y escalable, para trabajar con datos de gran complejidad, como son los biológicos. Cada “dato”, como una secuencia o entrada de una base de datos, tiene una gran cantidad de atributos asociados, todos los cuales pueden ser modelados de manera consistente usando un sistema de clases o módulos con interfaz dirigida a objetos.

Instalación de BioPerl

En sistemas Linux basados en Debian, la instalación de bioperl se puede hacer fácilmente usando el manejador de paquetes \( apt \)

• Lista los paquetes relacionados al proyecto BioPerl

  apt-cache search bioperl
  

• Instalación en los directorios pre-definidos (\( @INC \)) del sistema

  sudo apt install bioperl
  

Nota: Para las prácticas de este tema, es suficiente que instalen el conjunto de módulos básicos

• Instalación en tu $HOME con CPANminus

  cpanm Bio::Perl
  

Si tienes otros sistemas operativos, sigue las instrucciones en la página de installation, que encontrarás aquí: https://bioperl.org/INSTALL.html

Clases, un tipo especial de módulo

Como se mencionó en la sección anterior, los módulos del proyecto BioPerl tiene una interfaz dirigida a objetos. Las clases son un tipo de módulo que incluye siempre, por definición, una subrutina '\( new \)', en español 'constructor', que genera las instancias de la clase al hacer la llamada my $object = clase->new().

Desde un punto de vista conceptual, un objeto o instancia de una clase es un agregado de variables (atributos), con identidad propia, que además incorpora una serie de subrutinas (métodos) para interaccionar con el mundo exterior o acceder a los contenidos del mismo.

En la práctica, una clase es un tipo de módulo que crea copias temporales de sí mismo, que duran mientras corre el programa que la invoca, cada una con su propia identidad. A cada instancia o copia de la clase le llamamos 'objeto'.

Un ejemplo en Biología podría ser la clase de las proteínas, que son polímeros de aminoácidos, pero cada una, cada objeto perteneciente a la clase, como la calmodulina, tiene su propia función, secuencia y nombre.

En este tema no podemos describir cómo se programan clases, nos limitaremos a aprender a usarlas. Es importante señalar que un módulo con interfaz dirigida a objetos, generalmente no utiliza el mecanismo de imporación-exportación descrito arriba para módulos tradicionales.

La interfaz dirigida a objetos se basa en el uso de constructores, destructores, métodos, herencia y sobrecarga de operadores.

Desde la óptica del usuario, lo relevante es entender que tenemos que generar objetos e interactuar con ellos mediante los métodos asociados.

Creación de un objeto de secuencias de la clase \( Bio::Seq \)

Dado que BioPerl cuenta con mucho módulos, tenemos que llamar al módulo adecuado para resolver nuestro problema. \( Bio::Seq \) es la clase principal para trabajar con secuencias en BioPerl. \( Bio::SeqIO \) puede leer y escribir secuencias biológicas en múltiples formatos, y por tanto tambien transformarlas de un formato a otro.

1. Llamada al módulo/clase \( Bio::Seq \)

#!/usr/bin/env perl 

use strict;
use warnings;

use Bio::Seq;

1. Generación de un objeto o instancia específica (copia) de la clase \( Bio::Seq \) guardado en la variable $\\( seq\_obj \). Esta variable es una referencia escalar.

# generación de un objeto
my $seq_obj = Bio::Seq->new(-seq => 'ATGATTGCACATTATAA',
                         -alphabet => 'dna' );

1. Para poder acceder a nuestra secuencia, tenemos que usar el método \( seq() \) del objeto $\\( seq\_obj \) para acceder a ella:

print $seq_obj->seq;

Uso de la clase \( Bio::Seq \)

#!/usr/bin/env perl 

use strict;
use warnings;

use Bio::Seq;

my $seq_obj = Bio::Seq->new(-seq      => 'ATGATTGCACATTATAA',
                            -alphabet => 'dna' );

print $seq_obj->seq;

ATGATTGCACATTATAA

• La variable $seq_obj es un el objeto sequencia, un objeto muy simple ya que sólo contiene la secuencia, con un solo atributo: alfabeto tipo dna (también hay rna y proteína).

• Los objetos de secuencia pueden ser creados manualmente, como hicimos arriba, pero son creados automáticamente en muchas operaciones de BioPerl, como cuando trabajamos con alineamientos o con entradas de bases de datos como Blast o GenBank.

• El objeto $seq_obj no sólo es un “contenedor” para la secuencia. Los objetos de la clase \( Bio::Seq \) tienen el método \( seq() \), que imprime su contenido.

• Además, podemos añadir más atributos al objeto $\\( seq\_obj \), como se muestra seguidamente.

Adición de nuevos atributos al objeto de la clase \( Bio::Seq \)

Típicamente una secuencia tiene varios atributos adicionales, como un identificador y una descripción con metadatos adicionales. Vamos a crear un nueva instancia de la clase \( Bio::Seq \), es decir, un nuevo objeto $\\( seq\_obj \), con esta información adicional.

$seq_obj = Bio::Seq->new(-seq        => "ATGATTGCACATTATAA",
                         -display_id => "XXX_0123",
                         -desc       => "Bichococcus imaginarius, genX",
                         -alphabet   => "dna" );

Imprimamos ahora la secuencia, con sus nuevos atributos, en formato FASTA

#!/usr/bin/env perl 

use strict;
use warnings;

use Bio::Seq;

my $seq_obj = Bio::Seq->new(-seq        => "ATGATTGCACATTATAA",
                            -display_id => "XXX_0123",
                            -desc       => "Bichococcus imaginarius, genX",
                            -alphabet   => "dna" );

# imprimimos la secuencia en formato FASTA
print ">", $seq_obj->id,  "|", $seq_obj->desc, "\n", $seq_obj->seq, "\n";

>XXX_0123|Bichococcus imaginarius, genX
ATGATTGCACATTATAA

Leer e imprimir secuencias con métodos de la clase \( Bio::SeqIO \)

Si queremos imprimir a un archivo el contenido del objeto $\\( seq\_obj \) que generamos arriba, podemos hacer uso de un nuevo objeto, de la clase \( Bio::SeqIO \), especializado en leer desde y escribir en archivos. De ahi su nombre SeqIO(), que alude a “sequence Input/Output”

use Bio::SeqIO;
...
my $seqio_obj = Bio::SeqIO->new(-file   => '>sequence.fasta', 
                                -format => 'fasta' );

Lo interesante es que los objetos de las clases \( Bio::Seq \) \( Bio::SeqIO \) pueden trabajar de manera coordinada, como se muestra en el siguiente ejemplo:

#!/usr/bin/env perl 
use strict;
use warnings;

# 1. lamada a las clases
use Bio::Seq; 
use Bio::SeqIO;

# 2. generamos una instancia de la clase Bio::Seq
my $seq_obj = Bio::Seq->new(-seq        => "ATGATTGCACATTATAA",
                            -display_id => "XXX_0123",
                            -desc       => "Bichococcus imaginarius, genX",
                            -alphabet   => "dna" );

# 3. imprimimos la secuencia en formato FASTA usando una instancia de Bio::SeqIO
# 3.1 generamos un nuevo objeto de la clase  Bio::SeqIO
my $seqio_obj = Bio::SeqIO->new(-file   => '>Bichococcus_imaginarius_genX.fasta', 
                                -format => 'fasta' );

Leer e imprimir secuencias con métodos de la clase \( Bio::SeqIO \)

# 3.2 usamos el método write_seq() del objeto $seqio_obj para escribir el contenido del
# objeto $seq_obj a disco.
$seqio_obj->write_seq($seq_obj);

[ -s Bichococcus_imaginarius_genX.fasta ] && cat Bichococcus_imaginarius_genX.fasta

>XXX_0123 Bichococcus imaginarius, genX
ATGATTGCACATTATAA

Cambiar el formato de las secuencias con \( Bio::SeqIO \)

Una gran ventaja de usar clases bien diseñadas, es que es muy fácil para el usuario cambiar atributos de un objeto. Por ejemplo, podemos pedirle a \( Bio::SeqIO \) que nos imprima la secuencia ahora en formato GenBank.

Para ello sólo necesitamos cambiar el atributo format del objeto $\\( seqio \), y usar un nuevo nombre de archivo de salida, como se muestra en el siguiente bloque.

my $seqio_obj = Bio::SeqIO->new(-file   => '>Bichococcus_imaginarius_genX.gbk', 
                                -format => 'genbank' );

LOCUS       XXX_0123                  17 bp    dna     linear   UNK 
DEFINITION  Bichococcus imaginarius, genX
ACCESSION   unknown
FEATURES             Location/Qualifiers
ORIGIN      
        1 atgattgcac attataa
//

¡Ta'güeno!, ¿verdad? ;)

Importar secuencias desde un archivo.

Para leer un archivo que tenemos en disco, usamos nuevamente la clase \( Bio::SeqIO \)

Como ejemplo, leamos el archivo GenBank que acabamos de generar e imprimamos en formato GenBank

#!/usr/bin/env perl 
use strict;
use warnings;

# 1. llamada a las clases
use Bio::SeqIO;

# 2. leamos la la secuencia en formato GenBank usando una instancia de Bio::SeqIO
# 2.1  generando un nuevo objeto de la clase  Bio::SeqIO
my $seqio_obj = Bio::SeqIO->new(-file   => 'Bichococcus_imaginarius_genX.gbk', 
                                -format => 'genbank' );

# 2.2 usamos el método next_seq() del objeto $seqio_obj leer el contenido del
# objeto $seqio_obj y asignarlo a un nuevo objeto $seq_obj.

# para leer una sola secuencia usamos:
# my $seq_obj = $seqio_obj->next_seq();

# si el objeto contiene varias, ponemos la llamada a next_seq() dentro de un bucle while
while ( my $seq_obj = $seqio_obj->next_seq ) {
    # print the sequence
    print ">", $seq_obj->id, "\n", $seq_obj->seq,"\n";
}

>XXX_0123
ATGATTGCACATTATAA

Script \( convert\_seqFormats\_SeqIO.pl \) para interconvertir formatos de secuencias

Con lo aprendido en secciones anteriores, debería ser fácil escribir un interconvertidor de formatos de secuencias.

#!/usr/bin/env perl
use strict;
use warnings;
use File::Basename;
use Bio::SeqIO; 

my $progname = basename($0);
my $VERSION = 0.2; # November 21st, 2013

# get command-line arguments, or die with a usage statement
my $usage =<<"EOF";

USAGE:
$progname v.$VERSION infile infileformat outfile outfileformat

 accepted formats are:
   Fasta       FASTA format
   EMBL        EMBL format
   GenBank     GenBank format
   swiss       Swissprot format
   PIR         Protein Information Resource format
   GCG         GCG format
   raw         Raw format (one sequence per line, no ID)
   ace         ACeDB sequence format
   game        GAME XML format
   phd         phred output
   qual        Quality values (get a sequence of quality scores)
   Fastq       Fastq format
   SCF         SCF tracefile format
   ABI         ABI tracefile format
   ALF         ALF tracefile format
   CTF         CTF tracefile format
   ZTR         ZTR tracefile format
   PLN         Staden plain tracefile format
   EXP         Staden tagged experiment tracefile format
EOF

Script \( convert\_seqFormats\_SeqIO.pl \) para interconvertir formatos de secuencias

my $infile        = shift or die $usage;
my $infileformat  = shift or die $usage;
my $outfile       = shift or die $usage;
my $outfileformat = shift or die $usage;

# create one SeqIO object to read in,and another to write out
my $in = Bio::SeqIO->new(
                             -file   => "<$infile",
                             -format => $infileformat,
                             );
my $out = Bio::SeqIO->new(
                              -file   => ">$outfile",
                              -format => $outfileformat,
                              );
 
# write each entry in the input file to the output file
while (my $seq = $in->next_seq) {
       $out->write_seq($seq);
}

Script \( convert\_seqFormats\_SeqIO.pl \) para interconvertir formatos de secuencias

Una llamada al script, para convertir el archivo Bichococcus_imaginarius_genX.gbk en formato GenBank a formato EMBL, sería así:

./convert_seqFormats_SeqIO.pl Bichococcus_imaginarius_genX.gbk genbank Bichococcus_imaginarius_genX.embl EMBL

Y podemos ver el resultado así

[ -s Bichococcus_imaginarius_genX.embl ] && cat 'Bichococcus_imaginarius_genX.embl'

ID   XXX_0123   standard; dna; UNK; 17 BP.
XX
AC   unknown;
XX
DE   Bichococcus imaginarius, genX
XX
FH   Key             Location/Qualifiers
FH
XX
SQ   Sequence 17 BP; 7 A; 2 C; 2 G; 6 T; 0 other;
     atgattgcac attataa                                                       17
//

¿Dónde encuentro más información para aprender más \( Bio::Perl \)?

Para más ejemplos de uso de la clase \( Bio::SeqIO \), vean el siguiente HOWTO en la página de bioperl.org https://bioperl.org/howtos/

En la página de BioPerl HOWTOs encontrarás muchos ejemplos más, muy bien escritos. Para empezar, recomiendo estudiar los siguientes, en este órden

• Beginners HOWTO
• SeqIO HOWTO
• SearchIO HOWTO
• EUtilities Cookbook HOWTO
• Features and annotations HOWTO

Escribir \( run\_PhyloTools.pl \) que llama a PhyloTools.pm

En el repositorio está el script \( run\_phylo\_pipeline1.pl \) que han de modifcar para que:

1. Tenga una interfaz de usuario con Getopts::Std
2. Llame a archivos de secuencia de DNA o Proteína en formao FASTA
3. Haga alineaminetos múltiples con muscle o clustalo, y escriba los alineamientos en un formato especificado por el usuraio
4. Corra FastTree o Phyml para inferir filogenias de máxima verosimilitud, a partir de los alineamientos. Tomar en cuenta que FastTree requiere alineamientos en formato FASTA y PhyML en formato PHYLIP
5. En el caso de elegir PhyML, que el usuario le pueda pasar el modelos de sustitución o matriz empírica a usar para estimar la filogenia de alineamientos de nt|aa
6. Pueda hacer LRTs entre parejas de hipótesis anidadas, para seleccionar el mejor modelo
7. Llame a PhyloTools.pm, que tiene las subrutinas ya pre-programadas para los requisitos arriba listados

• Puden usar el script \( run\_phylo\_pipeline1.pl \) como punto de arranque, y modifíquenlo acorde a sus necesidades, o para hacer lo que arriba se sugiere.

Escribir \( run\_PhyloTools.pl \) que llama a PhyloTools.pm

• Va un ejemplo de corrida \( run\_phylo\_pipeline1.pl \)

. ./set_script_dir_by_hostname
$script_dir/run_phylo_pipeline1.pl

  USAGE:
    run_phylo_pipeline1.pl version 0.2 requires 2 arguments:

    - Required
        1. A fasta file name 
    2. Sequence type: <aa|nt>

    - Optional
        3. an alignment algorithm name <clustalo|muscle>; default: clustalo

  AIM:
    will run clustalo or muscle to align the input fasta file 
         and FastTree to infer phylogeny

  CWD:
    /home/vinuesa/Cursos/perl4bioinfo    

Escribir \( run\_PhyloTools.pl \) que llama a PhyloTools.pm

• Va un ejemplo de corrida \( run\_phylo\_pipeline1.pl \) con argumentos

. ./set_script_dir_by_hostname
$script_dir/run_phylo_pipeline1.pl GDP_12_prokEuc.faa aa clustalo

# sub run_clustalo is running /usr/local/bin/clustalo -i GDP_12_prokEuc.faa -o GDP_12_prokEuc_cluoAln.faa --force
# run_clustalo() returned aligned outfile GDP_12_prokEuc_cluoAln.faa
# sub run_FastTree is running: /usr/local/bin/FastTree -lg -gamma -quiet < GDP_12_prokEuc_cluoAln.faa > GDP_12_prokEuc_cluoAln_FTLGG.ph
# run_clustalo() returned aligned outfile GDP_12_prokEuc_cluoAln_FTLGG.ph

• pueden visualizar y editar el alineamiento GDP_12_prokEuc_cluoAln.faa con \( seaview \)
• pueden visualizar y editar el árbol GDP_12_prokEuc_cluoAln_FTLGG.ph con \( figtree \)

Es todo, les deseamos mucho éxito en su futuro como programadores de Perl!