domingo, 24 de mayo de 2015

Generación de árboles filogenéticos UPGMA a partir de un código de proteína, en Python


En este post se presenta una aplicación para la generación de árboles filogenéticos UPGMA a partir de un código de proteína en formato FASTA.
Los árboles filogenéticos son útiles para encontrar los familiares más cercanos a nuestro organismo de estudio. Se ha implementado mediante la creación de un script Python, que utiliza la biblioteca externa BioPython para la creación de árboles filogenéticos. La aplicación permite obtener secuencias similares a partir de una secuencia de proteína dada. Estas secuencias obtenidas por Blast son procesadas con un algoritmo de alineamiento múltiple (Clustal), cuyo resultado sirve para crear la matriz de distancias y el árbol filogenético UPGMA. La aplicación es multiplataforma, para sistemas Microsoft Windows y Linux.

La configuración permite seleccionar cuántas secuencias como máximo se quiere utilizar y el umbral de similitud para seleccionar las secuencias. El script se puede lanzar en modo atendido (menú de aplicación) y desatendido (para usar en procesos batch de procesamiento secuencial). Así tendremos las siguientes opciones de inicio de aplicación:

python filo_tree.py

Muestra ayuda de funcionamiento de la aplicación.

python filo_tree.py -g

Lanza el siguiente menú de aplicación:


Las opciones de aplicación son las siguientes:
  • CARGAR PROTEÍNA: Carga datos de proteína en formato FASTA a partir de un fichero.
  • VISUALIZAR PROTEÍNA: Visualiza la proteína cargada en formato FASTA.
  • BUSCAR SECUENCIAS SIMILARES: A partir de la proteína cargada se realiza una búsqueda de secuencias similares mediante el algoritmo blastp en una base de datos online. En esta opción se escoge tanto el número de secuencias a devolver como la base de datos a utilizar. Se genera un fichero XML con dichas secuencias.
  • VISUALIZAR SECUENCIAS SIMILARES: Se visualizan por pantalla las secuencias generadas por blastp que están guardadas en el fichero XML.
  • GENERAR ÁRBOL FILOGENÉTICO UPGMA: A partir de las secuencias guardadas en el fichero XML se generan las secuencias de alineamiento múltiple, con el fin de crear la matriz de distancias y poder crear el árbol filogenético UPGMA. En esta opción se escoge el umbral de similitud.
  • ACERCA DE: Créditos de la aplicación.
  • SALIR: Salir de la aplicación.

python filo_tree.py -p FICHERO.FASTA -nseq N_SEQ -e UMBRAL_SIMILITUD -bd N_BD

Lanza el script filo_tree.py, buscando secuencias similares a la secuencia proteica definida en formato FASTA en el fichero FICHERO.FASTA. El número de secuencias a buscar está definido por N_SEQ y el umbral de similitud por UMBRAL_SIMILITUD. Además se tiene que especificar la base de datos a la que conectarse para realizar la búsqueda, donde N_BD es un número entero correspondiente con:
  1. nr: Non-redundant GenBank CDS translations + RefSeq + PDB + SwissProt + PIR + PRF, excluding those in PAT, TSA, and evr_nr.
  2. refseq_protein: Protein sequences from NCBI Reference Sequence project.
  3. swissprot: Las major release of the UniProtKB/SWISS-PROT protein sequence database (no incremental updates).
  4. pat: Proteins from the Patent division of GenBank.
  5. pdb: Protein sequences from the 3-dimensional structure records from the Protein Data Bank.
  6. env_nr: Protein sequences translated from the CDS annotation of metagenomic nicleotide sequences.
  7. tsa_nr: Protein sequences translated from CDSs annotated on trascriptome shotgun assemblies

El script realiza alineamiento múltiple de las secuencias obtenidas mediante Clustal, generando la matriz de distancias y árbol filogenético UPGMA.

Ejemplo: Búsqueda de proteína contenida en H0XFT0.fasta en la base de datos de swissprot, recuperando 50 secuencias con e-value 1.

python filo_tree.py -p H0XFT0.fasta -nseq 50 -e 1 -bd 3

Devuelve lo siguiente:


Y el siguiente fichero que contiene la matriz de distancias y el árbol filogenético:


Esquema general de la aplicación


El esquema general de la aplicación consiste en:

  1. Pasar como argumento (parámetro) un fichero con una secuencia de proteína en formato FASTA, al algoritmo BLAST (implementado en una función dentro del paquete BioPython). Dicha función busca un número de secuencias determinado y las alinea por pares a partir de la secuencia FASTA dada. Es decir, alinea N secuencias a partir de una (no es alineamiento múltiple, es alineamiento por pares). Esta información (la secuencia FASTA original, la secuencia similar y el alineamiento, por cada secuencia similar encontrada) se guarda en un fichero XML.
  2. A partir del fichero BLAST_FILE.XML se obtienen las secuencias similares (no los alineamientos),y se guardan en un fichero en formato FASTA (FASTA_MULTIPLE.FASTA). Por tanto ahora tenemos las N secuencias similares en un formato adecuado, listas para un alineamiento múltiple.
  3. Se realiza un alineamiento múltiple a partir del algoritmo de Clustal por medio de la aplicación externa clustalw2, que se puede obtener de http://www.clustal.org/. El script Python busca la aplicación en el PATH y lanza el programa (algoritmo), guardando la salida (el alineamiento múltiple) en el fichero FASTA_MULTIPLE.aln.
  4. El fichero de alineamiento múltiple es el parámetro de entrada para el cálculo de la matriz de distancias (implementado en la clase DistanceCalculator de BioPython). Los datos de la matriz de distancia se vuelcan o a un fichero de salida TREE_UPGMA_fecha_hora.TXT o por pantalla. 
  5. La variable DISTANCE MATRIX (matriz_distancia) es a su vez el parámetro de entrada del método upgma del objeto producto de la instanciación de la clase DistanceTreeConstructor para la construcción del árbol UPGMA:
  6. Una vez calculado el árbol filogenético UPGMA ya se pueden obtener todos los datos, volcando a fichero (TREE_UPGMA_fecha_hora.TXT) o pantalla los datos propios del árbol, la representación del árbol en ASCII (mediante la función Phylo.draw_ascii de BioPython) y la generación de un fichero gráfico PNG (TREE_UPGMA_fecha_hora.PNG) a partir de Phylo.draw_graphviz de BioPython. Opcionalmente se puede visualizar por pantalla mediante Phylo.Show() de BioPython.
Parámetro de entrada - Fichero FASTA 

La aplicación tiene que tener como parámetro de entrada un fichero con un código de proteína en formato FASTA. Por ejemplo el fichero H0XFT0.fasta (adjunto con el script a modo de ejemplo):
Este fichero se puede obtener desde el portal web de UNIPROT http://www.uniprot.org/.

Saludos.

sábado, 11 de abril de 2015

Script para calcular tamaño de ficheros pasados como parámetros

Hola. En este post presento un script muy sencillo para calcular el total del tamaño de los ficheros pasados como parámetros. Se hace uso del módulo os para tratamiento de ficheros.

# -*- coding: utf-8 -*-
#-------------------------------------------------------------------------------
# Name:        sizefiles.py
# Purpose:
#
# Author:      Ángel luis
#
# Created:     10-04-2015
# Copyright:   (c) Ángel 2015
# Licence:     GPL 3
#-------------------------------------------------------------------------------

import os
import sys

def main():
    # Get files.
    if len(sys.argv) == 0:
        raw_input('I need some file names. Press any key...')
    else:
        # Arguments.
        files = [i for i in sys.argv]
        try:
            total_size = 0
            for i in files[1:]: total_size += os.path.getsize(os.path.realpath(i))
            raw_input("Total size: %d bytes. Press any key..." % total_size)
        except:
            raw_input("Files exists?. Press any key...")

if __name__ == '__main__':
    main()

La forma de ejecutarlo es python sizefiles.py fichero1 fichero2 ficheroN. Devuelve la suma del tamaño en bytes de los ficheros pasados como parámetros.

Saludos.

miércoles, 8 de abril de 2015

Compilar Tinker Molecular Modeling en Linux

Hola. En este post vamos a ver como compilar el software de modelización molecular Tinker.


Lo primero es descargarse el código fuente, que se puede encontrar en:

http://dasher.wustl.edu/tinker/

Instalamos la aplicación en Linux normalmente.

Para compilar Tinker podemos utilizar el compilador Fortran G95, que se obtiene de:

http://www.g95.org/downloads.shtml

Una vez descargado Fortran, se instala normalmente (es un paquete DEB).

Una vez instalado Tinker y el compilador Fortran G95, los pasos a seguir son los siguientes:

1) Nos vamos a la carpeta MolecularTools, que es donde está instalado Tinker.



2) En MolecularTools/tinker/linux/gfortran tenemos los ficheros make para compilar la aplicación. A estos ficheros les damos permisos de ejecución. He hecho chmod 777 sobre todos ellos.


3) Estos ficheros (compile.make, library.make, link.make, rename.make) se copian al directorio MolecularTools/tinker/linux/source.



4)  Nos vamos a este directorio donde se encuentra el código fuente y vamos ejecutando uno a uno los fichero make:


./compile.make
./library.make
./link.make
./rename.make


5) En la carpeta MolecularTools/tinker/bin aparecen nuestros ficheros compilados:


Ya podemos ejecutar perfectamente Tinker, en MolecularTools, mediante Force Field Explorer.


¿Por qué compilar? Hay veces que se instala Tinker con unos binarios por defecto, y a la hora de generar la molécula simplemente no la visualiza. Es necesario compilar la aplicación, y una forma de hacerlo es mediante el compilador Fortran G95.


Saludos.