Contar elementos en un fichero con fgrep + wc

¡Ajá! Muchos pensábais que estaba acabado, pero no va a pasar un año sin que publique nada por aquí.

Os presento un comando que sirve para contar cuántas veces se repite una cadena de caracteres en un fichero.

$ fgrep -o ">" file.fas | wc -l

Este ejemplo sirve para contar el número de secuencias de un fichero FASTA. fgrep imprime la cadena cada vez que la encuentra en el fichero y wc cuentra el número de veces que lo repite.

Hasta la próxima.

Correr scripts en Perl como un ejecutable

Un simple pero muy útil y eficiente truco es crear una carpeta ~/bin donde añadir binarios de programas y así ejecutarlos rápidamente sin tener que ir a la localización de los mismos como expliqué aquí.

En el caso de que se trate de un script en perl tenemos que modificar el fichero. Si se encuentra en la carpeta bin, debemos de comprobar que la primera línea resulta ser al como esto.

#! /usr/bin/perl

Puede ocurrir que no tenga permisos de ejecución, para dárselos bastaría con hacer en consola:

$ chmod +x script.pl

Así de fácil. Así de útil.

 

 

Compilando r8s en Ubuntu 8.04

No te preocupes, no has retrocedido 3 años en el tiempo. Este tutorial sirve para compilar un software de datación de filogenias llamado r8s. Hay que decir que la última versión 1.71 fue lanzada como en 2006 y que los compiladores que requiere están obsoletos y no se encuentran para versiones recientes, al menos, para Debian/Ubuntu. Por este motivo, probé a compilar en una versión antigua de Ubuntu, la 8.04 LTS “Hardy Heron”. Me bajé la imagen ISO de la distribución y la instalé en VirtualBox. Una vez dentro de Ubuntu 8.04 me descargué el código de r8s v1.71. Con esto sobre el escritorio, comienza el manual:

1. Instalar el paquete “atool” y descomprimir el fichero:

$ sudo aptitude install atool

$ aunpack r8s1.7.dist.tar.Z

2. Instalamos compiladores de Fortran. No sé cuáles de ellos son estrictamente necesarios, en mi caso instalé los siguientes paquetes:

$ sudo aptitude install libg2c0-dev cfortran fortran77-compiler fortran95-compiler g77

3. En la carpeta que anteriormente descomprimimos, nos dirigimos a el subdirectorio “src” y compilamos:

$ cd r8s1.71/src

$ make

Si no da ningún mensaje de error, debería de estar compilado. Para comprobar que el proceso ha transcurrido correctamente, podemos ejecutar r8s:

$ ./r8s

Y debería de aparecernos un prompt consistente en “r8s>”.

Bueno, ahora queda entender cómo funciona. pero eso ya será en otro episodio.

[Bioinformática] DendroPy: Librería para análisis filogenético. Árbol consenso MrBayes

La informática tiene sus misterios, o tal vez yo soy lo suficientemente avispado como para escudriñar hechos que aparentemente no tienen lógica. Resulta que tras correr una larga simulación para realizar una inferencia filogenética bayesiana con MrBayes, el software era incapaz de realizar la filogenia consenso con la misma orden que comúnmente he ido utilizando desde que aprendí a correr el programa.

Buscando soluciones, me topé con DendroPy, una librería escrita en Python para análisis filogenético.

Para instalarlo, puedes utilizar uno de estos dos comandos.

$ sudo pip install dendropy

con utilizando setuptools:

$ sudo easy_install -U dendropy

Si todo ha marchado correctamente, podremos utilizar SumTrees, el programa de la librería apropiado para realizar una filogenia consenso con los resultados de MrBayes.

Por ejemplo, para tener un resumen de la probabilidad posterior en la filogenia, la línea sería esta:

$ sumtrees.py --min-clade-freq=0.95 --burnin=200 --support-as-labels --output=result.tre treefile1.tre treefile2.tre treefile3.tre

Gracias a esto logré aprovechar ese análisis. Siempre Python.

[Bioinformática] Orthoptera Species File: Base de datos taxonómica de ortópteros

Amantes de los grillos y saltamontes, si queréis tener información taxonómica sobre vuestro orden de insectos preferido, no dudéis en visitar Orthoptera Species File en la siguiente dirección:

http://orthoptera.speciesfile.org/HomePage.aspx

Vuestros amigos os lo agradecerán.

Un cordial saludo.

[Bioinformática] Texshade: Alineamientos de nucleótidos o proteínas en LaTeX (LyX)

Si necesitas añadir alineamientos de nucleótidos o proteínas en tu documento LaTeX, TeXShade es tu paquete. La manera de trabajar con este paquete es muy sencilla:

En el preámbulo añadimos el paquete con esta línea:

\usepackage{texshade}

En el cuerpo del documento añadimos los comandos, donde introducimos la dirección del fichero del alineamiento en formato FASTA, y seguidamente algunas características que queremos que aparezcan en la edición del alineamiento:

\begin{texshade}{/direccion/del/fichero/gb.fas}
\feature{ttop}{1}{1..528}{bar:conservation}{}
\showfeaturestylename{ttop}{conservación}
\end{texshade}

Esto es parte del resultado, en mi caso lo he combinado con el paquete landscape para que las páginas que son ocupadas por el alineamiento estén apaisadas.

Más información | Manual TeXShade

[Bioinformática] abifpy: Módulo de Biopython para leer ficheros ABI *.ab1

Llevaba tiempo buscando alguna manera de introducir secuencias en formato abi en Biopython. Este formato contiene las lecturas del secuenciador así que no solo tenemos la información de la secuencia sino de la calidad o el nivel de confianza de dada una de las bases secuenciadas.

El módulo de Biopython se llama Abifpy y está disponible para descargar libremente desde el Github del autor. Una manera de instalarlo es introducir la ruta del módulo en el fichero ~/.bashrc:

# Add custom Python modules to the Python path.
PYTHONPATH=$PYTHONPATH:~/abifpy
PYTHONPATH=$PYTHONPATH:~/usr/lib/python2.6/site-packages
export PYTHONPATH

Ahora al iniciar  Python deberíamos tener el módulo operativo. La forma de uso es muy sencilla:

$ python
>>> import abifpy
>>> secuencia = abifpy.Trace(“file.ab11″) # abrir fichero
>>> secuencia = abifpy.Trace(“file.ab1″, trimming=True) # abrir fichero y quitar los extremos de la secuencia, de menor calidad
>>> secuencia.seq # ver secuencia
>>> secuencia.qualVal # ver valores de calidad de cada secuencia
>>> secuencia.name # nombre de la secuencia

Fácil, ¿verdad? Pues tan más útil resulta.

[Bioinformática] Trabajando con cantidades genómicas de secuencias mediante Galaxy

Actualmente estoy aprendiendo a trabajar con secuencias del genoma humano. Como son cantidades ingentes de secuencias es necesario un software específico y, a ser posible, que te permita trabajar desde la web para que nuestro ordenador no inicie una combustión espontánea por sobrecarga . Para esta tarea estoy utilizando el software Galaxy, desarrollado por la universidad de Pennsylvania (EE UU) que integra multitud de herramientas para el análisis de genomas. Además tienen el software ejecutable desde un servidor al que se accede a través de esta web:

http://main.g2.bx.psu.edu/

En este minitutorial, explicaré cómo extraer secuencias Alu del genoma en formato FASTA. Primeramente nos registramos en el servidor para que recuerde nuestro trabajo (opcional), para ello nos dirigimos a la pestaña “User”. Una vez logueados comenzamos a trabajar.

1. Obtener datos del genoma. Nos dirigimos a Get Data >> UCSC Main table browser, y nos abrirá el buscador de la Universidad de California Santa Cruz. Aquí seleccionamos Genoma Humano, Ensamblaje de 2006 (actualmente está mejor anotado que el de 2009). Ahora seleccionamos el tipo de secuencia con la que queremos trabajar, en este caso Variation and Repeats y en track seleccionamos RepeatMasker como software rastreador de las secuencias. Después restringimos la búsqueda en Filter >> Edit >> RepName >> Alu*, y aceptamos haciendo clic en “Submit”. Seleccionamos como formato de salida BED. Iniciamos la búsqueda con “get output” y seguidamente “Send query to Galaxy”. Pasado un tiempo, en la barra lateral de la derecha de Galaxy se cargarán nuestros datos y podremos editar algunos atributos como el nombre algo importante para distinguirlo de las salidas en posteriores pasos. La salida es una tabla donde para cada secuencia nos indica el cromosoma donde se sitúa, el comienzo, el final, el nombre y la cadena en la que se dispone.
2. Calcular tamaño de la secuencia. Hemos obtenido más de 1.100.000 secuencias, cifra que es muy grande por lo que es interesante reducirla por ejemplo, filtrando según el tamaño. Sin embargo, esa es una información que no tenemos en los datos obtenidos, pero que podemos calcular a partir de las columnas de comienzo y fin de la secuencia. Para ello nos dirigimos a Text Manipulation >> Compute. En “Add Expression” indicamos que queremos restar el número de la columna 2 al número de la columna 3 con expresiones regulares, en este caso escribiremos “c3-c2″. Selecionamos el conjuto de datos en el que queremos que se aplique la operación y pulsamos “Submit”. Tras este paso tendremos un nuevo conjunto de datos con una columna adicional que postrará el resultado de la resta para cada secuencia.
3. Filtrar secuencias por tamaño. Vamos a ordenar las secuencias obtenidas según tamaño, para ello en Filter and Sort >> Sort data in ascending or descending order, seleccionamos el conjunto de datos e indicamos que queremos ordenar según “c7″. Así podemos hacer una estimación visual de a partir de qué tamaño nos interesa filtar, acción que se realiza en Filter and Sort >> Filter data on any column using simple expressions, seleccionamos conjunto de datos y añadimos la condición “c7 > 300″. Esto nos dará una tabla con las secuencias de un tamaño superiora 300bp.
4. Obtener secuencia en formato FASTA.Para obtener la secuencia propiamente dicha y exportarla a un fichero formato FASTA nos vamos a Fetch Sequences >> Extract Genomic DNA using coordinates from assembled/unassembled genomes, y tras seleccionar datos y formato FASTA, ejecutamos la aplicación y al tiempo tendremos otra serie de datos en la barra derecha. Para exportarla a nuestro equipo hacemos clic en el nombre del conjunto de datos y presionamos sobre el diskette para guardar estos datos. Finalmente tendremos en nuestro disco duro un fichero FASTA con las secuencias Alu de una tamaño superior a 300bp anotadas en el genoma humano ensamblado en 2006.

Esto es todo, espero que os haya interesado tanto como a mí. Hasta el próximo tutorial.

[Bioinformática] Ejecutar BLAST de manera local con base de datos propia

Una de las primeras herramientas bioinformáticas que aprendemos a usar es BLAST. Es un algoritmo que realiza comparaciones de una secuencia problema con una base de datos y nos devuelve como salida un alineamiento local de aquellas secuencias de la base de datos que presentan mayor similitud con la problema.

Muchos de vosotros habréis ejecutado BLAST desde la página de NCBI, enfrentando vuestra secuencia contra la base de datos de alguno de los organismo modelo más estudiados en genética o, si como yo estudiáis bichos raros, contra alguna genérica. También existen otras vías para ejecutar en remoto estas BLAST contra estas bases de datos fuera del navegador como, por ejemplo, instalando BLAST en local o con BioPython.

Pues bien, podemos ejecutar localmente BLAST de manera simple con una base de datos propia. Para ello, necesitamos instalar BLAST. En mi caso los he instalado desde los repositorios de Bio-Linux:

$ sudo aptitude install bio-linux-blast bio-linux-blast+

Seguidamente, creamos nuestra base de datos a partir de un fichero con secuencias en formato FASTA de la siguiente manera:

$ formatdb -i grasshoppers.fasta -p F -n ghoppers

Y una vez creada la base de datos, ejecutamos BLAST para que nos devuelva la salida en HTML:

$ blastall -p blastn -d ghoppers -i miseq.fasta -T -o output.html

Esta orden nos generará un fichero HTML con las puntuaciones y los alineamientos que podemos abrir con nuestro navegador preferido.

Otra línea más compleja puede ser:

$ blastall -p blastn -d ghoppers -i miseq.fasta -v 10000 -e 0.0001 -m 7 -o output.xml

Que resulta un fichero XML con hasta 10000 hits con e-valor de hasta 0.0001.

Y esto es todo, espero que haya sido de vuestro interés.

Más información | Bio-Linux

[Bioinformática] UNAFold: Software para predecir estructura secundaria de secuencias nucleotídicas

Para mi trabajo necesito predecir el plegamiento de una secuencia de ARN y, para ello iba a utilizar UNAFold, una versión renovada actualizada de mfold que es mucho más fácil de instalar, que permite analizar múltiples secuencias y una mayor variedad de cálculos. Aunque para instalar UNAFold no es necesario mfold, requiere de mfol_util, pues contiene boxplot_ng y sir_graph, para obtener gráficos de las estructuras y de puntos de UNAFold.

Para tener una instalación completa del software, necesitamos instalar mfold_util. Para ello, desde la página de descarga de mfold, descargamos mfold_util. Descomprimimos y nos dirigimos a la carpeta para compilar e instalar, entonces:

$ cd mfold_util-4.6

$ ./configure

$ make

$ sudo make install

Listo. Ahora instalamos UNAFold. Descargamos de la página de descargas de UNAFold un fichero rpm que nos facilitan los desarrolladores. Es un fichero de instalación nativo de RedHat, por lo que para instalarlo en Ubuntu necesitamos convertirlo a deb con alien. Desde consola todo esto sería:

$ sudo aptitude install alien

$ sudo alien –to-deb unafold-3.8-1.i386.rpm

$ sudo gdebi unafold-3.8-1.i386.deb

Y ya podemos ejecutar UNAFold. Por ejemplo un alinemiento en formato fasta:

$ UNAFold.pl alineamiento.fas

Ahora llega lo peor, ponerse a interpretar tanto resultado.