Qué ha hecho la Bioinformática por nosotros?

Aunque esto quedaría mejor en el blog de mi colega José Maria, os paso este enlace en inglés que me ha llegado por la lista de correo del Protein Data Bank, donde puedes proponer y elegir la contribución más importante de la Bioinformática a la ciencia. Alguien ha propuesto ya la comparación entre el genoma humano y el del chimpancé, por si estabas pensando en esto, saludos.

Clusters de proteínas ortólogas mediante el algoritmo COGs

Tras leer el artículo que sale hoy en Nature sobre la expansión del espacio de secuencias de proteínas, que se asemeja a la expansión del Universo tal como la imaginaba Hubble, he vuelto a encontrarme con el algoritmo COGs (Clusters of Orthologous Groups), originalmente creado por Tatusov en 1997 y actualizado por David Kristensen hace unas semanas. En esencia este algoritmo permite definir grupos de secuencias (potencialmente) ortólogas entre 3 o más genomas a partir de los resultados de BLAST de todos contra todos. Para ello utiliza el concepto de los 'bidirectional best hits', que se dibuja como flechas bidireccionales en esta figura, tomada del paper de Kristensen:

El código que publico hoy permite calcular COGs a partir de un directorio con archivos FASTA con secuencias de proteínas, utilizando los binarios del algoritmo COGtriangles, que se pueden obtener de SourceForge, y los de BLAST, que podemos descargar del NCBI.

 #!/usr/bin/perl -w  

 # prog5.pl    
 # Script que ejecuta el algoritmo COGtriangle sobre los archivos FASTA contenidos  
 # en el directorio -d, basado en el código original de David Kristensen (Jun 2010).  
 # El programa asume que cada archivo .faa contiene CDS de genomas distintos   
 # y puede agrupar la misma secuencia en mas de un cluster, cuando sea multidominio.  
   
 use strict 'subs';  
 use Getopt::Std;  
 use File::Basename;  
   
 ## 0) global variables, including binaries to be installed   
 ##  and parameters that affect the whole algorithm  
 # URL: ftp://ftp.ncbi.nih.gov/blast/executables/  
 my $BLASTPATH  = '/home/whatever/blast-2.2.18/';  
 my $FORMATDBEXE = $BLASTPATH . 'bin/formatdb';  
 my $BLASTPGPEXE = $BLASTPATH . 'bin/blastpgp -I T -m 9 -b 1000 -v 1000 -z 100000000';  
 # URL: http://sourceforge.net/projects/cogtriangles/files/  
 my $COGPATH = '/home/soft/COGsoft/';  
 my $MAKEHASHEXE = $COGPATH . "COGmakehash/COGmakehash -s=',' -n=1";  
 my $READBLASTEXE= $COGPATH . 'COGreadblast/COGreadblast -e=10 -q=2 -t=2';  
 my $LSEEXE   = $COGPATH . 'COGlse/COGlse';  
 my $COGTRIANGLESEXE=   
    $COGPATH . "COGtriangles/COGtriangles -t=0.5 -e=10 -n='cluster' -s=1";  
 # the default COGtriangle algorithm requires two BLAST searches  
 my $TWOBLASTRUNS=1; # set to 0 to perform only a filtered BLAST search  
 # extension of valid files in -d directories  
 my $FAAEXT   = '.faa';   
   
   
 ## 1) process command-line arguments ##############################  
 my ($INP_FASTAdir,$INP_skipBLAST,$INP_BLASToptions,$newDIR,$id) = ('',0,'');  
 my ($n_of_sequences,$pwd,$genome,$clusterid,%sequences,%genomeused,%COGs) = (0);  
 my ($allfaafile,$p2ofilename,$lsefilename,$lseoutfilename,$coglogfilename) =   
    ('all.faa','all.p2o.csv','all.lsejob.csv','all.lse.csv','all.cog.clusters.log');  
 my ($cogedgesfilename,$edgesfilename) = ('cog-edges.txt','all-edges.txt');  
   
 getopts('hsb:d:', \%opts);  
 if(($opts{'h'})||(scalar(keys(%opts))==0))   
 {   
    print "usage: $0 [options]\n";  
    print " -h this message\n";  
    print " -d input directory with FASTA .faa files (required, example -d genomes)\n";  
    print " -s skip BLAST and re-use previous searches (optional)\n";  
    die  " -b extra blastpgp options (optional, example -b '-a 2' to use 2 CPUs)\n";     
 }  
   
 if(!defined($opts{'d'})   
    || !-d $opts{'d'}){ die "# $0 : need a valid FASTA directory, exit\n"; }  
 else  
 {   
    $INP_FASTAdir = $opts{'d'};   
    if(substr($INP_FASTAdir,length($INP_FASTAdir)-1,1) eq '/'){ chop $INP_FASTAdir }  
    $pwd = `pwd`; chomp $pwd; $pwd .= '/';  
     $newDIR = $pwd.basename($INP_FASTAdir)."_COGs";  
    mkdir($newDIR) if(!-d $newDIR);  
 }  
   
 if(defined($opts{'s'})){ $INP_skipBLAST = 1 }  
 if(defined($opts{'b'}) && !$INP_skipBLAST){ $INP_BLASToptions = $opts{'b'} }  
   
 print "# results_directory=$newDIR\n\n";  
   
   
 ## 2) concatenate $FAAEXT files and create csv file with protein2organism indexing   
 opendir(FAADIR,$INP_FASTAdir) || die "# $0 : cannot list $INP_FASTAdir, exit\n";  
 @faafiles = grep {/$FAAEXT/} readdir(FAADIR);  
 close(FAADIR);  
 die "# $0 : need at least two $FAAEXT files in $INP_FASTAdir\n"   
     if(scalar(@faafiles)<2);   
   
 open(ALLFAA,">$newDIR/$allfaafile")   
     || die "# $0 : cannot create $newDIR/$allfaafile\n";  
 open(PROTORG,">$newDIR/$p2ofilename")   
     || die "# $0 : cannot create $newDIR/$p2ofilename, exit\n";  
 foreach my $faafile (@faafiles)  
 {  
    print "# reading $INP_FASTAdir/$faafile\n";  
      
    my $rhash_FASTA = read_FASTA_sequence("$INP_FASTAdir/$faafile");  
      
    foreach my $seq (sort {$a<=>$b} (keys(%$rhash_FASTA)))  
   {  
      $n_of_sequences++;  
         
       print PROTORG "$n_of_sequences,$faafile\n"; # csv file  
       #concatenated FASTA with numbers as headers  
       print ALLFAA ">$n_of_sequences\n$rhash_FASTA->{$seq}{'SEQ'}\n";   
         
       #store number - real header equivalence, could be done with BerkeleyDB   
      $sequence_data{$n_of_sequences} = $rhash_FASTA->{$seq}{'NAME'};  
     $sequence_data{$n_of_sequences.'PROT'} = $rhash_FASTA->{$seq}{'SEQ'};  
   }  
      
    $genomeused{$faafile} = 1;  
 }  
 close(PROTORG);  
 close(ALLFAA);  
   
   
 ## 3) run BLASTs searches #########################################  
 if(!$INP_skipBLAST)  
 {  
    print "# Formatting BLAST database...\n";  
    system("$FORMATDBEXE -i $newDIR/$allfaafile -o T");  
      
    if(-s 'formatdb.log'){ system("mv -f formatdb.log $newDIR") }  
   
    if($TWOBLASTRUNS)  
    {  
       print "# Running unfiltered BLAST...\n";  
       mkdir($newDIR.'/BLASTno');  
       system("$BLASTPGPEXE -d $newDIR/$allfaafile -i $newDIR/$allfaafile -F F -t F ".  
          "-o $newDIR/BLASTno/all.tab $INP_BLASToptions 2> " .  
           "$newDIR/BLASTno/blastlog_unfiltered");  
    }  
      
    print "# Running filtered BLAST...\n"; # see PMID: 20439257  
    mkdir($newDIR.'/BLASTff');  
    system("$BLASTPGPEXE -d $newDIR/$allfaafile -i $newDIR/$allfaafile -F T -t T ".  
          "-o $newDIR/BLASTff/all.tab $INP_BLASToptions 2> " .   
             "$newDIR/BLASTff/blastlog_unfiltered");  
 }  
   
   
 ## 4) calculate COG triangles ######################################  
 print "\n# making COG hash...\n";  
 mkdir($newDIR.'/BLASTconv');  
 system("$MAKEHASHEXE -i=$newDIR/$p2ofilename -o=$newDIR/BLASTconv");  
   
 print "# reading BLAST files...\n";  
 if($TWOBLASTRUNS)  
 {  
    system("$READBLASTEXE -d=$newDIR/BLASTconv -u=$newDIR/BLASTno " .    
             " -f=$newDIR/BLASTff -s=$newDIR/BLASTno");   
 }  
 else{ system("$READBLASTEXE -d=$newDIR/BLASTconv -u=$newDIR/BLASTff ".   
             "-f=$newDIR/BLASTff -s=$newDIR/BLASTff"); }  
   
 print "# checking Lineage-specific expansions...\n";  
 open (LSEIN, ">$newDIR/$lsefilename")   
    || die "# $0 : cannot create $newDIR/$lsefilename, exit\n";  
 foreach $genome (sort(keys(%genomeused)))  
 {  
    foreach my $genome2 (sort(keys(%genomeused)))  
    {  
       next if($genome eq $genome2);  
       print LSEIN "$genome,$genome2\n";  
    }  
 }  
 close LSEIN;  
   
 mkdir($newDIR.'/tmp');  
 system("$LSEEXE -t=$newDIR/tmp -d=$newDIR/BLASTconv -j=$newDIR/$lsefilename ".  
        " -p=$newDIR/$p2ofilename -o=$newDIR/$lseoutfilename > /dev/null");  
   
 print "# making COGs...\n"; # genera all-edges.txt , cog-edges.txt  
 system("$COGTRIANGLESEXE -i=$newDIR/BLASTconv -q=$newDIR/$p2ofilename " .  
           " -l=$newDIR/$lseoutfilename -o=$newDIR/$coglogfilename");  
         
 # $COGTRIANGLESEXE puts them in cwd, Jun2010!  
 if(-e $cogedgesfilename){ system("mv -f $cogedgesfilename $newDIR") }   
 if(-e $edgesfilename){ system("mv -f $edgesfilename $newDIR") }  
   
   
 ## 5) group sequences sharing COG #################################  
 open(COGS,"$newDIR/$cogedgesfilename") ||   
    die "# $0 : cannot read $newDIR/$cogedgesfilename, the program failed\n";  
 while(<COGS>)  
 {  
    #cluster00332,912,Buch_aph_Bp.faa,405,Buch_aphid_Sg.faa  
    my @data = split(/,/,$_);  
    ($clusterid,$id,$genome) = @data[0,1,2];  
    push(@{$COGs{$clusterid}{'genomes'}},$genome)   
         if(!grep(/^$genome$/,@{$COGs{$clusterid}{'genomes'}}));  
    push(@{$COGs{$clusterid}{'ids'}},$id)   
         if(!grep(/^$id$/,@{$COGs{$clusterid}{'ids'}}));  
 }  
 close(COGS);  
   
 print "\n# cluster list:\n";  
 mkdir($newDIR.'/clusters');  
 foreach $clusterid (sort {substr($a,7)<=>substr($b,7)} (keys(%COGs)))  
 {  
    print ">$newDIR/clusters/$clusterid.faa sequences=".  
          scalar(@{$COGs{$clusterid}{'ids'}}).  
       " genomes=".scalar(@{$COGs{$clusterid}{'genomes'}})."\n";  
   
    open(COGFAA,">$newDIR/clusters/$clusterid.faa") ||   
       die "# $0 : cannot create $newDIR/cogs/$clusterid.faa, exit\n";  
    foreach $id (@{$COGs{$clusterid}{'ids'}})  
    {  
       print COGFAA ">$sequence_data{$id} ($clusterid)\n$sequence_data{$id.'PROT'}\n";  
    }  
    close(COGFAA);  
 }  
   
 print "\n# total COGs = ".scalar(keys(%COGs))." in folder $newDIR/clusters\n";  
   
 exit(0);  
   
 ################################################################  
 sub read_FASTA_sequence  
 {  
    my ($infile) = @_;  
    my ($n_of_sequences,%FASTA) = (0);  
      
   open(FASTA,$infile) || die "# read_FASTA_sequence: cannot read $infile $!:\n";  
   while(<FASTA>)  
   {  
      if(/^\>(.*?)\n/)  
     {  
       $n_of_sequences++;  
       $FASTA{$n_of_sequences}{'NAME'} = $1;  
       }  
     else  
     {  
        $FASTA{$n_of_sequences}{'SEQ'} .= $_;  
       $FASTA{$n_of_sequences}{'SEQ'} =~ s/[\s|\n|\d]//g;  
       }  
    }     
   close(FASTA);  
          
    return(\%FASTA);  
 }  
   

Algunos comandos de R útiles en ciencia e investigación

Inspirado en:

http://www.personality-project.net/r/r.commands.html

http://www.statmethods.net/graphs/scatterplot.html

http://www.ats.ucla.edu/stat/R/notes/

Básicos

• Ver el directorio de trabajo:
  

getwd()

• Cambiar el directorio de trabajo:
  

setwd("path")

• Salir de R:
  

q()

• Obtener ayuda de un comando (las comillas son importantes):
  

help("command")

?"command"

• Ver los objetos y variables del espacio de trabajo:
  

ls()

• Borrar un objeto del espacio de trabajo:
  

rm(x)

rm(list=ls()) Borrar todos los objetos

• Guardar el espacio de trabajo:
  

save.image() # en el archivo .RData en el directorio de trabajo

save(object list,file="myfile.RData") # los objetos deseados en el archivo elegido

• Cargar un espacio de trabajo:
  

load("myfile.RData")

• Salir de R:
  

q()

Vectores y matrices

• Crear un vector:
  

x=c(1,2,4,8,16)

y=c(1:10)

z=c(rnorm(n))

• Operaciones con vectores:
  

x=c(x)+n # suma 'n' a cada elemento

z=c(x,y) # combina 2 vectores creando un vector lineal

z=cbind(x,y) # combina 2 vectores creando una matriz de 2 columnas

z=rbind(x,y) # combina 2 vectores creando una matriz de 2 filas

replace(x,x==0,NA) # sustituye los ceros por NA

• Operaciones con matrices:
  

mat[4,2] # muestra la 4ª fila y la 2ª columna

mat[3,] # muestra la 3ª fila

mat[,2] # muestra la 2ª columna

mat[,-3] # borra la 3ª columna

mat[-2,] # borra la 2ª fila

mat[1:3,3:5] # muestra las filas 1 a 3 y las columnas 3 a 5

Frames

• Crear un frame a partir de vectores (cada vector una columna del frame):
  

data.fr=data.frame(x1,x2,x3 …)

• Crear un frame a partir de una matriz:
  

as.data.frame(mat)

• Conocer si un objeto es un frame:
  

is.data.frame(mat)

• Convertir un frame en matriz:
  

as.matrix(data.frame)

• Leer un fichero de datos separados por tabuladores con los nombres de sus columnas (1ª línea) y guardarlo en un frame:
  

data.fr<-read.table("ruta fichero",header=T,sep="\t")

• Obtener los nombres de las columnas de un frame:
  

names(data.fr)

• Obtener los nombres de las filas de un frame:
  

row.names(data.fr)

• Obtener el número de filas de un frame:
  

nrow(data.fr)

• Ver las primeras filas del frame:
  

head(data.fr, n=10) # las 10 primeras líneas

• Crear vectores con el nombre de cada columna del frame (y sus datos)
  

attach(data.fr)

detach(data.fr) # borrar los vectores

• Añadir una fila a un frame:
  

data.fr <- rbind(data.fr,data.frame(colA=1,colB="abc",colC=rnorm(1)))

• Añadir una columna a un frame:
  

data.fr$colC <- data.fr$colA + 5 * data.fr$colB

• Filtrar un frame para obtener un subconjunto de datos y guardarlo en otro frame:
  

subset.frame<-subset(data.fr,colA>=5 & !is.na(colB) | ColC='V')

subset.frame<-subset(data.fr,complete.cases(data.fr)) # Obtener únicamente las filas con datos en todas sus columnas

• Ordenar un frame:
  

data.fr[order(data.fr$colB),] # Ordena el frame con los valores de O la columna B

data.fr[rev(order(data.fr$colB)),] # Orden inverso

Estadística

• Calcular varios parámetros estadísticos de un conjunto de datos:
  

max(x)
min(x)

mean(x)

median(x)

sum(x)

var(x) # produces the variance covariance matrix

sd(x) # standard deviation

mad(x) # median absolute deviation

fivenum(x) # Tukey fivenumbers min, lowerhinge, median, upper hinge, max

Table(x) # Matriz de frecuencias

scale(x,scale=T) # centers around the mean and scales by the sd

cumsum(x) # cumulative sum

cumprod(x)

cummax(x)

cummin(x)

rev(x) # reverse the order of values in x

• Matriz de correlación:
  

cor(x,y,use="pair") #correlation matrix for pairwise complete data, use="complete" for complete cases

• Test de correlación:
  

cor.test(x,y,method=c('pearson'))

• ANOVA (Análisis de la varianza):
  

aov.1 <- aov(colA~colB,data.fr) # one way analysis of variance (colA contains values and colB contains classes)

aov.2 = aov(colA~colB*colC,data.fr) #do a two way analysis of variance

summary(aov.1) #show the summary table

print(model.tables(aov.1,"means"),digits=3) #report the means and the number of values/class

boxplot(colA~colB,data.fr) #graphical summary appears in graphics window

• Regresión lineal:
  

fit<-lm(colB~colA,data.fr) #basic linear model where x and y can be matrices (see plot.lm for plotting options)

summary(fit)

plot(data.fr$colA,data.fr$colB)

abline(fit)

• Student's t-Test:
  

t.test(x,g)

pairwise.t.test(x,g)

power.anova.test(groups = NULL, n = NULL, between.var = NULL,

within.var = NULL, sig.level = 0.05, power = NULL)

power.t.test(n = NULL, delta = NULL, sd = 1, sig.level = 0.05,

power = NULL, type = c("two.sample", "one.sample", "paired"),

alternative = c("two.sided", "one.sided"),strict = FALSE)

Gráficas

• Lista de colores disponibles en R:
  

colors()

• Parámetros gráficos:
  

help(par)

• Gráfico sencillo:
  

plot(x,y,type="p") # types: “p”: points, “l”: line, …

plot(x,y,col="red",lwd=2,type="l")

• Dibujar un punto:
  

points(2,5,pch=16,col="blue",cex=3) # Dibuja el punto (2,5) en azul y con un tamaño relativo de 3

• Etiquetar un punto:
  

text(2,5,label="Punto",pos=4,offset=-0.5,font=3,cex=0.75)

• Gráfico de puntos:
  

plot(data.fr$colA, data.fr$colB,

main="Column A vs. Column B", xlab="Column A", ylab="Column B",

cex.main=2, cex.lab=1.5, pch=1)

• Superponer gráfica:
  

lines(1-spec, sens, type='b', col=39, pch=7, lty=2)

plot(data.fr$colA, data.fr$colB,

main="Column A vs. Column C",

cex.main=2, cex.lab=1.5, pch=1, add=TRUE)

• Ejemplo de gráfica pequeña dentro de otra mayor:
  

year <- 1900:2000

dollars <- (year-1899)^2

plot.within.aplot <- function()

{

    par(pin=c(1.5,1.5)) ### set the plot region size in inches

    par(mai=c(2.8,1.2,1,2.6)) ### set the 4 margin sizes in inches

    par(ps=6) ### set the font size to 6 point

    plot(year,dollars,col="blue",type="l",xlab="",ylab="") ### no axis labels

    par(new=T) ### the next plot will plot over the previous one

    par(ps=12) ### set the font size to 12 point

    par(pin=c(4.14,3.57)) ### set the plot region size back to the default

    par(mai=c(.956,.769,.769,.394)) ### set the margins back to the default

    plot(year,dollars,col="green",type="l") ### plot

}

plot.within.aplot()

Salida de datos

• Opciones de ventanas gráficas:
  

X11() # abrir ventana en Linux

windows() # abrir ventana en Windows

dev.list() # mostrar ventanas abiertas

dev.cur() # mostrar ventana actual de trabajo

dev.off() # cerrar ventana actual de trabajo

dev.set(3) # seleccionar la ventana 3 como ventana de trabajo

dev.off(3) # cerrar la ventana 3

graphics.off() # cerrar todas las ventanas

• Abrir fichero PDF para salida de datos:
  

pdf(file=”outfile.pdf”)

Tutorial para escribir trabajos académicos y tesis doctorales

A continuación se describirán técnicas útiles y recomendadas para escribir trabajos académicos y tesis doctorales.

1. Aprender a usar un buen editor de textos 

Antes de empezar a escribir hay que elegir un buen editor de textos. Normalmente la universidad posee licencia para usar Microsoft© Word, si no, en internet podemos descargar OpenOffice Writer que es gratuito y muy similar.

El siguiente paso es aprender a usarlo y aprovechar todas las funcionalidades que ofrece, para ello es muy recomendable perder tiempo en leer un buen tutorial puesto que más tarde ganaremos con creces ese tiempo.

Tutorial de Microsoft© Word

Tutorial de OpenOffice Writer

2. Organizar la estructura del documento

Para establecer las secciones que va a contener el documento es recomendable leer y revisar otros documentos similares y adoptar su estructura o adaptarla a nuestras necesidades.

En el caso de un trabajo académico una estructura posible sería la siguiente:

1. Índice
2. Introducción y antecedentes
3. Objetivos
4. Materiales, métodos y resultados
5. Conclusiones
6. Bibliografía
7. Anexos

Se pueden consultar ejemplos de trabajos académicos en diversos sitios de internet: Rincón del Vago, Monografías ...

 

Para una tesis doctoral la estructura puede ser un poco más compleja:

1. Abreviaturas
2. Índice
3. Antecedentes y objetivos
4. Introducción
5. Materiales
6. Métodos
7. Resultados
8. Discusión
9. Conclusiones
10. Bibliografía
11. Anexos

Existen páginas web que recogen tesis doctorales de diferentes universidades, por ejemplo: Tesis Doctorales en Red, Teseo, Biblioteca Virtual Miguel de Cervantes ...

3. Introducir el diseño del documento en el editor de textos

Una vez tengamos claro el diseño de nuestro documento y sus secciones, introduciremos todos estos datos en el editor de textos elegido.

• Definiremos el tamaño, orientación y márgenes del documento, incluyendo los márgenes necesarios para encuadernarlo posteriormente.
• Daremos formato al encabezado y pie de página, incluyendo el título de las secciones del documento y los números de página.
• Fijaremos los estilos-formatos de los párrafos, los títulos de sección, las subsecciones, tablas, figuras y otros elementos del documento, así luego todo el documento tendrá una presentación uniforme. Esto supone definir el tipo de letra, tamaño, alineación, tabulaciones, interlineado... para cada uno de los elementos nombrados.
• Definir el estilo y numeración de títulos de tablas y figuras que se van a insertar en el documento. Se puede especificar qué tipo de título se utilizará automáticamente al insertar una imagen o una tabla y que éstos se autonumeren para no tener que modificar su numeración si se inserta una nueva imagen o tabla entre dos ya existentes.
• Es recomendable escribir los títulos de las secciones e insertar algún contenido (texto, tablas e imágenes) y aplicarles el formato adecuado para ir viendo como quedará el documento.
• Una vez insertados los títulos de las secciones y de algunas subsecciones de prueba se puede crear el Índice en la sección correspondiente del documente y fijar el estilo del mismo. Este índice se podrá actualizar automáticamente cada vez que se realicen modificaciones del documento.

4. Elegir y usar un gestor de bibliografía

El gestor de bibliografía es un programa de ordenador que nos ayudará a buscar, introducir, ordenar y guardar nuestras referencias bibliográficas. Además los gestores modernos incorporan la capacidad de insertar y dar formato automáticamente a la bibliografía para introducirla en nuestros documentos.

Actualmente existen varios gestores de pago: EndNote, ProCite, Reference Manager ...

Aunque también existen alternativas gratuitas: Scholar's Aid, Bibus ...

También existe Zotero que es una extensión del navegador Firefox para administrar bibliografía desde el mismo o CiteULike que es un gestor online de bibliografía, ambos son gratuitos.

Con estos programas podremos introducir de una manera rápida y sencillas las citas bibliográficas deseadas en nuestro documento de Microsoft© Word u OpenOffice Word, para más información de cómo hacerlo consultar sus manuales de uso.

 

5. Comenzar a escribir

Una vez tenemos todo preparado para comenzar a escribir. Es recomendable comenzar a escribir las secciones de las que más información tengamos recopilada y ordenada, cómo pueden ser 'Antecedentes y objetivos', 'Métodos' y 'Resultados', a su vez las secciones de 'Materiales', 'Bibliografía', 'Abreviaturas' y 'Anexos' se pueden ir completando según aparezcan sus contenidos al escribir los 'Antecedentes y objetivos', 'Métodos' y 'Resultados'.

Una vez completados los apartados anteriores y con todas las ideas expuestas y organizadas se escribirán los apartados de 'Discusión de los resultados' y 'Conclusiones'.

Es recomendable dejar para el final la 'Introducción', puesto que esta sección es de longitud muy variable y así la podemos adaptar al tiempo que nos quede tras completar el resto.

Durante el proceso es recomendable ir actualizando automáticamente el Índice para poder ver en pocas líneas el progreso de la escritura.

Por último queda dar el formato deseado a la Bibliografía, lo cual lo realiza automáticamente el programa gestor de la misma.

 

6. Diseñar una portada

La portada de nuestro trabajo será la que dé una primera impresión del mismo, por ello ha de ser acorde o incluso mejor que los contenidos.

Es aconsejable no usar únicamente texto para la misma, introducir imágenes que describan el contenido del trabajo en una buena composición gráfica ayudará a dar buena presencia a los contenidos. Para realizar las composiciones gráficas se pueden usar editores de imágenes como Photoshop, Paint Shop Pro o Gimp (éste último es gratuito) o incluso con Microsoft© Powerpoint.

Se deberán incluir en la portada: el título del trabajo, la institución a la que se pertenece (universidad, facultad, laboratorio), el nombre y apellidos del autor y del director o tutor del trabajo, la fecha y el motivo del trabajo. Los logos de la institución también se pueden incluir.

 

7. Otros 

En las tesis doctorales o trabajos que han significado gran exfuerzo y trabajo durante un largo período de tiempo se suelen incluir agradecimientos personales o un 'Prólogo' antes del contenido científico del trabajo.

También se incluye una página firmada por los directores del trabajo certificando la autenticidad y valía del mismo y se enuncian las becas o proyectos que han financiado el mismo.

 

8. Comprobación final

Tras concluir el documento, es aconsejable, si no necesario, releer varias veces el mismo entre varias personas para buscar erratas, errores ortográficos y gramaticales, u otros fallos que se pudieran haber cometido (sobretodo si se ha copiado texto de otras fuentes).

Revisar que los cambios de página no alteran las diversas secciones del documento, es aconsejable que las secciones comienzen en páginas impares para que aparezcan en hoja nueva a la derecha. Comprobar que los gráficos y tablas están correctamente posicionados y que todo tiene los formatos y estilos deseados.

Si el documento ha de ser revisado por terceras personas, esta comprobación final se puede dejar para después de las correciones.

 

9. Imprimir 

Una vez comprobado que está todo en orden llega la tan ansiada hora de imprimir el documento. Es aconsejable imprimirlo a doble cara (por motivos ecológicos) y en una impresora láser (una resolución de 300x300 puede ser más que adecuada) si el documento es en blanco y negro. Si el documento consta de páginas a color, estas pueden ser extraídas e impresas a parte en una impresora láser especial o en una de inyección de tinta. Tener en cuenta el establecer los márgenes apropiados para poder encuadernar porteriormente el documento.

Regresión lineal e intervalo de confianza

Una cuestión habitual en bioinformática es la de buscar e interpretar posibles correlaciones entre variables de un conjunto de datos. Cuando la correlación entre dos variables es significativa a menudo es interesante modelar la relación entre ellas por medio de una función lineal de regresión. Como todo modelo, este tipo de funciones tienen asociado un grado de error o incertidumbre, que podrá o no ser tolerable según el uso que queramos darle.

En este ejemplo mostramos gráficamente el error de una función de regresión, partiendo de un conjunto de datos, que guardaremos en un fichero llamado datos_regresion.tab, con datos separados con tabuladores:

days y2000 y2010
1 33.61 20.58
5 31.96 23.70
9 40.42 23.30
15 35.43 33.12
18 46.13 31.13
22 39.19 36.30
25 42.72 33.97
29 46.72 31.26
31 42.62 32.33
36 45.09 34.10 
39 55.05 38.17
43 46.89 37.46
46 54.98 40.08
51 51.83 40.68
54 54.64 40.68
57 49.87 39.47
60 51.33 39.57

De las tres variables vamos a trabajar con las dos primeras, days, que cuenta los días desde el inicio de lexperimento, y y2000, que se corresponde con las mediciones de la variableX en el año 2000. El siguiente código en R muestra como calcular la correlación, el coeficiente de determinación (y su P-valor) y el intervalo de confianza al 95% de la recta de regresión:

1:  # prog4.R , inspirado por ejemplos de   
2:  # http://wwwmaths.anu.edu.au/~johnm/r-book/2edn/scripts/exploit.R  
3:    
4:  # función: grafica el intervalo de confianza 0.95 en torno a la regresión  
5:  curvasIC <- function(form=varX$y2010~varX$days,data=varX,lty=2,  
6:     col=3,lwd=2,newdata=data.frame(days=varX$days))  
7:  {   
8:    varX.lm <- lm(form,data=data)   
9:    predmod <- predict(varX.lm,newdata=newdata,interval="confidence")   
10:    lines(spline(newdata$days,predmod[,"fit"]),lwd=3)   
11:    lines(spline(newdata$days,predmod[,"lwr"]),lty=lty,lwd=lwd,col=col)   
12:    lines(spline(newdata$days,predmod[,"upr"]),lty=lty,lwd=lwd,col=col)   
13:  }  
14:    
15:  varX = read.table("datos_regresion.tab",header=T)  
16:    
17:  ## nos centramos en la serie 2000   
18:  lm2000 = lm(varX$y2000~varX$days) # calcula regresión lineal  
19:  sum = summary(lm2000)  
20:    
21:  ## imprime diagrama XY   
22:  title = sprintf("year 2000 (R^2=%.2f)",sum$r.squared)  
23:  plot(varX$days,varX$y2000,ylab="variableX [arbitrary units]",  
24:     xlab="(days from experiment start)",main=title )  
25:    
26:  ## imprime regresión, intervalos de confianza y función ajustada  
27:  curvasIC(form=varX$y2000~varX$days)  
28:  text(40,35,sprintf("variableX (days) = %1.2f + %1.2f days",  
29:     lm2000$coefficients[1],lm2000$coefficients[2]))  
30:  text(40,32,sprintf("P-value = %1.3g",sum$coefficients[2,4]))  
31:    
32:  ## crea copia en formato PDF  
33:  dev.copy(pdf, file="regresion2000.pdf")  
34:  dev.off()  

Evaluación de un alineamiento de proteínas

Un alineamiento de secuencias proteicas se puede evaluar numéricamente comparando las identidades entre los aminoácidos alineados. Se puede asignar una puntuación a cada par de aminoácidos usando una matriz de sustitución, que describe el ritmo al que un aminoácido en una secuencia cambia a otro a lo largo de la evolución natural.

Las matrices de sustitución más utilizadas son las PAM y las BLOSUM. Éstas matrices son las que utiliza el famoso programa de alineamiento BLAST.

En el siguiente ejemplo se mostrará una sencilla  subrutina en Perl que evalúa un alineamiento de proteínas utilizando la matriz BLOSUM62. Las dos secuencias alineadas se pasan como sendas cadenas de caracteres:

1:  # prog3.pl Evalúa un alineamiento
2:  sub score_alignment {  
3:       my ($type,$score_matrix,$seq1,$seq2)=@_;  
4:       my (@indices,@matrix,@score_results);  
5:       if ($score_matrix eq 'blosum62'){  
6:            @indices = ('C', 'S', 'T', 'P', 'A', 'G', 'N', 'D', 'E', 'Q', 'H', 'R', 'K', 'M', 'I', 'L', 'V', 'F', 'Y', 'W');  
7:            @matrix = (  
8:                 [9, -1, -1, -3, 0, -3, -3, -3, -4, -3, -3, -3, -3, -1, -1, -1, -1, -2, -2, -2],  
9:                 [-1, 4, 1, -1, 1, 0, 1, 0, 0, 0, -1, -1, 0, -1, -2, -2, -2, -2, -2, -3],  
10:                 [-1, 1, 4, 1, -1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -2, -2, -2, -2, -2, -3],  
11:                 [-3, -1, 1, 7, -1, -2, -1, -1, -1, -1, -2, -2, -1, -2, -3, -3, -2, -4, -3, -4],  
12:                 [0, 1, -1, -1, 4, 0, -1, -2, -1, -1, -2, -1, -1, -1, -1, -1, -2, -2, -2, -3],  
13:                 [-3, 0, 1, -2, 0, 6, -2, -1, -2, -2, -2, -2, -2, -3, -4, -4, 0, -3, -3, -2],  
14:                 [-3, 1, 0, -2, -2, 0, 6, 1, 0, 0, -1, 0, 0, -2, -3, -3, -3, -3, -2, -4],  
15:                 [-3, 0, 1, -1, -2, -1, 1, 6, 2, 0, -1, -2, -1, -3, -3, -4, -3, -3, -3, -4],  
16:                 [-4, 0, 0, -1, -1, -2, 0, 2, 5, 2, 0, 0, 1, -2, -3, -3, -3, -3, -2, -3],  
17:                 [-3, 0, 0, -1, -1, -2, 0, 0, 2, 5, 0, 1, 1, 0, -3, -2, -2, -3, -1, -2],  
18:                 [-3, -1, 0, -2, -2, -2, 1, 1, 0, 0, 8, 0, -1, -2, -3, -3, -2, -1, 2, -2],  
19:                 [-3, -1, -1, -2, -1, -2, 0, -2, 0, 1, 0, 5, 2, -1, -3, -2, -3, -3, -2, -3],  
20:                 [-3, 0, 0, -1, -1, -2, 0, -1, 1, 1, -1, 2, 5, -1, -3, -2, -3, -3, -2, -3],  
21:                 [-1, -1, -1, -2, -1, -3, -2, -3, -2, 0, -2, -1, -1, 5, 1, 2, -2, 0, -1, -1],  
22:                 [-1, -2, -2, -3, -1, -4, -3, -3, -3, -3, -3, -3, -3, 1, 4, 2, 1, 0, -1, -3],  
23:                 [-1, -2, -2, -3, -1, -4, -3, -4, -3, -2, -3, -2, -2, 2, 2, 4, 3, 0, -1, -2],  
24:                 [-1, -2, -2, -2, 0, -3, -3, -3, -2, -2, -3, -3, -2, 1, 3, 1, 4, -1, -1, -3],  
25:                 [-2, -2, -2, -4, -2, -3, -3, -3, -3, -3, -1, -3, -3, 0, 0, 0, -1, 6, 3, 1],  
26:                 [-2, -2, -2, -3, -2, -3, -2, -3, -2, -1, 2, -2, -2, -1, -1, -1, -1, 3, 7, 2],  
27:                 [-2, -3, -3, -4, -3, -2, -4, -4, -3, -2, -2, -3, -3, -1, -3, -2, -3, 1, 2, 11]  
28:            );  
29:       }  
30:       if (length($seq1)==length($seq2)){  
31:            my @seq1 = split(//,$seq1);  
32:            my @seq2 = split(//,$seq2);  
33:            for (my $i=0; $i<=$#seq1; $i++) {  
34:                 # Find letters position in matrix  
35:                 my ($pos1) = grep { $indices[$_] eq $seq1[$i] } 0 .. $#indices;  
36:                 my ($pos2) = grep { $indices[$_] eq $seq2[$i] } 0 .. $#indices;  
37:                 if (defined($pos1) && defined($pos2)){  
38:                      push (@score_results, $matrix[$pos1][$pos2]);  
39:                 } else {  
40:                      push (@score_results, '-');  
41:                 }  
42:            }  
43:       } else {  
44:            die "# Both sequences must have the same length to score the alignment\n";  
45:       }  
46:       return @score_results;  
47:  }  

Zotero - Gestor de bibliografía web 2.0

Zotero es un gestor de bibliografía gratuito que funciona como extensión de Firefox. Su funcionamiento es similar al de otros gestores de pago como EndNote o ProCite y además tiene todas las ventajas de poder ser consultado online como CiteULike
Permite su uso tanto online como offline, además permite importar y exportar referencias en múltiples formatos (RIS, TeX...). Y lo más importante, sus plugins permiten insertar referencias bibliográficas tanto en Microsoft Word como en OppenOffice con gran facilidad y con cientos de estilos diferentes para elegir.
Otra gran característica de este gestor bibliográfico es su posibilidad de leer citas bibliográficas en páginas web como Amazon, Google Scholar o Pubmed e incorporarlas a nuestra propia librería.
En definitiva, Zotero es una de las opciones más interesantes en gestores bibliográficos en la actualidad.
Más información:

• http://www.zotero.org
• http://es.wikipedia.org/wiki/Zotero
• http://www.blogoff.es

Gráfico sencillo de dominios sobre secuencia (glifo)

El problema que nos ocupa hoy es el de representar gráficamente dominios o regiones sobre una secuencia de referencia. Quizás el ejemplo más típico sea el de mapear los resultados de una búsqueda BLAST sobre la secuencia de partida, que no siempre alinean contra la secuencia completa.
Por ejemplo, dada una secuencia de 935 aminoácidos, BLAST me devuelve un alineamiento que solapa desde la posición 23 a la 512. Cómo podemos representar esto gŕaficamente de manera sencilla? Es posible usar BioPerl para este fin, como se explica en el taller de (bio)perl, pero si queremos un diagrama como este,


que representa la longitud total de la secuencia y marca el segmento alineado a escala, es mejor recurrir al módulo GD, como se muestra en este programa:

1:  #!/usr/bin/perl -w
2:  # prog2.pl: crea barra PNG con uno o más dominios sobre secuencia de referencia
3:
4:  use strict;
5:  use GD;
6:
7:  # variables globales que definen el tamaño de la figura rectangular
8:  my $ANCHOFIGURA = 100;
9:  my $ALTOFIGURA = 8;
10:  my @COLORDOMINIO= (150,150,150); # en 3 canales RGB
11:
12:  # ejemplo: secuencia de 935 aminoácidos, con un dominio reconocido del 23 al 512
13:  if(pinta_dominios_secuencia(935,'23:512','panel.png'))
14:  {
15:    print "# archivo de salida: panel.png\n";
16:  }
17:
18:  sub pinta_dominios_secuencia
19:  {
20:    # argumentos:
21:    # 1) tamaño de secuencia de referencia
22:    # 2) cadena con dominios separados por compas, como: '23:78,123:189'
23:    # 3) nombre de archivo PNG de salida
24:    my ($length,$domain_definitions,$outPNGfilename) = @_;
25:
26:    my ($n_of_domains,$dom,@domains,$startd,$endd) = (0);
27:
28:    # comprueba los dominios
29:    foreach $dom (split(/,/,$domain_definitions))
30:    {
31:      if($dom =~ /(\d+):(\d+)/)
32:      {
33:        ($startd,$endd) = ($1,$2);
34:        if($startd<1 || $endd > $length)
35:        {
36:          print "# pinta_dominios_secuencia: dominio incorrecto ($dom)\n";
37:          next;
38:        }
39:
40:        push(@domains,[$startd,$endd]);
41:        $n_of_domains++;
42:      }
43:      else{ print "# pinta_dominios_secuencia: dominio incorrecto ($dom)\n"; }
44:    }
45:
46:    if(@domains)
47:    {
48:      my $imagen = GD::Image->new($ANCHOFIGURA,$ALTOFIGURA);
49:
50:      # define colores básicos en 3 canales RGB
51:      my $blanco= $imagen->colorAllocate(255,255,255);
52:      my $negro = $imagen->colorAllocate(0,0,0);
53:      my $color = $imagen->colorAllocate(@COLORDOMINIO);
54:
55:      # crea panel rectangular con dimensiones globales $ANCHOFIGURA y $ALTOFIGURA
56:      $imagen->filledRectangle(0, 0, $ANCHOFIGURA, $ALTOFIGURA, $blanco); # fondo
57:      $imagen->rectangle(0, 0, $ANCHOFIGURA-1, $ALTOFIGURA-1, $negro); # borde
58:
59:      # añade dominios
60:      foreach $dom (@domains)
61:      {
62:          # calcula dimensiones del dominio a escala, eje X
63:          $startd = int( $ANCHOFIGURA * ( $dom->[0] / $length ) );
64:          if($startd < 1){ $startd = 1 } # evita pisar el borde
65:          $endd  = int( $ANCHOFIGURA * ( $dom->[1] / $length ) );
66:          if($endd > $ANCHOFIGURA-2){ $endd = $ANCHOFIGURA-2 } # evita pisar el borde
67:          # con 1 y -2 evitamos el borde en eje Y
68:          $imagen->filledRectangle($startd, 1, $endd, $ALTOFIGURA-2, $color);
69:      }
70:
71:      open(PNG,">$outPNGfilename")
               || die "# pinta_dominios_secuencia: no puedo crear $outPNGfilename\n";
72:      binmode PNG;
73:      print PNG $imagen->png();
74:      close(PNG);
75:    }
76:
77:    return $n_of_domains;
78:  }

Extraer secuencias intergénicas de archivos GenBank

La tarea de hoy consiste en leer un archivo GenBank, por ejemplo de un genoma circular bacteriano, con el fin de extraer las secuencias intergénicas. Para esta tarea vamos a utilizar código de BioPerl y antes necesitaremos descargar al menos un fichero GenBank, como se explica por ejemplo en el taller de (bio)perl. Una vez tengamos los archivos de entrada necesarios, podemos resolver el problema de las secuencias intergénicas con una subrutina como ésta:

1:  use Bio::SeqIO;
2:
3:  # prog1.pl
4:  # recibe hasta 4 argumentos:
5:  # 1) archivo entrada en formato GenBank
6:  # 2) nombre del archivo de salida en formato FASTA
7:  # 3) longitud mínima de las regiones intergénicas (opcional, por defecto toma todas)
8:  # 4) longitud máxima de las regiones intergénicas (opcional)
9:
10:  sub extract_intergenic_from_genbank
11:  {
12:    my ($infile,$out_intergenic_file,$min_intergenic_size,$max_intergenic_size) = @_;
13:
14:    my ($n_of_intergenic,$gi,$start,$end,$length,$strand,$taxon) = (0);
15:    my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank' );
16:
17:    open(FNA,">$out_intergenic_file") ||
18:      die "# extract_intergenic_from_genbank : cannot create $out_intergenic_file\n";
19:
20:    while( my $seq = $in->next_seq) # scan all sequences found in $input
21:    {
22:      my ($gbaccession,$sequence,$gen,@genes) = ( $seq->accession() );
23:      $sequence = $seq->primary_seq()->seq() ||
24:         'no encuentro la secuencia primaria, necesito un fichero GenBank completo!';
25:      $taxon = '';
26:
27:      for my $f ($seq->get_SeqFeatures)
28:      {
29:        if($f->primary_tag =~ /CDS|rRNA|tRNA/) # campos de 'genes'
30:        {
31:          $gi = '';
32:          if($f->has_tag('db_xref'))
33:          {
34:            my $crossrefs = join(',',sort $f->each_tag_value('db_xref'));
35:            if($crossrefs =~ /(GI\:\d+)/){ $gi = $1 }
36:          }
37:          elsif($f->has_tag('locus_tag'))
38:          {
39:            if($gi eq '' && $f->has_tag('locus_tag'))
40:            {
41:              $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag'));
42:            }
43:          }
44:
45:          next if($gi eq '');
46:
47:          $start = $f->start();
48:          $end  = $f->end();
49:          $strand = $f->location()->strand();
50:          push(@genes,[$start,$end,$gi,$strand]);
51:        }
52:        elsif($f->primary_tag() =~ /source/)
53:        {
54:          if($f->has_tag('organism'))
55:          {
56:            foreach my $element ($f->each_tag_value('organism'))
57:            {
58:              $taxon .= "[$element],";
59:            }
60:            chop $taxon;
61:          }
62:        }
63:      }
64:
65:      # calcular ORFs vecinos y corta las secuencias intergenicas
66:      for($gen=1;$gen<scalar(@genes);$gen++)
67:      {
68:        next if($genes[$gen-1][1] >= $genes[$gen][0]); #no solapa asume genes en orden
69:        next if( defined($min_intergenic_size) && $min_intergenic_size > 0 &&
70:          $genes[$gen][0]-$genes[$gen-1][1]+1 < $min_intergenic_size); # evita cortas
71:        next if( defined($max_intergenic_size) && $max_intergenic_size > 0 &&
72:          $genes[$gen][0]-$genes[$gen-1][1]+1 > $max_intergenic_size); # evita largas
73:
74:        $n_of_intergenic++;
75:
76:        # calcula bordes de intergenes
77:        $start = $genes[$gen-1][1] + 1;
78:        $end  = $genes[$gen][0] - 1;
79:        $length = $end - $start + 1;
80:
81:        print FNA ">intergenic$n_of_intergenic|$gbaccession|coords:$start..$end|".
82:          "length:$length|neighbours:$genes[$gen-1][2]($genes[$gen-1][3]),".
83:          "$genes[$gen][2]($genes[$gen][3])|$taxon\n";
84:        print FNA substr($sequence,$start-1,$length)."\n";
85:      }
86:    }
87:
88:    close(FNA);
89:
90:    return $n_of_intergenic;
91:  }