Genómica evolutiva 2021


En bioinformática es muy común tener los datos guardados en un archivo simple de texto plano. Estos archivos toman distintos formatos según el programa que lo genere, por ejemplo fasta, FASTQ, GFF, GTF, SAM, bed, y un largo etc, pero todos ellos son archivos simples de texto plano. Transformar archivos de texto de un formato a otro suele ser una tarea de todos los días en bioinfo.

En este práctico nos vamos a concentrar en el uso de los comandos grep, sed, y awk, y en las posibilidades que nos brindan para el procesamiento de archivos de texto en los típicos formatos de salida que manejan los programas bioinformáticos.


Antes de empezar con el contenido del práctico, descarguen a sus máquinas virtuales los archivos que vamos a utilizar durante el práctico, ubicados en el servidor rnaseq, en /home/estudiantes/archivos_p2.

usuario@equipo:~$ scp -r estudiantes@rnaseq.fcien.edu.uy:/home/estudiantes/archivos_p2 .

Para aquellos que estén trabajando en rnaseq, copien los archivos a su directorio personal con el siguiente comando:

usuario@equipo:~$ cp -r /home/estudiantes/archivos_p2 .


grep (Global Regular Expression Print)

grep busca un patrón o una expresión en los archivos de entrada e imprime a la salida estándar las líneas que lo contengan. Comenzando en la primera línea del archivo, grep la compara con la expresión de búsqueda, y si la comparación da resultado positivo imprime la línea. El proceso se repite línea a línea hasta que llega al final del archivo.

El formato del comando es:

usuario@equipo:~$ grep patrón file

o

usuario@equipo:~$ {comando} | grep patrón

Esto devuelve todas las líneas que contienen una cadena de texto que coincide con la expresión patrón en la entrada, sea un archivo o la salida de un comando anterior. Por ahora pensaremos en las expresiones como cadenas de texto, entonces grep devuelve todas las líneas que contienen patrón como una subcadena.

Ejemplo:

usuario@equipo:~$ cat ex1.fa #inspeccionamos el contenido del archivo ex1.fa

usuario@equipo:~$ grep ">" ex1.fa #buscamos las líneas que contienen el caracter ">"
>seq1
>seq2


Opciones útiles del comando grep


Ejemplos:

# buscamos las líneas que NO contienen el caracter ">"
usuario@equipo:~$ grep -v ">" ex1.fa

# contamos la cantidad de veces que aparece el carácter ">
usuario@equipo:~$ grep -c ">" ex1.fa


Expresiones regulares en grep


Las expresiones regulares son cadenas de caracteres basadas en reglas sintácticas que permiten definir de manera abstracta patrones de búsqueda. Cada carácter en una expresión regular es o bien un meta-carácter, que tiene un significado especial, o un carácter regular que tiene un significado literal. El meta-carácter por excelencia es el punto . , que coincide con cualquier carácter. Veamos un ejemplo:

usuario@equipo:~$ grep seq. ex1.fa
>seq1
>seq2

Notar que el “.” coincide tanto con “1” como con “2”.


Cuando trabajamos con expresiones regulares en grep, en general nos va a convenir usar la opción -E, que nos da un soporte extendido para expresiones regulares.


Operadores de repetición 


#buscamos las líneas que contienen las secuencias "AAATCC" o "AAACC", es decir que la "T" puede estar o no
usuario@equipo:~$ grep -E AAAT?CC ex1.fa



#similar a la anterior, pero ahora puede haber un número indeterminado de "T"s
usuario@equipo:~$ grep -E AAAT*CC ex1.fa


Listas de caracteres y rangos


Para hacer coincidir con un carácter de una lista de caracteres, usamos [ ]. Dentro de los corchetes van todos los caracteres admitidos. Por ejemplo:

#buscamos las líneas que contienen las secuencias "AAATCC" o "AAAGCC"
usuario@equipo:~$ grep -E AAA[TG]CC ex1.f

También se permiten los rangos de caracteres: [0-9] encuentra cualquier número y [a-zA-Z] encuentra cualquier letra del alfabeto.

Los paréntesis rectos también pueden ser usados para encontrar no-coincidencias, usando el símbolo ^ como el primer carácter dentro de los paréntesis. Por ejemplo, [^actgACTG] puede ser usado para encontrar las líneas que contengan caracteres que no sean uno de los cuatro nucleótidos.


Ancla al comienzo y final de línea


Cuando el carácter ^ se encuentra al principio de la expresión simboliza el comienzo de línea; cuando el carácter $ se encuentra al final de la expresión simboliza el final de línea. Estos operadores, conocidos como anclas, son útiles de muchas maneras. Supongamos que queremos encontrar una línea que no sólo contenga nuestra expresión sino que la misma se corresponda con la linea completa, o que se encuentre al comienzo o al final de la línea. Veamos un ejemplo con un archivo en formato FASTQ.

El formato FASTQ es un formato basado en texto para almacenar tanto una secuencia biológica (generalmente una secuencia de nucleótidos) como sus puntajes de calidad correspondientes. El formato usa 4 líneas por secuencia:

+ La línea 1 comienza con un `@` seguido del identificador de la secuencia.

+ La línea 2 tiene la secuencia en sí.

+ La línea 3 comienza con un `+` opcionalmente seguido del mismo identificador de la línea 1.

+ La línea 4 codifica los valores de calidad para la secuencia en la línea 2 y debe contener el mismo número de símbolos que caracteres en la secuencia.

Para contar cuantas secuencias hay en un archivo fastq, busquemos las líneas correspondientes al identificador de la secuencia.

usuario@equipo:~$ head ex2.fastq #examinamos el archivo fastq
usuario@equipo:~$ grep -c ^@SRR ex2.fastq #contamos los encabezados
50


Caracteres escapados


Son caracteres precedidos de una barra invertida. Dicha barra invertida realiza una de las dos siguientes acciones: (a) elimina un significado especial implícito de un carácter, o (b) agrega un significado especial a un carácter “no especial”.

Por ejemplo, si queremos buscar una expresión que contiene espacios, podemos o bien encerrar nuestra expresión entera entre comillas, o “escapar” el espacio, anteponiendo al carácter " " una barra invertida. Ejemplo:


Ejercicios:


sed (Stream Editor)


sed es un editor de flujos de texto que lee y modifica un texto línea a línea de acuerdo a un script. Los comandos de edición de sed permiten entre otras cosas sustituir o borrar. En general, la forma en que usamos sed es la siguiente:

usuario@equipo:~$ sed 'comando_sed' file
usuario@equipo:~$ {comando} | sed 'comando_sed'

Por defecto, el resultado se escribe a la pantalla. Puede redirigir la salida a un nuevo archivo con el operador >, o si desea editar el archivo existente en su lugar, debe usar la opción -i. Tener cuidado con esta opción, el resultado es una modificación irreversible.

Las expresiones regulares de sed son esencialmente las mismas que las de grep. Los caracteres especiales también son esencialmente los mismos, pero con algunas diferencias clave como la barra / y el &, que son caracteres especiales en sed.


El comando sustituir


En su forma más simple, el formato del comando es sustituir es:

usuario@equipo:~$ sed 's/patrón/texto_nuevo/flag' archivo

No vamos a discutir los distintos flags, por ahora sólo vamos a utilizar el flag g (por global), que reemplaza todas las coincidencias que encuentra. Veamos un ejemplo sencillo.

Supongamos que queremos cambiar los encabezados de nuestro archivo fasta, cambiando la palabra “seq” por “yeast”. El comando de sed sería:

usuario@equipo:~$ sed 's/seq/yeast/g' ex1.fa > yeast.fa
usuario@equipo:~$ head yeast.fa #examinamos el archivo con los encabezados editados


El ampersand (&)


El ampersand es usado para representar el patrón de búsqueda en el comando sustituir. Cualquiera que sea el texto que coincida con el patrón definido, puede utilizar el ampersand para recuperarlo en el patrón de reemplazo. Ejemplo:

usuario@equipo:~$ sed 's/seq/&_yeast/g' ex1.fa > seq_yeast.fa
usuario@equipo:~$ head seq_yeast.fa #examinamos el archivo con los encabezados editados

En este ejemplo el “&” representa el patrón de búsqueda “seq”, con lo que "&_yeast" se expande a “seq_yeast”.


El comando delete

La sintaxis del comando delete de sed es muy sencilla: '/patrón/d'. Por ejemplo:


# eliminamos las líneas que contienen "seq", es decir los encabezados
usuario@equipo:~$ sed '/seq/d' ex1.fa


Ejercicio:


El formato tabular de blastn


Antes de introducir el siguiente comando vamos a ver un formato de salida de blastn con el que vamos a trabajar durante el resto del curso. blastn se puede correr localmente (en nuestras máquinas o en el servidor remoto) por línea de comando. La salida puede tener distintos formatos. Uno de los más útiles para trabajar con herramientas de línea de comando como awk es el formato tabular 6, donde cada columna corresponde a un dato del mapeo de blast. Por ejemplo, la columna 1 corresponde al query ID, la columna 2 al subject ID, la 3 al porcentaje de identidad, y la 4 al largo del alineamiento. Puede a qué corresponde el resto de las columnas en http://www.metagenomics.wiki/tools/blast/blastn-output-format-6.


awk


awk a la vez es un lenguaje de scripting para escaneo y procesamiento de patrones de texto, creado por Aho, Weinberger y Kernighan, y es también el programa que interpreta ese lenguaje, al que invocamos con el comando awk. El programa opera línea por línea en el archivo de entrada y permite el uso de variables, funciones numéricas, operaciones lógicas, edición de texto, etc. Puede ser bastante sofisticado, por lo que esta guía es apenas un comienzo, pero debería darle una idea de lo que awk puede hacer.

awk permite escribir scripts en forma de declaraciones que definen condiciones que se deben evaluar en cada línea, y acciones que deben ser ejecutadas cuando la condición evalúa positivamente. Si no hay condiciones, las acciones se ejecutan en todas las líneas. Las condiciones pueden incluir operaciones lógicas, expresiones relacionales o ‘coincidencia de patrones’ (incluyendo expresiones regulares).


Flujo de trabajo de un script awk

Un script awk tiene una sección opcional BEGIN { }, en la cual las acciones se ejecutan antes de procesar cualquier contenido del archivo; luego la sección principal con las condiciones a evaluar y acciones a ejecutar en cada línea del archivo, y finalmente hay una sección END { }, que también es opcional y sus acciones suceden después de que la lectura del archivo haya terminado:

BEGIN {… comandos opcionales de inicialización awk …}
condición {… acción …}
condición {… acción …}
END { … finalización comandos awk …}


El flujo de trabajo de un script de awk es el siguiente:

De forma general, la forma en que ingresamos un script de awk por línea de comando es:

usuario@equipo:~$ awk 'condiciones {acciones}' file

o

usuario@equipo:~$ {…comandos…} | awk 'condiciones {acciones}'


Campos


awk considera que cada línea está compuesta por varios campos, cada uno separado por un “separador de campo” (que por defecto es uno o más espacios, o un tabulador). Tomemos como ejemplo la cadena “campo1 campo2 campo3”: awk ve esta cadena como compuesta por tres campos. Dentro de awk, el primer campo es referido como $1, el segundo $2, etc. Entonces, en nuestro ejemplo, $1 corresponde a la cadena “campo1”.

usuario@equipo:~$ echo "campo1 campo2 campo3" | awk '{print $1}'
campo1

El separador de campos puede ser cambiado mediante la opción -F. Por ejemplo, si seteamos -F"," el separador va a ser la coma.

usuario@equipo:~$ echo "campo1,campo2,campo3" | awk -F"," '{print $2}'
campo2


Variables internas predefinidas de awk



Operadores

Las condiciones de un script de awk pueden incluir operaciones lógicas, expresiones relacionales, operaciones aritméticas, o coincidencia de patrones. Para esto se usan una serie de operadores que se describen a continuación:

OPERADORES RELACIONALES

OPERADORES LÓGICOS (BOOLEANOS)

OPERADORES ARITMÉTICOS

BÚSQUEDA DE PATRONES


Ejemplos:

# imprimimos los hits con %ID > 99
usuario@equipo:~$ awk '$3 >= 99 {print $0}' blast.bl

# lo mismo pero agregando la condición de que el largo del alineamiento sea mayor a 10000
usuario@equipo:~$ awk '$3 >= 99 && $4 > 10000 {print $0}' blast.bl


Vamos a obtener el promedio del largo de los alineamientos para aquellos hits cuya identidad es mayor a 99% y el largo del alineamiento es mayor a 10000.

usuario@equipo:~$ awk '$3>99' blast.bl | awk 'BEGIN {sum=0} {sum=sum+$4} END {print sum/NR}'
12861.1

La variable sum que declaramos en la sección BEGIN va guardando la suma de los valores del campo 4, que corresponde al largo del alineamientos. Cuando se terminan de procesar las líneas, en la sección END imprimimos el valor de la suma dividido el número de líneas, que está dado por la variable interna NR. Esto nos da el promedio.

¿Por qué piensa que hicimos un pipeline con dos comandos awk en lugar de poner el $3>99 como condición dentro del mismo comando awk?


Ejercicio: