1.10: Filas y Columnas

Volvamos a la salida de las proteínas de levadura versus proteínas de levadura Autoexplosión que realizamos previamente, desde el archivo yeast_blastp_yeast_top2.txt.

Preguntemos lo siguiente: ¿cuántas secuencias en este conjunto de datos tenían una coincidencia con alguna otra secuencia? Para empezar, probablemente usaríamos un grep -v '#' para eliminar todas las líneas de comentarios, pero ¿entonces qué? Podríamos intentar usar wc para contar las líneas, pero solo después de eliminar también los auto-hits, donde el ID en la primera columna es igual al ID en la segunda columna. Ninguna de las utilidades que hemos visto hasta ahora —grep, sort, head o tail — puede realizar esta tarea. Necesitamos una nueva herramienta, awk, que es una herramienta de procesamiento línea por línea y columna por columna para archivos de texto: awk '<program><file>o... | awk' '<program>.

Escrita a finales de la década de 1970 y lleva el nombre de sus autores (Alfred Aho, Peter Weinberger y Brian Kernigan), awk proporciona un sofisticado lenguaje de programación que facilita el análisis de datos tabulares como los resultados de BLAST anteriores. La sintaxis para awk puede ser bastante compleja, pero gran parte de la complejidad se puede ignorar en el uso regular.

Primero, respondamos a nuestra pregunta específica, de cuántas secuencias tenían coincidencias con otras secuencias, y luego veremos alguna sintaxis awk de manera más general. El comando awk que queremos, imprimiendo solo aquellas líneas donde las dos primeras columnas no son iguales, es awk '{if (¡$1! = $2) {print $0}} '.

Desglosando el comando awk, el “programa” que se ejecuta está delimitado por las comillas simples (que cotejan su contenido en un único parámetro de línea de comandos que se envía al programa awk). El código dentro del par exterior de corchetes se ejecuta para cada línea. Por cada línea, si el contenido de la primera columna (representado por $1) no es igual a (! =) la segunda columna ($2), luego se ejecuta el código en el siguiente par de llaves, y se imprime la línea ($0) (a salida estándar).

En teoría, en este punto deberíamos poder sustituir al menos -S por un wc para contar las líneas, y así contar el número de secuencias que tenían coincidencias con otras secuencias. Desafortunadamente, en este caso, la teoría no se alinea con la realidad: inspeccionar la salida anterior revela que el ID YDR545W todavía está representado por dos líneas, por lo que esta secuencia se contaría dos veces.

¿Por qué? En el comando BLAST, solicitamos los dos primeros HSP por consulta con -max_target_seqs 2 y -max_hsps 1, así que esperábamos que el mejor HSP sería para la secuencia misma, con el segundo mejor (si existiera) para ser a un no-auto-hit. Pero en este caso, blastx decidió reportar dos no-auto-hits. De hecho, si tuviéramos que inspeccionar YDR545W, YOR396W y YLR467W, encontraríamos que sus secuencias son idénticas, por lo que BLAST necesitaba elegir dos HSP del empate de tres vías.

Para obtener el número correcto de secuencias que tenían coincidencias con otras, necesitamos eliminar los duplicados que aún se puedan encontrar en la primera columna. Podemos hacer esto agregando una especie -k1,1d -u, para una respuesta final de 2,884.

Para cualquier conjunto de datos suficientemente complejo, es una buena idea verificar tantas suposiciones al respecto como sea posible al realizar un análisis. En lo anterior, usar wc en las líneas para contar el número de secuencias que tenían aciertos a otras implicaba que en la primera columna no se listaba ningún ID dos veces. En este caso, comparar los recuentos con y sin el tipo -k1,1d -u serviría para verificar o rechazar esta suposición. En capítulos posteriores, aprenderemos más técnicas para este tipo de “comprobación de la cordura”.

Sintaxis básica para awk

Aunque awk se puede utilizar como un lenguaje de programación con todas las funciones (completo con bucles, matrices, etc.), para necesidades de programación sofisticadas, otros lenguajes como Python y R suelen ser más adecuados. Echemos un vistazo a un subconjunto práctico de su sintaxis.

Las sentencias en el bloque BEGIN se ejecutan antes de que se procese cualquier línea de la entrada, las sentencias en el bloque medio sin adornos se ejecutan para cada línea y las sentencias en el bloque END se ejecutan después de procesar la última línea. Cada uno de estos bloques es opcional, aunque el bloque medio “para cada línea” rara vez se omite en la práctica.

Al procesar una línea, hay varias variables disponibles para su uso. Aunque muchos de ellos comienzan con un $, no son variables de entorno (porque hemos colocado el “programa” dentro de comillas simples, se envían a awk como cadenas inalteradas, signos $ intactos).

• $0
  • Sostiene el contenido de toda la línea que se está procesando.
• $1, $2, etc.
  • $1 contiene el contenido de la primera columna de la línea actual, $2 contiene el contenido de la segunda columna de la línea actual, y así sucesivamente. Al igual que sort, awk por defecto asume que las columnas están separadas por espacios en blanco.
• NF
  • La variable especial NF contiene el número de columnas (también conocidas como campos) en la línea actual (awk no requiere que todas las líneas contengan el mismo número de columnas).
• NR
  • NR contiene el número de líneas que se han procesado hasta el momento, incluida la línea actual. Así, en la línea INICIO, NR tiene 0; durante el procesamiento de la primera línea, NR tiene 1; y en el bloque END, NR contiene el número total de líneas en la entrada.

Tenga en cuenta que tanto el NF como el NR carecen de un prefijo $. El texto colocado en los tres bloques puede constar de una amplia variedad de cosas, incluyendo condicionales, declaraciones impresas, manipulaciones lógicas y matemáticas, y así sucesivamente. Quizás una colección de ejemplos ayude a ilustrar la utilidad de awk. Todos estos ejemplos se basan en la salida BLAST anterior, después de filtrar las líneas de comentario con grep -v '#'.

Este comando imprime solo las dos primeras columnas de la tabla, separadas por un espacio (el valor predeterminado cuando se incluye una coma en una instrucción print):

En lugar de separar las dos columnas de salida por un espacio, podemos separarlas por una cadena como:::, produciendo solo una sola columna conglomerada de salida.

Si queremos agregar una nueva primera columna que simplemente contenga el número de línea, podemos usar la variable NR junto con la variable $0:

Las sentencias IF permiten a awk ejecutar otras sentencias condicionalmente; la sintaxis es if () {} <logical expression><statements to execute>. Además, si una columna contiene valores numéricos, awk puede trabajar con ellos como tales, y awk incluso entiende la notación científica. Aquí hay un ejemplo donde solo se imprimen líneas con valores HSP E (la décima columna en nuestro ejemplo) de menos de 1e-10.

Observe que la organización de los corchetes produce una estructura de bloques anidada; aunque para este simple caso se podría omitir el conjunto de corchetes internos, suele ser la mejor práctica incluirlos, ya que ilustran exactamente qué sentencia es controlada por el if anterior. ^[1]

Las declaraciones IF pueden controlar múltiples condiciones y, a veces, ayuda a romper los programas awk en varias líneas para ayudar con su legibilidad, especialmente cuando se incluyen en scripts ejecutables. Esta sofisticada declaración agrega una nueva primera columna que categoriza cada HSP como “genial”, “bueno” o “ok”, dependiendo del valor E, imprimiendo solo los dos ID y el valor E (columnas 1, 2 y 10):

Es bastante fácil determinar si una columna en particular es igual a una cadena dada, por ejemplo, para sacar todas las líneas donde la primera columna es YAL054C:

Los cálculos matemáticos son una buena característica de awk. Por ejemplo, las columnas 4 y 5 contienen la longitud total de la secuencia de consulta y la secuencia del sujeto, respectivamente, por lo que podríamos desear imprimir la relación de estas dos como una columna adicional al final.

Luego podríamos canalizar el resultado a una especie -k3,3g | tail -n 5 para ver los cinco HSP con las proporciones más grandes. Tenga en cuenta, sin embargo, que al realizar operaciones matemáticas o comparaciones con columnas, cualquier contenido que no pueda ser analizado como un número (1.5 puede ser, como pueden 2 y 4e-4, pero no i5 o NA) puede ser truncado (por ejemplo, 10x1 se trata como solo 10) o tratados como 0. El uso de sort on column con -g puede revelar tales problemas potenciales, ya que se usa el mismo método subyacente para analizar.

Hay una variedad de funciones matemáticas integradas en awk. Aquí hay una muestra:

• registro ()
  • Devuelve el logaritmo natural de su argumento, como en la impresión $10 * log ($3 * $4) para imprimir el log de la multiplicación de la tercera y cuarta columnas por la décima columna. ^[2]
• longitud ()
  • La función length () devuelve el número de caracteres en su argumento, como en length ($1) para la longitud de carácter de la primera columna, y length ($0) para la longitud de caracteres de toda la línea (espacios y caracteres de tabulación incluidos).
• **
  • Este operador devuelve el lado izquierdo elevado a la potencia del lado derecho, como en $1**2 para el cuadrado de la primera columna.
• %
  • El operador de módulo, devolviendo el resto después de dividir el lado izquierdo por el lado derecho. Por ejemplo, NR%4 será 1 en la primera línea, 2 en la segunda, 3 en la tercera, 0 en la cuarta, 1 en la quinta, y así sucesivamente.
• exp ()
  • Esta función devuelve su argumento planteado al poder natural e. Por ejemplo, log (exp ($1)) devuelve el valor de la primera columna.
• int ()
  • Devuelve la parte entera de su argumento. Por ejemplo, int (6.8) devuelve 6 y int (-3.6) devuelve -3.
• rand ()
  • Cuando no se le dan argumentos, devuelve un número aleatorio entre 0 (inclusive) y 1 (exclusivo). Por defecto, cada vez que se ejecuta awk, la serie de números aleatorios producidos por múltiples llamadas a rand () es la misma. Para obtener una serie aleatoria “aleatoria”, ejecute srand () (que “siembra” el generador de números aleatorios) en el bloque BEGIN, como en BEGIN {srand ()} {print rand (), $0}.

Las expresiones lógicas se pueden combinar con operadores booleanos, incluyendo && para “y” y || para “or” (lo que produce true si uno o ambos lados son verdaderos), y la agrupación se puede lograr con paréntesis. Por ejemplo, podríamos desear imprimir solo aquellas líneas donde la primera columna no es igual a la segunda, y o bien la décima columna es menor que 1e-30 o la segunda columna es YAL044C.

Hasta ahora, no hemos hecho mucho uso de los bloques BEGIN o END, que son especialmente útiles cuando definimos y actualizamos nuestras propias variables. Podemos realizar esta tarea con una = asignación (no confundir con la comparación ==). Este comando imprime los valores promedio de E en nuestro archivo de resultado BLAST de ejemplo.

Este comando funciona porque el lado derecho de una asignación a una variable con = se evalúa antes de que ocurra la asignación. Así, en el bloque BEGIN, la variable sumeval se inicializa a 0, luego para cada línea se suma el valor de sumeval al contenido de la décima columna (el valor E de esa línea), y el resultado se almacena en sumeval. Finalmente, en el bloque END, sumeval contiene la suma total de los valores E, y podemos dividir este resultado por el número de líneas procesadas, NR.

Podemos ejecutar múltiples sentencias dentro de un solo bloque si las separamos con punto y coma. En el ejemplo anterior, el valor promedio de E calculado incluye autoaciertos. Podemos filtrarlos con una sentencia if antes de modificar sumeval, pero luego no vamos a querer dividir el resultado por NR, porque eso también incluirá los recuentos de autoaciertos. Para resolver este problema, tendremos que mantener dos variables.

Como antes, algunos ID todavía están presentes más de una vez en la primera columna con esta solución, por lo que puede tener más sentido filtrar primero las líneas deseadas usando la solución awk y sort -k1,1d -u desde arriba, y luego usar otro awk para el cálculo promedio.

Ejercicios

1. En el archivo pz_blastx_yeast_top10.txt, ¿cuántos HSP (líneas) tienen un valor E que es menor que 1e-30 o tienen un valor de identidad mayor al 50%? Use awk, wc y grep si es necesario para calcular la respuesta.
2. El archivo contig_stats.txt describe estadísticas para cóntigos de un ensamblaje de genoma de novo (un cóntig es una “pieza” ensamblada del genoma). En la cuarta columna se describe el contenido de GC de los diversos cóntigos. Otros análisis indican que la mayoría de los genes correctamente ensamblados tienen una cobertura promedio (segunda columna) de entre 80.0 y 150.0. Use awk para determinar el contenido promedio de GC para cóntigos con cobertura entre 80.0 y 150.0. Luego use otra invocación de awk para determinar el contenido promedio de GC para todos los demás cóntigos. (No cuente la línea de cabecera en el archivo en sus cálculos).
3. El archivo pz.annot.txt es el resultado de un análisis de ontología génica (GO) para el conjunto completo de secuencias de ADNc de Papilio zelicaon ensambladas. Debido a la naturaleza incompleta del proceso de anotación, no a todas las secuencias se les asignaron términos GO. ¿Cuántos ID de secuencia diferentes se representan en este archivo?
4. Algunas versiones del programa de clasificación pueden ordenar líneas en orden “aleatorio” usando el indicador -R. Sin embargo, en lugar de usar esta bandera, use grep, awk (con la característica rand ()) sort (sin el indicador -R) y diríjase para seleccionar cinco ID aleatorios de Pz_CDNAS.fasta. Un resultado de ejemplo podría verse así:
   El mismo comando debería producir una lista diferente de cinco ID cada vez que se ejecuta.

---

1. Esta construcción anidada, un bloque controlado dentro de otro bloque que se ejecuta para cada elemento de un conjunto (en este caso, para cada línea), es uno de nuestros primeros ejemplos de programación! Un indicio para reducir la confusión al producir tales estructuras es rellenar su estructura desde el exterior hacia adentro, agregando pares de símbolos y luego “retrocediendo” en ellos según sea necesario. En este ejemplo, podríamos haber comenzado con awk ", y luego agregamos los corchetes para producir awk '{}', siguiente awk '{if () {}}', y finalmente rellenamos la lógica con awk '{if ($10 < 1e-10) {print $0}}'.
   
2. Si le preocupa dónde se permiten espacios en las declaraciones awk, trate de no estarlo: en su mayor parte, están permitidos en cualquier lugar, y puede sentirse libre de usarlos para mejorar la legibilidad de sus declaraciones. No están permitidos en palabras clave y variables: i f ($1 > $2) {print N R} sería una expresión no válida porque i f, $1 y N R tienen espacios erróneos.