Sesión 1. Introducción al biocómputo en sistemas GNU/Linux
Pablo Vinuesa, Centro de Ciencias Genómicas - UNAM
2024-09-02
- 1 Presentación
- 2 Navegación del sistema y operaciones
básicas - primer contacto con un sistema Linux
- 2.1 Conexión a un servidor y
exploración de sus características básicas
- 2.1.1 ssh establecer sesion remota encriptada (segura) via ssh al servidor con número dado de IP
- 2.1.2 hostname muestra el nombre del host (la máquina a la que estoy conectado) y la IP
- 2.1.3 uname muestra el sistema operativo del host
- 2.1.4 top o htop muestran los procesos en ejecución y los recursos que consumen
- 2.2 Exploración del sistema de archivos mediante la combinación de 3 comandos y 6 símbolos básicos
- 2.3 Permisos (modo de un archivo)
- 2.4 Moviéndonos por el sistema de
archivos: comando cd y haciendo uso de ruta absoluta
/full/path/to/fileor relativa../file - 2.5 Generación de directorios: comando mkdir
- 2.6 Copiar, mover, renombrar y borrar archivos con: cp, mv, rm, scp y sftp
- 2.7 Descarga de archivos desde la línea de comandos usando wget
- 2.8 Derterminar el tipo de archivo con file
- 2.9 Generación de ligas simbólicas a archivos: comando ln -s /ruta/al/archivo/fuente nombre_la_liga
- 2.10 Visualización de contenidos de archivos: comando head, tail, cat, less, more
- 2.11 Edición de archivos con los editores vim o [g|n]edit
- 2.12 Edición de archivos con el editor de flujo sed (stream editor)
- 2.13 Uso de tuberias de herramientas UNIX/Linux para filtrado de texto con cut, grep, sort, uniq, wc y |
- 2.14 redireccionado de la salida STOUT a un archivo con el comando >
- 2.15 Manual de cada comando: man command
- 2.16 Ayuda de comandos del \(Shell\): command –help
- 2.17 El poderoso comando \(find\)
- 2.18 ¿Cuánto espacio queda en disco? - comandos df y du
- 2.1 Conexión a un servidor y
exploración de sus características básicas
- 3 Inicios de programación \(Shell\) en \(Bash\)
- 3.1 Tipos, asignación y uso de variables
- 3.2 Condicionales
- 3.3 Bucles for
- 3.4 Bucles while
- 3.5 \(Bash\) scripts - primeros pasos
- 3.6 Funciones en Bash
- 3.7 Machote de un Bash script que llama a funciones
- 3.8 Alineamiento múltiple de secuencias e interconversión de formatos con align_seqs_with_clustal_or_muscle.sh y convert_alnFormats_using_clustalw.sh
- 4 El lenguaje de procesamiento de patrones AWK
- 5 Ejercicios integrativos de uso de herramientas de filtrado del Shell
- 6 Reto de programación - ejercicio de parseo de archivos FASTA
- 7 Programando un pipeline en \(Bash\)
- 8 Consideraciones finales y referencias recomendadas para continuar aprendiendo programación \(Shell\)
1 Presentación
Este apunte fue creado para el Taller de Pagenómica y Filogenómica Microbiana de los Talleres Internacionales de Bioinformática Talleres NNB, celebrados en el Centro de Ciencias Genómicas de la Universidad Nacional Autónoma de México, del 22 al 26 de enero de 2024 por Pablo Vinuesa, CCG-UNAM @pvinmex.
Version: 2024-09-02
Si éste es tu primer contacto con Linux, te recomiendo que leas primero esta presentación sobre introducción al biocómputo en sistemas Linux - PDF.
1.1 Licencia y términos de uso
El material del Taller de Pagenómica y Filogenómica Microbiana lo distribuyo públicamente a través de este repositorio GitHub bajo la Licencia No Comercial Creative Commons 4.0
This
work is licensed under a
Creative
Commons Attribution-NonCommercial 4.0
1.2 Referencias adicionales
Una vez que domines los comandos básicos que se presentarán seguidamente, recomiendo revisar tutoriales mucho más detallados y completos como los siguientes:
- The Linux Command Line - a complete introduction. William E. Shotts, Jr. No Starch Press
- Bash Reference Manual
- Advanced Bash Scripting Guide
- Bioinformatics Data Skills: Reproducible and Robust Research with Open Source Tools. Vince Buffalo. O’Reilly Media 2014
- Taller de introducción al biocómputo en sistemas UNIX y GNU/Linux para filoinformática
3 Inicios de programación \(Shell\) en \(Bash\)
Vamos a aprender algunos de los elementos sintácticos básicos de programación Shell, implementados en el lenguaje \(Bash\).
3.1 Tipos, asignación y uso de variables
El primer paso es la asignación de variables. Todos los lenguajes hacen uso de ellas. Se usan para guardar valores sencillos (variables escalares) numéricos o de cadenas de caracteres, o listas y estructuras de datos más complejas, como \(arreglos\) indexados por índices posicionales, o \(hashes\), que son arreglos indexados por llaves especificas.
Independientemente del tipo de variable, los nombres que les damos deben inciar con un caracter alfabético o guión bajo, y el resto de caracteres deben ser alfanuméricos, es decir de la clase[a-zA-Z0-9_]
y no contener espacios.
Las variables de \(Shell\) (\(Bash\) inclusive) son globales y no necesitan ser declaradas (salvo \(hashes\)).
Según su uso, las variables temporales que se usan como aliases temporales por ejemplo en bucles, suelen ser nombradas con una sola letra, por conveniencia y por ser obvio lo que contienen
for f in *.fna; do echo $f; donePara otras variables que se usan en diferentes puntos de un programa, conviene hacer uso de nombres cortos pero informativos, como
formato_de_entrada=''
formato_de_salida=''
runmode=''Es una buena práctica mantener un estilo consistente para los nombres de variables, evitando usar variables en puras mayúsculas, para minimizar el riesgo de interferencia con variables de ambiente, como USER, HOME, SHELL, PATH, que por convención se nombran en mayúsculas.
3.1.1 Variables escalares
La sintaxis básica de asignación de un valor simple a una variable escalar y uso de comillas
varName=VALUE- para recuperar el valor de una variable, le añadimos el
prefijo $. Para imprimir el valor asignado a la variable,
usamos
echo $varName
archivo_de_comandos_linux=linux_commands.tab
echo "$archivo_de_comandos_linux"
# asignación de cadena de texto, con espacios y otros símbolos, entre comillas sencillas
var2='Cadena con espacios, entre comillas sencillas'
echo "var2: $var2"
# asignación de cadena de texto, con espacios y otros símbolos, incluyendo variables (en este caso de ambiente) entre comillas dobles, para interpolación de variables
saludo_inicial="Bienvenido a $HOSTNAME, $USER! Te recuerdo que hoy es $(date)"
echo "$saludo_inicial"
# Llamado a variable entre comillas sencillas NO INTERPOLA!!!
echo ' >>> Ojo, VARIABLE ENTRE COMILLAS SENCILLAS NO INTERPOLA: $saludo_inicial; se imprime literalmente'
echo## linux_commands.tab
## var2: Cadena con espacios, entre comillas sencillas
## Bienvenido a puma, vinuesa! Te recuerdo que hoy es lun 02 sep 2024 14:31:17 CST
## >>> Ojo, VARIABLE ENTRE COMILLAS SENCILLAS NO INTERPOLA: $saludo_inicial; se imprime literalmenteNo es necesario hacer llamados a las variables entre comillas dobles, pero es una buena práctica, como se muestra en los ejemplos del bloque anterior.
3.1.2 Captura de la salida de un comando en una variable con var=$(comando)
wkdir=$(pwd)
date=$(date | awk '{print $3,$2,$6}' | sed 's/ //g')
h=$HOSTNAME
echo ">>> working in: $wkdir at <$h> on <$date>"## >>> working in: /home/vinuesa/Cursos/intro_biocomp4genomics_TFC/sesion1_intro2linux at <puma> on <sep02CST>3.1.3 Modificación de variables y operaciones con ellas
wkdir=$(pwd)
echo "wkdir: $wkdir"
# 1. cortemos caracteres por la izquierda (todos los caracteres por la izquierda, hasta llegar a último /)
basedir=${wkdir##*/}
echo "basedir: $basedir # \${wkdir##*/}"
# 2. cortemos caracteres por la derecha (cualqier caracter hasta llegar a /)
echo "path to basedir: ${wkdir%/*} # \${wkdir%/*}"
# 3. contar el número de caracteres (longitud) de la variable
echo "basedir has ${#basedir} characters # \${#basedir}"
## wkdir: /home/vinuesa/Cursos/intro_biocomp4genomics_TFC/sesion1_intro2linux
## basedir: sesion1_intro2linux # ${wkdir##*/}
## path to basedir: /home/vinuesa/Cursos/intro_biocomp4genomics_TFC # ${wkdir%/*}
## basedir has 19 characters # ${#basedir}3.2 Condicionales
- La sintaxis básica de un condicional simple, llamado directamente desde la línea de comando (en formato de una línea) es así:
if [ condición ]; then comando1; comando2 ...; fi
- también hay una versión más corta para tests simples:
[ condición ] && comando1 &&
comando2
- En un script generalmente escibimos los condicionales como un bloque indentado, para mejor legibilidad, asi como el que sigue:
if [ condición ]; then
comando1
comando2
fiNoten que las líneas de comandos (comando1, comando2) no necesitan terminarse con ‘;’
3.2.1 Comparación de íntegros en condicionales
i=5
j=3
if [ "$i" -lt "$j" ]; then
echo "$i < $j"
elif [ "$i" -gt "$j" ]; then
echo "$i > $j "
fi## 5 > 33.2.2 Comparación de cadenas de caracteres en condicionales
c=carla
j=juan
if [ "$c" == "$j" ]; then
echo "$c = $j"
elif [ "$c" != "$j" ]; then
echo "c:$c != j:$j "
fi## c:carla != j:juan3.2.3 Comprobación de la existencia ([ -e file ]) de un archivo de tamaño > 0 bytes
touch empty_file
ls -l empty_file
ls -l assembly*gz
f=$(ls assembly*gz)
if [ -e empty_file ]; then
echo "empty_file file exists"
fi
if [ ! -s empty_file ]; then
echo "empty_file file exists but is empty"
fi
if [ -s "$f" ]; then
size=$(du -h assembly_summary.txt.gz | cut -f1)
# o tambien
# size=$(ls -hs assembly_summary.txt.gz | cut -d' ' -f1)
echo "$f exists and has size: $size"
fi
## -rw-rw-r-- 1 vinuesa vinuesa 0 sep 2 14:31 empty_file
## -rw-r--r-- 1 vinuesa vinuesa 6780296 ago 17 2023 assembly_summary.txt.gz
## empty_file file exists
## empty_file file exists but is empty
## assembly_summary.txt.gz exists and has size: 6.5M3.2.4 La versión corta de test [ condición ] && comando1 && comando2 …
# también podemos usar la versión corta del test:
f=$(ls *gz) # <<< peligroso: usar sólo si estás seguro que existe un solo archivo *gz en el directorio
[ -s "$f" ] && echo "$f exists and is non-empty"## assembly_summary.txt.gz exists and is non-empty3.2.5 if; elif; else
if [[ "$OSTYPE" == "linux-gnu" ]]
then
OS='linux'
no_cores=$(awk '/^processor/{n+=1}END{print n}' /proc/cpuinfo)
host=$(hostname)
echo "running on $host under $OS with $no_cores cores :)"
elif [[ "$OSTYPE" == "darwin"* ]]
then
OS='darwin'
no_cores=$(sysctl -n hw.ncpu)
host=$(hostname)
echo "running on $host under $OS with $no_cores cores :)"
else
OS='windows'
echo "oh no! another windows box :( ... you should better change to linux :) "
fi
## running on puma under linux with 4 cores :)3.3 Bucles for
la sintaxis general de un bucle for en \(Bash\) es:
for ALIAS in LIST; do CMD1; CMD2; done
donde el usuario tiene que cambiar los términos en mayúsculas por opciones concretas. ALIAS es el nombre de una variable temporal a la que se asigna secuencialmente cada valor de la lista LIST. Cómo generar la lista LIST dependerá de la situación.
Así por ejemplo, si tuviéramos muchos archivos de secuencias homólogas con la extensión *.faa en un directorio, podríamos alinearlas secuencialmente con \(clustalo\) usando un comando como el siguiente:
for file in *.faa; do clustalo -i $file -o ${file%.*}_cluoAln.faa; done
donde ALIAS=file, LIST=*.faa y CMD1 es una llamada al programa de alineamientos múltiples \(clustalo\) que veremos más adelante en este taller.
3.3.1 Ejemplo de bucle for, acoplado a las herramientas de filtrado y de manipulación de variables
La idea del ejercicio es generar archivos a partir de linux_basic_commands.tab que contengan sólo los comandos de cada clase, nombrando a los archivo resultantes con el valor de dicha clase, almacenados en la segunda columna de la tabla
# veamos la cabecera y cola del archivo linux_basic_commands.tab
head linux_basic_commands.tab
echo '--------------------------------------------------'
tail linux_basic_commands.tab## IEEE Std 1003.1-2008 utilities Name Category Description First appeared
## admin SCCS Create and administer SCCS files PWB UNIX
## alias Misc Define or display aliases
## ar Misc Create and maintain library archives Version 1 AT&T UNIX
## asa Text processing Interpret carriage-control characters System V
## at Process management Execute commands at a later time Version 7 AT&T UNIX
## awk Text processing Pattern scanning and processing language Version 7 AT&T UNIX
## basename Filesystem Return non-directory portion of a pathname; see also dirname Version 7 AT&T UNIX
## batch Process management Schedule commands to be executed in a batch queue
## bc Misc Arbitrary-precision arithmetic language Version 6 AT&T UNIX
## --------------------------------------------------
## val SCCS Validate SCCS files System III
## vi Text processing Screen-oriented (visual) display editor 1BSD
## wait Process management Await process completion Version 4 AT&T UNIX
## wc Text processing Line, word and byte or character count Version 1 AT&T UNIX
## what SCCS Identify SCCS files PWB UNIX
## who System administration Display who is on the system Version 1 AT&T UNIX
## write Misc Write to another user's terminal Version 1 AT&T UNIX
## xargs Shell programming Construct argument lists and invoke utility PWB UNIX
## yacc C programming Yet another compiler compiler PWB UNIX
## zcat Text processing Expand and concatenate data 4.3BSDAntes de correr el bucle, lista los archivos en el directorio de trabajo
# veamos el contenido del directorio antes de correr el bucle
lsAhora el bucle. En este caso ALIAS=type y LIST corresponde a la lista
de valores únicos almacenados en la segunda columna de la tabla:
$(cut -f2 linux_basic_commands.tab | sort -u)
#>>> Ejemplo integrativo: usa un bucle for, acoplado a las herramientas de filtrado arriba mostradas,
# para generar archivos que contengan solo los comandos de las diferentes categorias
# nombrando a los archivos por estas
# for type in $(cut -f2 linux_basic_commands.tab | sort -u); do grep "$type" linux_basic_commands.tab > ${type}.cmds; done
for type in $(cut -f2 linux_basic_commands.tab | cut -d_ -f1 | sort -u | grep -v Category); do
grep -w "$type" linux_basic_commands.tab > ${type}.cmds
done
Y voilà:
# veamos el contenido del directorio después de correr el bucle
ls *.cmds## administration.cmds
## Batch.cmds
## C.cmds
## Filesystem.cmds
## FORTRAN77.cmds
## management.cmds
## Misc.cmds
## Network.cmds
## Process.cmds
## processing.cmds
## programming.cmds
## Programming.cmds
## SCCS.cmds
## Shell.cmds
## Std.cmds
## System.cmds
## Text.cmds
## utilities.cmds# veamos el contenido de uno de los nuevos archivos generados
cat programming.cmds## cc/c99 C programming Compile standard C programs IEEE Std 1003.1-2001
## cflow C programming Generate a C-language call graph System V
## command Shell programming Execute a simple command
## ctags C programming Create a tags file 3BSD
## cxref C programming Generate a C-language program cross-reference table System V
## echo Shell programming Write arguments to standard output Version 2 AT&T UNIX
## expr Shell programming Evaluate arguments as an expression Version 7 AT&T UNIX
## false Shell programming Return false value Version 7 AT&T UNIX
## fort77 FORTRAN77 programming FORTRAN compiler XPG4
## getopts Shell programming Parse utility options
## lex C programming Generate programs for lexical tasks Version 7 AT&T UNIX
## logger Shell programming Log messages 4.3BSD
## nm C programming Write the name list of an object file Version 1 AT&T UNIX
## printf Shell programming Write formatted output 4.3BSD-Reno
## read Shell programming Read a line from standard input
## sh Shell programming Shell, the standard command language interpreter Version 7 AT&T UNIX (in earlier versions, sh was either the Thompson shell or the PWB shell)
## sleep Shell programming Suspend execution for an interval Version 4 AT&T UNIX
## strings C programming Find printable strings in files 2BSD
## strip C programming Remove unnecessary information from executable files Version 1 AT&T UNIX
## tee Shell programming Duplicate the standard output Version 5 AT&T UNIX
## test Shell programming Evaluate expression Version 7 AT&T UNIX
## true Shell programming Return true value Version 7 AT&T UNIX
## xargs Shell programming Construct argument lists and invoke utility PWB UNIX
## yacc C programming Yet another compiler compiler PWB UNIX# finalmente borremos los nuevos archivos generados
rm *.cmds3.4 Bucles while
la sintaxis general de un bucle while en \(Bash\) es:
while <CRITERIOS>; do <BLOQUE>; done
donde el usuario tiene que cambiar los términos en <MAYÚSCULAS> por opciones concretas. El <BLOQUE> de comandos se ejecutará mientras que el|los <CRITERIOS> se cumplan, es decir, que tengan un estatus de salida (exit status) igual a 0.
Van unos ejemplos genéricos:
# cuenta de 1..3; noten el uso del operador (( )), usado para operaciones aritméticas; ((count++)) autoincrementa la variable count
count=1
while (( count <= 3)); do
echo $count
((count++))
done## 1
## 2
## 33.4.1 Procesamiento dea archivos con bucles while read
Los bucles \(while\) se usan principalmente para procesar archivos. De hecho, la manera más robusta y canónica de procesar archivos en \(BASH\) es mediante el uso de bucles \(while\).
Reescribamos haciendo uso de un bucle \(while\) el ejemplo anterior que usaba un bucle for para procesar el archivo linux_basic_commands.tab.
El bucle \(while\) procesará línea a línea el contenido del archivo linux_basic_commands.tab
Con read -r asignamos a las variables name
category descr y appeared el contenido de cada uno
de los cuatro campos de la tabla linux_basic_commands.tab, leída
por el bucle con done < file. Cada línea es cortada por
espacios, tabuladores o saltos de línea, acorde a la variable de
ambiente \(\$IFS\).
#>>> Ejemplo integrativo: usa un bucle while, acoplado a grep para generar archivos que contengan solo los comandos de las diferentes categorias
# nombrando a los archivos por estas
while read -r name category descr appeared; do
# La variable category contendrá sólo la primera palabra del campo correspondiente.
grep -w "$category" linux_basic_commands.tab > "${category}".cmds
done < linux_basic_commands.tabY voilà:
# veamos el contenido del directorio después de correr el bucle
ls *.cmds## Batch.cmds
## C.cmds
## Filesystem.cmds
## FORTRAN77.cmds
## Misc.cmds
## Network.cmds
## Process.cmds
## Programming.cmds
## SCCS.cmds
## Shell.cmds
## Std.cmds
## System.cmds
## Text.cmds# veamos el contenido de uno de los nuevos archivos generados
cat Shell.cmds## command Shell programming Execute a simple command
## echo Shell programming Write arguments to standard output Version 2 AT&T UNIX
## expr Shell programming Evaluate arguments as an expression Version 7 AT&T UNIX
## false Shell programming Return false value Version 7 AT&T UNIX
## getopts Shell programming Parse utility options
## logger Shell programming Log messages 4.3BSD
## printf Shell programming Write formatted output 4.3BSD-Reno
## read Shell programming Read a line from standard input
## sh Shell programming Shell, the standard command language interpreter Version 7 AT&T UNIX (in earlier versions, sh was either the Thompson shell or the PWB shell)
## sleep Shell programming Suspend execution for an interval Version 4 AT&T UNIX
## tee Shell programming Duplicate the standard output Version 5 AT&T UNIX
## test Shell programming Evaluate expression Version 7 AT&T UNIX
## true Shell programming Return true value Version 7 AT&T UNIX
## xargs Shell programming Construct argument lists and invoke utility PWB UNIX# finalmente borremos los nuevos archivos generados
rm *.cmds3.4.1.1 Eficiencia de programación: optimización de bucles para evitar trabajo redundante
El bucle \(while\) mostrado arriba, tal y como está programado, es altamente ineficiente. ¿Porqué?
Revisa el código y el formato de la tabla. Recuerda: \(read\) lee línea a línea los contenidos de la tabla y ejecuta \(grep\) en cada una de ellas. Por tanto la llamada a \(grep\) por cada línea de la tabla es redundante y sobreescribe el archivo *.cmds correspondiente por cada nueva instancia que encuentra de “category”. Esto es claramente ineficiente, desperdiciando recursos de cómputo.
Un código eficiente debe ejecutar \(grep\) y escribir el archivo correspondiente sólo con instancias únicas del patrón guardado en la variable $category, tal y como hicimos en el ejemplo del bucle \(for\) mostrado arriba.
3.4.1.1.1 Optimización
del bucle while haciendo uso de \(sustitución\
de\ procesos\): <(comando)
Una opción muy conveniente es hacer uso de \(sustitución\ de\ procesos\) mediante el uso
de la siguiente construcción sintáctica: <(comando). La
salida del comando en <(comando) es interpretada como si
fuera un archivo. Esto lo podemos utilizar para alimentar al bucle \(wile\) sólo con instancias únicas de la
clase “categoría” recuperadas a partir de la columna #2 de la tabla,
reescribiendo el código como sigue:
# generación de una lísta de categorías no redundantes (únicas) a partir de las listadas en la columna 2 de la tabla
sed '1d' linux_basic_commands.tab | cut -f 2 | cut -d' ' -f1 | sort -u## Batch
## C
## Filesystem
## FORTRAN77
## Misc
## Network
## Process
## Programming
## SCCS
## Shell
## System
## TextUna vez resuelto el problema de cómo generar una lísta de categorías
no redundantes a partir de las listadas en la columna #2 de la tabla,
hagamos uso de \(sustitución\ de\
procesos\) agregando el código mostrado arriba a nuestro bucle
\(while\) mediante la sintaxis
<(comando). Recuerda, el bucle consumirá la salida de
<(comando) como si fuera un archivo.
#>>> Ejemplo de bucle while optimizado: uso de "sustitución de procesos" con la sintaxis '<(comando)' para leer sólo instancias únicas de categorías
while read -r category; do
# La variable category contiene las instancias únicas de las categorías listadas en la columna #2 de la tabla
grep -w "$category" linux_basic_commands.tab > "${category}".cmds
done < <( sed '1d' linux_basic_commands.tab | cut -f 2 | cut -d' ' -f1 | sort -u)3.4.1.1.2 Optimización del bucle while mediante uso de \(hashes\) y \(continue\)
Otra opción para evitar el trabajo redundante en el bucle \(while\) original es el uso combinado de \(hashes\) acoplado a la estructura de control \(continue\) específica de bucles \(for\) y \(while\).
Los hashes (también conocidos en otros lenguajes como diccionarios o arreglos asociativos) son estructuras de datos muy útiles en programación que permiten establecer asociaciones entre un valor único conocido como llave y otra variable.
Un hash de BASH tiene la siguiente estructura genérica:mi_hash=( [llave1]=valor1 [llave2]=valor2 )
Para crear un hash o arreglo asociativo en BASH es
neceasrio declararlo como tal mediante declare -A
mi_hash.
# declaración del hash mi_hash
declare -A mi_hashEn el siguiente ejemplo haremos uso del hash ‘seen’ (visto) para guardar como llaves los valores únicos del primer campo de la columna 2 (Category) de la tabla linux_basic_commands.tab procesada línea a línea por el bucle.
Haremos uso de la sintaxis (( seen[“$category”]++ ))
para autoincrementar el valor asociado a la llave cada vez que
encontremos una nueva línea que contenga una instancia de una Categoría
particular.
#>>> Ejemplo de bucle while optimizado: uso de hashes y continue
# declaramos el hash seen
declare -A seen
# Noten que abajo usamos un nombre arbitrario de variable '_' para almacenar en ella
# lo que no nos interesa de la línea que viene después del segundo campo (category)
while read -r name category _; do
# La variable category contendrá sólo la primera palabra del campo correspondiente
# (recuerda, por defecto read delimita campos por espacio, tabulador y salto de línea).
# y la guardamos como llave del hash seen, autoincrementando su cuenta.
# Los dobles paréntesis de requiren para conferir del contexto de operación aritmética
(( seen["$category"]++ ))
# comprobamos que sólo hayamos visto una instancia de $category;
# si son más, continue ignorará la líne actual y continua con la siguiente
# Noten que para recuperar el valor asociado a la llave en el hash, usamos la sintaxis ${mi_hash[mi_llave]}.
# Usamos los dobles paréntesis para hacer un test de comparación entre íntegros.
# - Si hemos 'visto' category más de una vez, pasamos a la siguiente iteración del bucle,
# evitando correr grep y escribir a disco de manera redundante
(( ${seen["$category"]} > 1 )) && continue
# grep recibe instancias únicas de $category
grep -w "$category" linux_basic_commands.tab > "${category}".cmds
done < linux_basic_commands.tab
#---------------------------------------
# Impresión de contenido del hash 'seen'
#---------------------------------------
# Iteraremos sobre las llaves del hash mediante un bucle for
# para imprimir el valor (conteo) asociado al mismo
echo "# Impresión de contenido del hash 'seen: categorías y conteo de sus ocurrencias'"
for key in "${!seen[@]}"; do
echo "$key => ${seen[$key]}"
done## # Impresión de contenido del hash 'seen: categorías y conteo de sus ocurrencias'
## Programming => 1
## Shell => 14
## FORTRAN77 => 1
## Std => 1
## SCCS => 10
## Process => 14
## C => 9
## Text => 29
## System => 1
## Batch => 11
## Misc => 38
## Network => 4
## Filesystem => 28La programación del bucle \(while\) mostrado arriba es eficiente al ejecutar la llamada a \(grep\) y escritura a disco del archivo correspondiente sólo una vez para cada instancia (única) de categoría.
El bucle for al término del bloque de código mostrado arriba nos permite ver exactamente cuántas operaciones de llamada a \(grep\) y escritura a disco nos hemos ahorrado, ya que nos muestra el conteo de cada ‘categoría’ usada como llave del hash. Vean en el código del bulce \(for\) mostrado arriba cómo imprimir el contenido del hash \(seen\) iterando sobre las llaves del hash.
…
3.5 \(Bash\) scripts - primeros pasos
Para finalizar esta sección, voy a mostrarles varios ejemplos muy sencillos de Shell scripts escritos en \(Bash\) que integran varios aspectos de la sintaxis básica del lenguaje arriba descritas, así como extensiones adicionales como son argumentos posicionales y funciones
3.5.1 la “shebang line” (#!/usr/bin/env bash) y permisos de ejecución - el script ls_dir
Se conocen como scripts a archivos de texto plano (codificación ASCII) que contienen los comandos a ser ejecutados secuencialmente por un intérprete de comandos particular, como \(bash\), \(perl\), \(R\) …
Les muestro abajo el código del script ls_dir.
#!/usr/bin/env bash
# the 1st line in a script is the so-called shebang line
# which indicates the system which command interpreter
# should read and execute the following code, bash in this case
# The shebang line shown above, is a portable version for bash scripts
# AUTHOR: Pablo Vinuesa
# AIM: learning basic BASH-programming constructs for TIB-filoinfo
# https://github.com/vinuesa/TIB-filoinfo
for file in *
do
if [ -d "$file" ]
then
echo "$file"
fi
donePuedes copiar el código a un archivo que vas a nombrar como ls_dir. Usa para ello el comando \(cat\), de la manera abajo indicada:
cat > ls_dir
PEGA AQUÍ EL CÓDIGO ARRIBA MOSTRADO
CTRL-DEl script ls_dir busca los directorios presentes en el
directorio actual usando un bucle \(for\) para analizar cada archivo econtrado
por \(ls\), evaluando seguidamente si
el archivo en turno es un directorio con un condicional if [ -d
“$file” ]
Noten que la primera línea inicia con lo que se conoce como una shebang line:
#!/usr/bin/env bash
Esta línea, en la cabecera del archivo (¡sin dejar especios a la izquierda o arriba de la línea!), le indica al sistema operativo qué intérprete de comandos usar para la ejecución del código que sigue. En este caso se llama al intérprete de comandos \(bash\)
Noten también que cualquier línea que inicie con #, después de la shebang line, representa un comentario, es decir, que el texto que sigue al gato no es interpretado por \(bash\)
Para ejecutar el script como si fuera un comando cualquiera de \(Linux\), le otorgamos permisos de ejecución:
chmod 755 ls_dirComprueba los permisos:
$ ls -l ls_dir
-rwxr-xr-x 1 vinuesa vinuesa 493 ago 11 18:37 ls_dirComo puedes ver el usuario, el grupo al que pertenece y todos los demás pueden ejecutarlo (\(x\))
Finalmente copia o mueve el script a un directorio que esté en el \(PATH\), típicamente a \(HOME/bin\)
mv ls_dir ~/biny ya puedes usarlo como un comando cualquiera del sistema:
ls_dir## bin
## data
## MobaXterm_screen_shots
## pics
## recA_seq_data
## tareas
## TIB2019-T3
## TIB-filo
## working_with_linux_commands_longNota: el script ls_dir funciona perfectamente, pero no es la manera más eficiente de buscar directorios. La función \(find\) de \(Linux\) es mucho más eficiente.
3.5.2 Paso de argumentos posicionales a un script o función - el script find_dir
Los scripts (o funciones) pueden recibir opciones dadas por los usuarios para comportarse acorde a ellas. Esto les da gran flexibilidad.
La manera más sencialla de pasarle opciones al script es mediante argumentos posicionales, como los usados por find_dir. El nombre de argumentos posicionales alude al hecho que el script of función recibe los argumentos en el orden en el que listan en la línea de comandos, o en el orden en el que el script que llama a una función, le pasa los argumentos.
Detro del script o función, los argumentos posicionales quedan guardados en las variables$1, $2 ... $9
, respectivamente.
Los scripts y funciones pueden recibir opciones dadas por los usuarios para comportarse acorde a ellas. Esto les confiere gran versatilidad.
La sintaxis de paso de argumentos (arg) a un script o función es muy sencilla:
mi_script.sh arg1 arg2 ... funcion_x arg1 arg2
Veamos un ejemplo real de paso de opciones a un script o a una función es mediante argumentos posicionales, como los usados por find_dir
#!/usr/bin/env bash
# The 1st line in a script is the so-called shebang line
# which indicates the system which command interpreter
# should read and execute the following code, bash in this case
# The shebang line shown above, is a portable version for bash scripts
# AUTHOR: Pablo Vinuesa
# AIM: learning basic BASH-programming constructs for TIB2019-T3
# https://github.com/vinuesa/TIB-filoinfo
# global variables
progname=${0##*/} # find_dir; elimina la ruta absoluta que precede a la localización del script
VERSION='0.1_29Jul19'
# Function definition
function print_help()
{
# this is a so-called HERE-DOC, to easily print out formatted text
cat << HELP
$progname v$VERSION usage synopsis:
$progname <int [maximal search depth for (sub)directories, e.g.:1>
AIM: find directories below the current one with the desired max_depth
USAGE EXAMPLES:
$progname 1
$progname 2
HELP
exit 1
}
# capture user-provided arguments in variables $1, $2 ... $9
# and save them to named variables, for better readability
max_depth=$1
type=${2:-d} # if not provided, the second argument is set to 'd' by default.
# check arguments
[ -z $max_depth ] && print_help
# execute find command with the provided arguments
find . -maxdepth $max_depth -type $type$HOME/bin
, como se mostró en el ejemplo anterior.
chmod 755 find_dir && mv find_dir $HOME/binfind_dir hace la siguiente llamada a \(find\)
find . -maxdepth $max_depth -type $type
para encontrar directorios y subdirectorios por debajo del acutal, al nivel de profundidad indicado por el usuario, quien pasa un íntegro al script como único argumento
find_dir introduce la función print_help(),
function print_help(){
YOUR CODE GOES INHERE ...
}que simplemente imprime un mensaje de ayuda si es llamado sin argumentos
[ -z $max_depth ] && print_help
Probemos el script:
- primero llamando al script sin argumento:
find_dir
find_dir v0.1_29Jul19 usage synopsis:
find_dir <int [maximal search depth for (sub)directories, e.g.:1>
AIM: find directories below the current one at the desired max_depth
USAGE EXAMPLES
find_dir 1
find_dir 2
- ahora pasándole un argumento:
find_dir 13.6 Funciones en Bash
A medida que los programas se hacen más grandes y complejos, se vuelven más difíciles de diseñar, escribir y mantener.
Escribir funciones permite reducir la extensión, redundancia y complejidad de un programa complejo, facilitando su mantenimiento.
Muchos programas tienen que ejecutar múltiples veces una misma acción, como por ejemplo verificar que se ha escrito a disco un archivo con resultados generados por el programa. Escribir una función que verifique la existencia de un archivo pasado a la misma como argumento permite reducir redundancia y extensión del código, ya que en el lugar indicado del script principal, se llama a la función en vez de repetir el cógido encapsulado en la misma.
Escribir funciones además ayuda al programador a mantener código, ya que si hubiera un error en la función, sólo hay que corregirlo en un bloque de código.
Por tanto una práctica fundamental en programación es la modularización del cógido en funciones y librerías de funciones relacionadas.
La sintaxis básica de una función es la siguiente:
function funname() {
command1
command2
...
}
Una buena función debe estar especializada en realizar una sola tarea, eso si, bajo diferentes condiciones u opciones, que se le pasan como argumentos posicionales.
Las funciones típicamente deben ir al inicio del script principal, antes de que sean llamadas por el mismo. También se pueden colectar conjuntos de funciones relacionadas en un archivo o librería de funciones que se llaman al inicio del script principal mediante el comando \(source\), antes de que éste haga llamadas a las funciones, como se muestra en el “machote” de script que se presenta más adelante.
3.6.1 Llamada a la función \(print\_numbered\_table\_header\_fields\) con uno o dos argumentos posicionales
Si usas con frecuencia algunas funciones, puedes ponerlas también el archivo de inicialización \(\$HOME/.bashrc\) para sesiones locales en tu máquina, o en \(\$HOME/.bash\_profile\), en una máquina remota. De esta manera, cada vez que inicias una sesión local o remota, las funciones serán exportadas al ambiente y las podrás usar como si fuera cualquier otro comando de Linux.
Como ejemplo, veamo la función \(print\_numbered\_table\_header\_fields\), que tengo en mi \(\$HOME/.bashrc\).
Los argumentos posicionales
$1, $2 …son capturados en la función como variables localizadas o locales usandolocal. Localizar una variable a una función quiere decir que sólo es visible denro de ésta. Esto es importante para evitar que interfieran con otras variables definidas en el script principalEl número total de argumentos recibidos queda guardado en la variable
$#Examina la función ¿Qué crees que hace?
function print_numbered_table_header_fields()
{
#: AUTHOR: Pablo Vinuesa, @pvinmex, CCG-UNAM
# provide table name to parse (tsv format expected; can skip a certain number of comment lines)
local table=$1
local skip_lines=$2
[ $# -lt 1 ] && echo "$FUNCNAME usage: <table_name> [<number_of_top_comment_lines_to_skip>]"
if [ $# -eq 1 ]; then
head -1 "$table" | sed 's/\t/\n/g' | nl
elif [ $# -eq 2 ]; then
tail -n +"$((skip_lines+1))" "$table" | sed 's/\t/\n/g' | nl
fi
}Copia la función al archivo \(\$HOME/.bash\_profile\) en buluc y ejecuta
source $HOME/.bash_profilepara releer el archivo y que se exporte la función al ambiente.llamado a la función sin argumentos, imprime la ayuda
print_numbered_table_header_fieldsprint_numbered_table_header_fields usage: <table_name> [<number_of_top_comment_lines_to_skip>]- llamado a la función con el nombre de la tabla a procesar como único argumento
print_numbered_table_header_fields linux_basic_commands.tab
1 IEEE Std 1003.1-2008 utilities Name
2 Category
3 Description
4 First appeared3.7 Machote de un Bash script que llama a funciones
El siguiente bloque muestra un machote para un \(script\) básico, que llama a funciones y que recibe argumentos posicionales desde la línea de comandos. Puedes usarlo para facilitarte la escritura de tus propios scripts.
#!/usr/bin/env bash #: PROGRAM:#: AUTHOR: # #: PROJECT START: # #: AIM: progname=${0##*/} # vers='0.1' # GLOBALS #DATEFORMAT_SHORT="%d%b%y" #TIMESTAMP_SHORT=$(date +${DATEFORMAT_SHORT}) #date_F=$(date +%F |sed 's/-/_/g')- #date_T=$(date +%T |sed 's/:/./g') #start_time="$date_F$date_T" #current_year=$(date +%Y) #----------------------------------------------------------------------------# #>>>>>>>>>>>>>>>>>>>>>>>>>>>> FUNCTION DEFINITIONS <<<<<<<<<<<<<<<<<<<<<<<<<<# #----------------------------------------------------------------------------# # import functions defined in the aux_fun_lib file source $HOME/bin/lib/aux_fun_lib function check_dependencies() { for programname in prog1 prog2 do bin=$(type -P $programname) if [ -z "$bin" ]; then echo echo "# ERROR: $programname not in place!" echo "# ... you will need to install \"$programname\" first or include it in \$PATH" echo "# ... exiting" exit 1 fi done echo echo '# Run check_dependencies() ... looks good: all required binaries and perl scripts are in place.' echo } #----------------------------------------------------------------------------------------- function print_help() { cat << HELP $progname v$vers usage synopsis: $progname AIM: OUTPUT: HELP check_dependencies exit 0 } #----------------------------------------------------------------------------------------- function check_output() { #>>> set color in bash # SEE: echo http://stackoverflow.com/questions/5947742/how-to-change-the-output-color-of-echo-in-linux # And very detailed: http://misc.flogisoft.com/bash/tip_colors_and_formatting # ANSI escape codes # Black 0;30 Dark Gray 1;30 # Red 0;31 Light Red 1;31 # Green 0;32 Light Green 1;32 # Brown/Orange 0;33 Yellow 1;33 # Blue 0;34 Light Blue 1;34 # Purple 0;35 Light Purple 1;35 # Cyan 0;36 Light Cyan 1;36 # Light Gray 0;37 White 1;37 RED='\033[0;31m' GREEN='\033[0;32m' #YELLOW='\033[1;33m' #BLUE='\033[0;34m' #CYAN='\033[0;36m' NC='\033[0m' # No Color => end color #printf "I ${RED}love${NC} ${GREEN}Stack Overflow${NC}\n" outfile=$1 if [ -s "$outfile" ] then echo -e "${GREEN} >>> wrote file $outfile ...${NC}" else echo echo -e "${RED} >>> ERROR! The expected output file $outfile was not produced, will exit now!${NC}" echo exit 1 fi } #----------------------------------------------------------------------------------------- # # >>>> MAIN CODE <<<< # var1=$1 var2=${2:-10} [ -z "$var1" ] && print_help cat << PARAMS $progname $var1 $var2 PARAMS prog1 $var1 $var2
3.8 Alineamiento múltiple de secuencias e interconversión de formatos con align_seqs_with_clustal_or_muscle.sh y convert_alnFormats_using_clustalw.sh
Para finalizar esta sección de introducción a la programación en \(Bash\), explora el código del script align_seqs_with_clustal_or_muscle.sh, que permite hacer alineamientos múltiples con \(clustalw\) o con \(muscle\), dos programas muy populares para este fin.
Explora también el código del script convert_alnFormats_using_clustalw.sh, que permite hacer inteconversión de formatos con clustalw.
Los alineamientos múltiples y las llamadas a \(clustalw\) las explicaremos en la sesión4_alineamientos.
Con lo aprendido hasta ahora y los comentarios que documentan los scripts, deberías poder entender lo que hacen. Si persiste alguna duda me preguntas a mí o a \(google\).
Vuelve a explorar ambos scripts después de haber revisado la sesión4_alineamientos.
4 El lenguaje de procesamiento de patrones AWK
AWK es un lenguaje de programación diseñado para procesar de manera fácil y versátil textos con cierta estructura, ya sean ficheros o flujos de datos.
El nombre AWK deriva de las iniciales de los apellidos de sus autores: Alfred Aho, Peter Weinberger, y Brian Kernighan.
\(awk\), cuando está escrito todo en minúsculas, hace referencia al intérprete de comandos, es decir, el programa de Unix que interpreta programas escritos en el lenguaje de programación AWK. Por tanto AWK es un lenguaje interpretado por el intérprete de comando \(awk\).
AWK fue creado como un reemplazo a los algoritmos escritos en C para análisis de texto. Fue una de las primeras herramientas en aparecer en Unix (en la versión 3). Ganó popularidad rápidamente por la gran funcionalidad que permite añadir a las tuberías de comandos de Unix. Por ello se considera como una de las utilidades necesarias de este sistema operativo y de Linux.
Debido a su densa notación, todos estos lenguajes son frecuentemente usados para escribir programas de una línea, como veremos seguidamente.
4.1 Estructura de los programas AWK
En general, a \(awk\) se le dan dos piezas de datos o información: un fichero de órdenes y un archivo de entrada.
Un fichero de órdenes (que puede ser un fichero real, o puede ser incluido en la invocación de \(awk\) desde la línea de comandos) contiene una serie de sentencias o instrucciones que le indican a \(awk\) cómo procesar el fichero de entrada.
El fichero primario de entrada es normalmente texto estructurado con un formato particular, como archivos con campos separados por tabuladores (tablas). En vez de leer un archivo, \(awk\) puede procesar el flujo recibido de un comando anterior mediante el uso de pipes ‘|’.
4.1.1 Estructura básica - ejemplos genéricos
- Un programa \(AWK\) típico consiste en una serie de líneas, típicamente de la forma:
awk ‘/patrón/ { acción }’ archivo(s), o awk
‘condicione(s) { acción }’ archivo(s), o awk
‘condicione(s)’ archivo(s), o awk ‘{acción}’
archivo(s), o
donde:
- la acción por defecto es imprimir
{print} - patrón es una expresión regular
- acción es un código a ejecutar si la línea actualmente en memoria contiene el patrón y/o cumple la(s) condición(es).
La mayoría de las implementaciones de \(AWK\) usan expresiones regulares extendidas (EREs) por defecto. \(AWK\) lee línea por línea el fichero de entrada. Cuando encuentra una línea que coincide con el “patrón”, ejecuta la(s) orden(es) indicadas en “acción”.
- Para llamar a \(awk\) desde la línea de comandos, usaremos la siguiente sintaxis:
awk ‘CODIGO AWK’ ARCHIVO_A_PROCESAR
- para usarlo en una tubería de UNIX, conecta el STDOUT de un programa al STDIN de \(awk\) mediante \(|\):
STDOUT_programaX | awk ‘CODIGO AWK’ > output_file.txt
4.1.2 Formas alternativas del código AWK:
BEGIN{acción}/patrón/ && condición1 { acción }
Ejecuta las órdenes de acción al comienzo de la ejecución, antes de que
los datos comiencen a ser procesados.
/patrón/ && condiciín1 { acción1 }END { acción2
} Similar a la forma previa, pero ejecuta las órdenes de acción2
después de que todos los datos sean procesados.
/patrón/ Imprime las líneas que contienen al patrón.
{ acción } Ejecuta acción por cada línea en la
entrada.
Cada una de estas formas pueden ser incluidas varias veces en un archivo o script de AWK. El script es procesado de manera progresiva, línea por línea, de izquierda a dercha. Entonces, si hubiera dos declaraciones “BEGIN”, sus contenidos serán ejecutados en orden de aparición. Las declaraciones “BEGIN” y “END” no necesitan estar en forma ordenada.
4.1.3 Sintaxis condensada de AWK
- AWK, al ser un lenguaje de programación completo, contiene sintaxis
para escribir:
- condicionales y bucles
for$ y $wile - operadores aritméticos
+, -, *, /, %, =, ++, –, +=, -=, …) - operadores boleanos
||, && - operadores relacionales
<, <=, == !=, >=, > - funciones integradas:
length(str); int(num); index(str1, str2); split(str,arr,del); sprintf(fmt,args); substr(str,pos,len); tolower(str); toupper(str) - funciones escritas por el usuario
function FUNNAME (arg1, arg1){code} - estructuras de datos como arreglos asociativos (hashes o
diccionarios):
array[string]=value, entre otros.
- condicionales y bucles
- AWK Maneja también una serie de variables propias, de las que les resalto sólo las más usadas:
$0 guarda el valor de la fila actual en memoria de un archivo de entrada
$1, $2 ... guarda los contenidos de los campos de una fila
FILENAME nombre del archivo de entrada actualmente en procesamiento
FS separador de campos (por defecto SPACE or TAB)
NR guarda el número de campos delimitados por FS en registro o fila actual
OFS separador de campo de la salida (SPACE por defecto)
ORS separador de registro de la salida (\n por defecto)4.2 Ejemplos básicos pero muy útiles de uso de AWK
No podemos aprender aquí más que algunos idiomas de AWK muy útiles, como veremos seguidamente
Para empezar, usaremos \(AWK\) para
procesar esta simple tabla, que te pido copies al arhivo
tabla.csv haciendo uso de
nombre,salario,hrs Pepe,21,0 Lola,19,0 Carla,15.50,10 Yadira,25,20 Mary,22.50,22 Susana,17,18
Te recuerdo la secuencia de comandos para generar un archivo con cat.
cat > tabla.csv
nombre,salario,hrs
Pepe,21,0
Lola,19,0
Carla,15.50,10
Yadira,25,20
Mary,22.50,22
Susana,17,18
CTRL+dConfirma que se generó el archivo
# confirma que tienes la tabla correctamente guardada
cat tabla.csv## nombre,salario,hrs
## Pepe,21,0
## Lola,19,0
## Carla,15.50,10
## Yadira,25,20
## Mary,22.50,22
## Susana,17,18La tabla contiene los \(salarios/hr\) y \(número\ de\ horas\) trabajadas por cada trabajador/a.
Usaremos \(AWK\) para imprimir los nombres de trabajadorxs y pagos correspondientes.
\(AWK\) procesa los archivos línea a línea y separa automáticamente los campos delimitados por espacios o tabuladores. Esta es la opción por defecto.
En nuestro caso, los campos están delimitados por comas ‘,’. Podemos hacer una de dos cosas:
- Convertir el archivo csv a uno delimitado por tabuladores
- Usar las variables internas de \(AWK\) \(FS\) y \(OFS\).
Veamos ejemplos de ambos casos, empezando por el archivo tabla.tsv:
# conversión de csv a tsv usando sed
[ -s tabla.csv ] && sed 's/,/\t/g' tabla.csv > tabla.tsv
cat tabla.tsv## nombre salario hrs
## Pepe 21 0
## Lola 19 0
## Carla 15.50 10
## Yadira 25 20
## Mary 22.50 22
## Susana 17 18- Llamada de \(AWK\) sin código de filtrado o condición; sólo acción de imprimir
# Llamada de AWK sin código de filtrado o condición; sólo acción de imprimir
awk '{print}' tabla.tsv## nombre salario hrs
## Pepe 21 0
## Lola 19 0
## Carla 15.50 10
## Yadira 25 20
## Mary 22.50 22
## Susana 17 18- Llamada de \(AWK\) con código de
filtrado o condición
NR > 1para imprimir las líneas de registro (NR) mayores a 1 (elimina cabecera)
# Llamada de AWK con código de filtrado o condición;
awk 'NR > 1' tabla.tsv## Pepe 21 0
## Lola 19 0
## Carla 15.50 10
## Yadira 25 20
## Mary 22.50 22
## Susana 17 18- Llamada de \(AWK\) con código de
filtrado o condición
NR > 1para imprimir las líneas de registro (NR) mayores a 1 (elimina cabecera) y horas trabajadas > 0.
# Llamada de AWK varias condiciones de filrado;
awk 'NR > 1 && $3 > 0' tabla.tsv## Carla 15.50 10
## Yadira 25 20
## Mary 22.50 22
## Susana 17 18- Cálculo de las pagas de cada trabajador que ha trabajado > 0 hrs.
# Llamada de AWK varias condiciones de filrado y acción a ejecutar para las líneas que pasan los filtros
awk 'NR > 1 && $3 > 0 { print $1, $2* $3 }' tabla.tsv## Carla 155
## Yadira 500
## Mary 495
## Susana 306- Cálculo de las pagas de cada trabajador que ha trabajado > 0 hrs, con formato de salida tabular e impresión de cabecera en bloque BEGIN.
# Llamada de AWK varias condiciones de filrado y acción a ejecutar para las líneas que pasan los filtros
awk 'BEGIN{OFS="\t"; print "Nombre\tpaga"}NR > 1 && $3 > 0 { print $1, $2* $3 }' tabla.tsv## Nombre paga
## Carla 155
## Yadira 500
## Mary 495
## Susana 306- Lista los nombres de trabajadores que no han trabajado en el último periodo.
# Llamada de AWK varias condiciones de filrado y acción a ejecutar para las líneas que pasan los filtros
awk 'NR > 1 && $3 == 0 {print $1}' tabla.tsv## Pepe
## Lola- Imprime los nombres y números de horas trabajadas, haciendo uso del archivo tabla.csv, con formato de salida tabular, e impresión de cabecera en bloque BEGIN.
# Llamada de AWK varias condiciones de filrado y acción a ejecutar para las líneas que pasan los filtros
awk 'BEGIN{FS=","; OFS="\t"} { print $1, $3 }' tabla.csv | column -t## nombre hrs
## Pepe 0
## Lola 0
## Carla 10
## Yadira 20
## Mary 22
## Susana 18Estos simples ejemplos dan una primera idea de la flexibilidad y elegancia de \(AWK\) para procesar textos con cierto formato. Veremos más ejemplos a lo largo del curso.
\(AWK\) es mucho más poderoso de lo que muestran estos ejemplos básicos. Es un lenguaje de programación completo. Puedes ver mis apuntes del taller de biocómputo en Linux que contiene una extensa sección sobre programación en \(AWK\) para procesamiento de datos bioinformáticos.
5 Ejercicios integrativos de uso de herramientas de filtrado del Shell
5.1 Filtrado de archivos separados por tabuladores (tablas) con AWK y su graficado con R
- Asocia cada nombre de columna de la tabla assembly_summary.txt.gz con su número de campo
# asocia cada nombre de columna de la cabecera con el número de la columna correspondiente
zcat assembly_summary.txt.gz | head -2 | sed '1d; s/\t/\n/g' | cat -n## 1 # assembly_accession
## 2 bioproject
## 3 biosample
## 4 wgs_master
## 5 refseq_category
## 6 taxid
## 7 species_taxid
## 8 organism_name
## 9 infraspecific_name
## 10 isolate
## 11 version_status
## 12 assembly_level
## 13 release_type
## 14 genome_rep
## 15 seq_rel_date
## 16 asm_name
## 17 submitter
## 18 gbrs_paired_asm
## 19 paired_asm_comp
## 20 ftp_path
## 21 excluded_from_refseq
## 22 relation_to_type_material- cuenta aquellas entradas de la tabla que tienen un número de accesión revisado v2
# >>> ojo, es importante definer FS="\t", para que tome como campos sólo a aquellos separados por tabulador, y no considere los espacios <<<
# ejempo de código AWK con la estructura: BEGIN_BLOCK, condición, acción, END_BLOCK
# recuerden, como assembly_summary.txt.gz está comprimido, necesitamos zcat para poderlo leer y enviar su STDOUT a awk mediante el pipe '|'
zcat assembly_summary.txt.gz | awk 'BEGIN{FS="\t"} $16 ~ /v2$/ {n++} END{print n}' ## 4488- cuenta aquellas entradas de la tabla que tienen un número de accesión revisado v2, publicados en 2019 para genomas en estado Scaffold
# ejempo de código AWK con la estructura: BEGIN_BLOCK, condición1 && condición2 && condición3, acción, END_BLOCK
zcat assembly_summary.txt.gz | awk 'BEGIN{FS="\t"} $16 ~ /v2$/ && $15 ~ /2019/ && $12 == "Scaffold" {n++} END{print n}' ## 31- veamos las entradas de la tabla que tienen un número de accesión revisado v2, publicados en 2019 para genomas en estado Scaffold, pero imprime sólo los campos organism_name y ftp_path en formato tabla (OFS=“), imprimiendo sólo las primeras 3 líneas
# ejempo de código AWK con la estructura: BEGIN_BLOCK, condición1 && condición2 && condición3, acción
zcat assembly_summary.txt.gz | awk 'BEGIN{FS="\t"; OFS="\t"} $16 ~ /v2$/ && $15 ~ /201./ && $12 == "Scaffold" {print $8, $20}' | head -3## Leptospira interrogans ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/370/085/GCF_002370085.2_ASM237008v2
## Helicobacter pylori GAM100Ai ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/310/005/GCF_000310005.2_ASM31000v2
## Helicobacter pylori GAM101Biv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/344/945/GCF_000344945.2_ASM34494v2- veamos las entradas de la tabla que tienen un número de accesión revisado v2, publicados en 2019 para genomas en estado Scaffold, pero imprime sólo los campos organism_name y ftp_path separados por ’ ~~~ ’, imprimiendo sólo las primeras 3 líneas
# ejempo de código AWK con la estructura: BEGIN_BLOCK, condición1 && condición2 && condición3, acción
zcat assembly_summary.txt.gz | awk 'BEGIN{FS="\t"; OFS=" ~~~ "} $16 ~ /v2$/ && $15 ~ /201./ && $12 == "Scaffold" {print $8, $20}' | head -3## Leptospira interrogans ~~~ ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/370/085/GCF_002370085.2_ASM237008v2
## Helicobacter pylori GAM100Ai ~~~ ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/310/005/GCF_000310005.2_ASM31000v2
## Helicobacter pylori GAM101Biv ~~~ ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/344/945/GCF_000344945.2_ASM34494v2- genera una estadística del número de genomas por especie (columna # 8) del género Pseudomonas en formato tabular [EspecieTABnum_genomas], y muestra sólo las especies con al menos 20 genomas secuenciados. Añade una cabecera a la salida.
# genera una estadística del número de genomas por especie (columna # 8) del género Pseudomonas en formato tabular [Especie\tnum_genomas], y muestra sólo las especies con al menos 20 genomas secuenciados!
echo -e "Especie\tnum_genomas"
zcat assembly_summary.txt.gz | grep Pseudomonas | cut -f8 | sort | uniq -c | sort -nrk1 | sed 's/Pseudomonas[[:space:]]/Pseudomonas_/' | awk 'BEGIN{FS=" "; OFS="\t"}{if($1 >= 20) print $2, $1}'## Especie num_genomas
## Pseudomonas_aeruginosa 4358
## Pseudomonas_sp. 1499
## Pseudomonas_fluorescens 113
## Pseudomonas_syringae 87
## Pseudomonas_putida 86
## Pseudomonas_stutzeri 69
## Pseudomonas_syringae 64
## Pseudomonas_syringae 53
## Pseudomonas_coronafaciens 47
## Pseudomonas_protegens 37
## Pseudomonas_savastanoi 36
## Pseudomonas_savastanoi 36
## Pseudomonas_chlororaphis 24
## Pseudomonas_lundensis 21- Repitamos el ejercicio anterior, generando un archivo con campos separados por comas (csv), que se puede leer en R para genarar un \(data frame\) de \(R\) y generar una gráfica fácilmente. Para ello lo guardamos en un archivo
# Noten que primero imprimimos una cabecera al archivo Pseudomonas_species_with_gt_20_genomes.csv, y seguidamente le concatenamos con '>>' el resultado del parseo del archivo assembly_summary.txt.gz con nuestra tubería de comandos
echo -e "Especie,num_genomas" > Pseudomonas_species_with_gt_20_genomes.csv
zcat assembly_summary.txt.gz | grep Pseudomonas | cut -f8 | sort | uniq -c | sort -nrk1 | sed 's/Pseudomonas[[:space:]]/Pseudomonas_/' | awk 'BEGIN{FS=" "; OFS=","}{if($1 >= 20) print $2, $1}' >> Pseudomonas_species_with_gt_20_genomes.csv- Visualiza la cabecera del archivo que acabamos de escribir
head -5 Pseudomonas_species_with_gt_20_genomes.csv## Especie,num_genomas
## Pseudomonas_aeruginosa,4358
## Pseudomonas_sp.,1499
## Pseudomonas_fluorescens,113
## Pseudomonas_syringae,87- Llamemos a R para hacer una gráfica
$ R# 1. leemos el archivo a una estructura de datos tipo dataframe
dfr <- read.csv(file = "Pseudomonas_species_with_gt_20_genomes.csv", header = TRUE)
str(dfr)## 'data.frame': 14 obs. of 2 variables:
## $ Especie : chr "Pseudomonas_aeruginosa" "Pseudomonas_sp." "Pseudomonas_fluorescens" "Pseudomonas_syringae" ...
## $ num_genomas: int 4358 1499 113 87 86 69 64 53 47 37 ...head(dfr)## Especie num_genomas
## 1 Pseudomonas_aeruginosa 4358
## 2 Pseudomonas_sp. 1499
## 3 Pseudomonas_fluorescens 113
## 4 Pseudomonas_syringae 87
## 5 Pseudomonas_putida 86
## 6 Pseudomonas_stutzeri 69dotchart(dfr$num_genomas, labels = dfr$Especie, main = "Número de genomas por especie")
5.1.1 Reto de programación \(awk\) y \(R\)
Repite el ejercicio anterior, incluyendo el graficado del número de genomas por especie para el género Acinetobacter, pero graficando sólo aquellas especies con mínimo 5 y máximo 100 genomas
¡Felicidades, ya estás aprendiendo a programar! No es tan difícil, ¿verdad?
5.2 Ejercicios de exploración y parseo de archivos FASTA
Te propongo el siguiente ejercicio con un archivo de secuencias de DNA en formato FASTA para practicar algunos aspectos de lo aprendido en esta primera sesión.
Para correr los ejercicios, asegúrate de tener el archivo recA_Bradyrhizobium_vinuesa.fna en el directorio actual de trabajo.
El archivo recA_Bradyrhizobium_vinuesa.fna contiene secuencias del gen recA de bacterias del género Bradyrhizobium depositadas en GenBank por P. Vinuesa.
5.2.1 Búsqueda y descarga de secuencias en GenBank usando el sistema ENTREZ
Este bloque muestra el comando que usé para descargarlas usando el sistema ENTREZ de NCBI. El comando debe pegarse en la ventana superior del sistema ENTREZ.
# pega esta sentencia en la ventana de captura para interrogar la base de datos
# de nucleótidos de NCBI mediante el sistema ENTREZ
'Bradyrhizobium[orgn] AND vinuesa[auth] AND recA[gene]'
No hace falta que las descarges de NCBI. Para facilitar el acceso a las mismas, usa el siguiente código
5.2.2 Acceso a las secuencias
En primer lugar, debes estar en tu directorio de en el que deseas procesar las secuecnias y desde ahí generar descargar el archivo FASTA con las secuencias desde el repositorio GitHub usando \(wget\):
# descarga el archivo recA_Bradyrhizobium_vinuesa.fna
wget -c https://raw.githubusercontent.com/vinuesa/TIB-filoinfo/master/docs/sesion1_intro2linux/data/recA_Bradyrhizobium_vinuesa.fna
# verifica que lo tienes en tu directorio
ls recA_Bradyrhizobium_vinuesa.fna5.2.3 Inspección y estadísticas básicas de las secuencias descargadas
- ¿Cuántas secuencias hay en el archivo recA_Bradyrhizobium_vinuesa.fna?
- Explora la cabecera y cola del archivo con head y tail
- Despliega las 5 primeras lineas de cabeceras fasta usando grep y head para explorar su estructura en detalle
- Calcula el número de generos que contiene el archivo FASTA
- Calcula el número de especies que contiene el archivo FASTA
- Imprime una lista ordenada de mayor a menor, del numero de especies que contiene el archivo FASTA
5.2.4 Edición de las cabeceras FASTA mediante herramientas de filtrado de UNIX
- Explora nuevamente todas las cabeceras FASTA del archivo recA_Bradyrhizobium_vinuesa.fna usando grep y less
- Simplifica las cabeceras FASTA usando el comando sed (stream editor)
El objetivo es eliminar redundancia y los campos gb|no.de.acceso, así como todos los caracteres ‘( , ; : )’ que impedirían el despliegue de un árbol filogenético, al tratarse de caracteres reservados del formato NEWICK. Dejar solo el numero GI, así como el género, especie y cepa indicados entre corchetes.
Es decir vamos a: - reducir Bradyrhizobium a ‘B’ - eliminar ’ recombinase …’ y reemplazarlo por ‘]’ - eliminar ‘genosp.’ - sustituir espacios por guiones bajos
Nota: hagan uso de expresiones regulares como ’.*’ y ‘[[:space:]]’
- Cuando estén satisfechos con el resultado, guarden la salida del comando en un archivo llamado recA_Bradyrhizobium_vinuesa.fnaed
5.3 Solución a la práctica y un ejercicio adicional
Este ejercicio está basado en un capítulo que escribí para el manual de Sistemática Molecular y Bioinformática. Guía práctica, editado por la Facultad de Ciencias.
5.3.1 Inspección y estadísticas básicas de las secuencias descargadas
- ¿Cuántas secuencias hay en el archivo recA_Bradyrhizobium_vinuesa.fna?
grep -c '^>' recA_Bradyrhizobium_vinuesa.fna ## 125- Veamos las 5 primeras lineas de cabeceras fasta usando grep y head
grep '^>' recA_Bradyrhizobium_vinuesa.fna | head -5## >EU574327.1 Bradyrhizobium liaoningense strain ViHaR5 recombination protein A (recA) gene, partial cds
## >EU574326.1 Bradyrhizobium liaoningense strain ViHaR4 recombination protein A (recA) gene, partial cds
## >EU574325.1 Bradyrhizobium liaoningense strain ViHaR3 recombination protein A (recA) gene, partial cds
## >EU574324.1 Bradyrhizobium liaoningense strain ViHaR2 recombination protein A (recA) gene, partial cds
## >EU574323.1 Bradyrhizobium liaoningense strain ViHaR1 recombination protein A (recA) gene, partial cds- Cuenta el numero de generos y especies que contiene el archivo FASTA
grep '^>' recA_Bradyrhizobium_vinuesa.fna | cut -d' ' -f2,3 | sort | uniq -c## 18 Bradyrhizobium canariense
## 18 Bradyrhizobium elkanii
## 6 Bradyrhizobium genosp.
## 28 Bradyrhizobium japonicum
## 15 Bradyrhizobium liaoningense
## 8 Bradyrhizobium sp.
## 32 Bradyrhizobium yuanmingense- Imprime una lista ordenada de mayor a menor, del numero de especies que contiene el archivo FASTA
grep '^>' recA_Bradyrhizobium_vinuesa.fna | cut -d' ' -f2,3 | sort | uniq -c | sort -nrk1## 32 Bradyrhizobium yuanmingense
## 28 Bradyrhizobium japonicum
## 18 Bradyrhizobium elkanii
## 18 Bradyrhizobium canariense
## 15 Bradyrhizobium liaoningense
## 8 Bradyrhizobium sp.
## 6 Bradyrhizobium genosp.5.3.2 Edición de las cabeceras FASTA mediante herramientas de filtrado de UNIX
- Exploremos todas las cabeceras FASTA del archivo recA_Bradyrhizobium_vinuesa.fna usando grep
# grep '^>' recA_Bradyrhizobium_vinuesa.fna | less # para verlas por página
grep '^>' recA_Bradyrhizobium_vinuesa.fna | head # para no hacer muy extensa la salida## >EU574327.1 Bradyrhizobium liaoningense strain ViHaR5 recombination protein A (recA) gene, partial cds
## >EU574326.1 Bradyrhizobium liaoningense strain ViHaR4 recombination protein A (recA) gene, partial cds
## >EU574325.1 Bradyrhizobium liaoningense strain ViHaR3 recombination protein A (recA) gene, partial cds
## >EU574324.1 Bradyrhizobium liaoningense strain ViHaR2 recombination protein A (recA) gene, partial cds
## >EU574323.1 Bradyrhizobium liaoningense strain ViHaR1 recombination protein A (recA) gene, partial cds
## >EU574322.1 Bradyrhizobium liaoningense strain ViHaG8 recombination protein A (recA) gene, partial cds
## >EU574321.1 Bradyrhizobium liaoningense strain ViHaG7 recombination protein A (recA) gene, partial cds
## >EU574320.1 Bradyrhizobium liaoningense strain ViHaG6 recombination protein A (recA) gene, partial cds
## >EU574319.1 Bradyrhizobium yuanmingense strain ViHaG5 recombination protein A (recA) gene, partial cds
## >EU574318.1 Bradyrhizobium yuanmingense strain ViHaG4 recombination protein A (recA) gene, partial cds- simplifiquemos las cabeceras FASTA usando el comando sed (stream editor)
El objetivo es eliminar redundancia y los campos gb|no.de.acceso, así como todos los caracteres ‘( , ; : )’ que impedirían el despliegue de un árbol filogenético, al tratarse de caracteres reservados del formato NEWICK. Dejar solo el numero de accesión, así como el género, especie y cepa indicados entre corchetes.
Es decir vamos a: - reducir Bradyrhizobium a ‘B.’ - eliminar ’ recombination …’ y reemplazarlo por ‘]’ - eliminar ‘genosp.’ - sustituir espacios por guiones bajos
Noten el uso de expresiones regulares como ’.*’ y ‘[[:space:]]’
sed 's/ Bra/ [Bra/; s/|gb.*| /|/; s/Bradyrhizobium /B/; s/genosp\. //; s/ recomb.*/]/; s/[[:space:]]/_/g;' recA_Bradyrhizobium_vinuesa.fna | grep '>' | head -5
sed 's/ Bra/ [Bra/; s/|gb.*| /|/; s/Bradyrhizobium /B/; s/genosp\. //; s/ recomb.*/]/; s/[[:space:]]/_/g;' recA_Bradyrhizobium_vinuesa.fna | grep '>' | tail -5## >EU574327.1_[Bliaoningense_strain_ViHaR5]
## >EU574326.1_[Bliaoningense_strain_ViHaR4]
## >EU574325.1_[Bliaoningense_strain_ViHaR3]
## >EU574324.1_[Bliaoningense_strain_ViHaR2]
## >EU574323.1_[Bliaoningense_strain_ViHaR1]
## >AY591544.1_[Bjaponicum_bv._genistearum_strain_BC-P14]
## >AY591543.1_[Bbeta_strain_BC-P6]
## >AY591542.1_[Bcanariense_bv._genistearum_strain_BC-P5]
## >AY591541.1_[Bcanariense_bv._genistearum_strain_BC-C2]
## >AY591540.1_[Balpha_bv._genistearum_strain_BC-C1]- Cuando estamos satisfechos con el resultado, guardamos la salida del comando en un archivo usando ‘>’ para redirigir el flujo de STDOUT a un archivo de texto
sed 's/ Bra/ [Bra/; s/|gb.*| /|/; s/Bradyrhizobium /B/; s/genosp\. //; s/ recomb.*/]/; s/[[:space:]]/_/g;' recA_Bradyrhizobium_vinuesa.fna > recA_Bradyrhizobium_vinuesa.fnaed
5.3.3 Generación automática de archivos FASTA especie-específicos (avanzado)
- Convertir archivos FASTA a formato “FASTAB” usando perl 1-liners.
Vamos a transformar los FASTAS de tal manera que las secuencias queden en la misma línea que su cabecera, separada de ésta por un tabulador. Esto puede ser muy útil para filtrar el archivo resultante con grep. Veamos un ejemplo:
perl -pe 'unless(/^>/){s/\n//g}; if(/>/){s/\n/\t/g}; s/>/\n>/' recA_Bradyrhizobium_vinuesa.fnaed | head##
## >EU574327.1_[Bliaoningense_strain_ViHaR5] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGACTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCAGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574326.1_[Bliaoningense_strain_ViHaR4] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574325.1_[Bliaoningense_strain_ViHaR3] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574324.1_[Bliaoningense_strain_ViHaR2] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574323.1_[Bliaoningense_strain_ViHaR1] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574322.1_[Bliaoningense_strain_ViHaG8] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574321.1_[Bliaoningense_strain_ViHaG7] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574320.1_[Bliaoningense_strain_ViHaG6] ATGAAGCTCGGCAAGAACGACCGGTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGACATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGGCGTATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACGCTGGTGCGCTCGGGCGCTGTCGATGTGCTGGTGATCGATTCGGTTGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574319.1_[Byuanmingense_strain_ViHaG5] ATGAAGCTCGGCAAGAACGACCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCTTGCCCAAGGGGCGCATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCCGTCGACGTTCTCGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGACTTCAAGCACGGTTGATGAGCCAGGCGCTGCGCAAGCTCACGGCCTCCATCAACAAATCCAACACGATGGTGATCTTCATCAACCAGATCperl -pe 'unless(/^>/){s/\n//g}; if(/>/){s/\n/\t/g}; s/>/\n>/' recA_Bradyrhizobium_vinuesa.fnaed > recA_Bradyrhizobium_vinuesa.fnaedtab- Filtrar el archivo fnaedtab generado en 9 para obtener solo las secuencias de B._yuanmingense del mismo, guardarlo en un archivo y convertirlo de nuevo a formato FASTA.
grep yuanmingense recA_Bradyrhizobium_vinuesa.fnaedtab | head -5## >EU574319.1_[Byuanmingense_strain_ViHaG5] ATGAAGCTCGGCAAGAACGACCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCTTGCCCAAGGGGCGCATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCCGTCGACGTTCTCGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGACTTCAAGCACGGTTGATGAGCCAGGCGCTGCGCAAGCTCACGGCCTCCATCAACAAATCCAACACGATGGTGATCTTCATCAACCAGATC
## >EU574318.1_[Byuanmingense_strain_ViHaG4] ATGAAGCTCGGCAAGAACGACCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCTTGCCCAAGGGGCGCATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCCGTCGACGTTCTCGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGACTTCAAGCACGGTTGATGAGCCAGGCGCTGCGCAAGCTCACGGCCTCCATCAACAAATCCAACACGATGGTGATCTTCATCAACCAGATC
## >EU574297.1_[Byuanmingense_strain_InRo02] ATGAAGCTCGGCAAGAACGACCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTGGGCATCGGCGGCCTGCCCAAGGGACGCATCGTCGAGATCTACGGGCCGGAATCTTCGGGCAAGACCACGCTCGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGTGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCCGTCGATGTTCTGGTGGTCGATTCGGTGGCGGCGCTGGTGCCGAAGGCCGAGCTCGAAGGCGAGATGGGTGACGCGCTGCCTGGCCTGCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCATCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574296.1_[Byuanmingense_strain_InKo02] ATGAAGCTCGGCAAGAACGATCGCTCCATGGACATCGAGGCGGTCTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGACGCATCGTCGAGATCTACGGGCCGGAATCCTCCGGCAAGACCACGCTCGCGCTGCATACGGTGGCGGAGGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGTTCGGGCGCCGTCGACGTTCTCGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGCCTGCAGGCGCGGTTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATC
## >EU574295.1_[Byuanmingense_strain_InKo01] ATGAAGCTCGGCAAGAACGATCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCCTGCCCAAGGGACGCATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCCCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCGGTCGATGTTCTCGTGGTCGATTCGGTGGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGCTTGCAGGCCCGTCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCAACACCATGGTGATCTTCATCAACCAGATCgrep yuanmingense recA_Bradyrhizobium_vinuesa.fnaedtab > recA_Byuanmingense.fnaedtab- Estas dos lineas no contienen nada nuevo en cuanto a sintaxis. Simplemente llamamos a perl para sustituir los tabuladores por saltos de linea y asi reconstituir el FASTA.
perl -pe 'if(/^>/){s/\t/\n/}' recA_Byuanmingense.fnaedtab | head -5## >EU574319.1_[Byuanmingense_strain_ViHaG5]
## ATGAAGCTCGGCAAGAACGACCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCTTGCCCAAGGGGCGCATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCCGTCGACGTTCTCGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGACTTCAAGCACGGTTGATGAGCCAGGCGCTGCGCAAGCTCACGGCCTCCATCAACAAATCCAACACGATGGTGATCTTCATCAACCAGATC
## >EU574318.1_[Byuanmingense_strain_ViHaG4]
## ATGAAGCTCGGCAAGAACGACCGCTCCATGGACATCGAGGCGGTGTCCTCCGGCTCGCTCGGGCTCGATATCGCGCTCGGCATCGGCGGCTTGCCCAAGGGGCGCATCGTCGAGATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCATACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGCGCCTTCATCGACGCCGAGCACGCGCTCGATCCGGTCTATGCCCGCAAGCTCGGCGTCAACATCGACGAGCTCCTGATCTCGCAGCCCGACACCGGCGAGCAGGCGCTGGAGATCTGCGACACCCTGGTGCGCTCGGGCGCCGTCGACGTTCTCGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAGCTCGAGGGCGAGATGGGCGATGCGCTGCCGGGACTTCAAGCACGGTTGATGAGCCAGGCGCTGCGCAAGCTCACGGCCTCCATCAACAAATCCAACACGATGGTGATCTTCATCAACCAGATC
## >EU574297.1_[Byuanmingense_strain_InRo02]perl -pe 'if(/^>/){s/\t/\n/}' recA_Byuanmingense.fnaedtab > recA_Byuanmingense.fna- Llamar a un bucle for de shell para generar archivos fastab para todas las especies
for sp in $(cut -d_ -f2 recA_Bradyrhizobium_vinuesa.fnaedtab | sort -u | sed 's/\[//'); do
grep "$sp" recA_Bradyrhizobium_vinuesa.fnaedtab > "recA_${sp}.fnaedtab"
done- Veamos el resultado
ls *fnaedtab## recA_Balpha.fnaedtab
## recA_Bbeta.fnaedtab
## recA_Bcanariense.fnaedtab
## recA_Belkanii.fnaedtab
## recA_Bjaponicum.fnaedtab
## recA_Bliaoningense.fnaedtab
## recA_Bradyrhizobium_vinuesa.fnaedtab
## recA_Bsp..fnaedtab
## recA_Byuanmingense.fnaedtabhead -5 recA_Bjaponicum.fnaedtab## >EU574316.1_[Bjaponicum_strain_NeRa16] ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCGGGTTCTCTCGGGCTCGACATTGCACTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCACGCAAGCTCGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACCGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC
## >EU574315.1_[Bjaponicum_strain_NeRa15] ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCCGGTTCTCTCGGGCTCGACATTGCACTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCACGCAAGCTCGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACCGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC
## >EU574314.1_[Bjaponicum_strain_NeRa14] ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCCGGTTCTCTCGGGCTCGACATTGCGCTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGAATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCGCGCAAGCTGGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACGGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGYCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC
## >EU574313.1_[Bjaponicum_strain_NeRa12] ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCCGGTTCTCTCGGGCTCGACATTGCACTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCACGCAAGCTCGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACCGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC
## >EU574312.1_[Bjaponicum_strain_NeRa11] ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCCGGTTCTCTCGGGCTCGACATTGCACTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCACGCAAGCTCGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACCGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC- Finalmente convertimos todos los archivos fnatabed a FASTA con el siguiente bucle for:
for file in $(ls *fnaedtab | grep -v vinuesa); do perl -pe 'if(/^>/){s/\t/\n/}' $file > ${file%.*}.fas; done- Visualizemos las cabeceras de dos archivos FASTA especie-específicos
grep '>' recA_Bjaponicum.fas | head -5## >EU574316.1_[Bjaponicum_strain_NeRa16]
## >EU574315.1_[Bjaponicum_strain_NeRa15]
## >EU574314.1_[Bjaponicum_strain_NeRa14]
## >EU574313.1_[Bjaponicum_strain_NeRa12]
## >EU574312.1_[Bjaponicum_strain_NeRa11]- y confirmemos que son fastas regulares
head -6 recA_Bjaponicum.fas## >EU574316.1_[Bjaponicum_strain_NeRa16]
## ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCGGGTTCTCTCGGGCTCGACATTGCACTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCACGCAAGCTCGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACCGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC
## >EU574315.1_[Bjaponicum_strain_NeRa15]
## ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCCGGTTCTCTCGGGCTCGACATTGCACTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGCATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCACGCAAGCTCGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACCGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGCCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC
## >EU574314.1_[Bjaponicum_strain_NeRa14]
## ATGAAGCTCGGCAAGAACGACCGGTCGATGGATGTCGAGGCGGTGTCCTCCGGTTCTCTCGGGCTCGACATTGCGCTGGGGATCGGCGGTCTGCCCAAGGGGCGCGTCGTGGAAATCTACGGGCCGGAATCCTCGGGCAAGACCACGCTGGCGCTGCACACGGTGGCGGAAGCGCAGAAGAAGGGCGGAATCTGTGCCTTCATCGACGCCGAGCACGCGCTCGACCCGGTCTATGCGCGCAAGCTGGGGGTCAATATCGACGAGCTGCTGATTTCGCAGCCGGACACGGGCGAGCAGGCGCTGGAAATCTGCGACACGCTGGTGCGCTCGGGTGCGGTGGACGTTCTGGTGGTCGATTCGGTCGCGGCGCTGGTGCCGAAGGCCGAACTCGAGGGCGAGATGGGCGATGCGCTGCCGGGTCTCCAGGCCCGYCTGATGAGCCAGGCGCTGCGCAAGCTGACGGCCTCCATCAACAAGTCCCACACCATGGTGATCTTCATCAACCAGATC- si quieren, borren todos los archivos generados, para empezar con un directorio de trabajo limpio, para repetir el ejercicio ;)
rm *fnaed *fnaedtab *fas Stenotrophomonas_complete_genomes_and_ftp_paths.txt empty_file Pseudomonas_species*6 Reto de programación - ejercicio de parseo de archivos FASTA
Como ejercicio, para repasar lo que hemos aprendido en esta sesión les propongo repetir el ejercicio de parseo de archivos FASTA pero con secuencias del gen rpoB de Bradyrhizobium
Los que no tengan instalado MobaXterm, tendrán el reto adicional de instalarlo y familiarizarse con él.
6.1 Inspección y estadísticas básicas de las secuencias descargadas
- Descargar las secuencias de NCBI usando el portal ENTREZ nucleotides con: ‘Bradyrhizobium[orgn] AND vinuesa[auth] AND rpoB[gene]’
- Renombra el archivo descargado a rpoB_Bradyrhizobium_vinuesa.fna
- ¿Cuántas secuencias hay en el archivo rpoB_Bradyrhizobium_vinuesa.fna?
- Explora la cabecera y cola del archivo con head y tail
- Despliega las 5 primeras lineas de cabeceras fasta usando grep y head para explorar su estructura en detalle
- Calcula el número de generos que contiene el archivo FASTA
- Calcula el número de especies que contiene el archivo FASTA
- Imprime una lista ordenada de mayor a menor, del numero de especies que contiene el archivo FASTA
6.2 Edición de las cabeceras FASTA mediante herramientas de filtrado de UNIX
- Explora nuevamente todas las cabeceras FASTA del archivo rpoB_Bradyrhizobium_vinuesa.fna usando grep y less
- Simplifica las cabeceras FASTA usando el comando sed (stream editor)
El objetivo es eliminar redundancia y longitud excesiva de las cabeceras FASTA, así como todos los caracteres ‘( , ; : )’ que impedirían el despliegue de un árbol filogenético, al tratarse de caracteres reservados del formato NEWICK. Dejar solo el numero de accesión, así como el género, especie y cepa indicados entre corchetes.
- Cuando estén satisfechos con el resultado, guarden la salida del comando en un archivo llamado rpoB_Bradyrhizobium_vinuesa.fnaed
7 Programando un pipeline en \(Bash\)
Para finalizar este tutorial, suguiero que vean el el \(script\) run_phylip.sh que puedes descargar desde aqui. Es mucho más largo y complejo que los anteriores, e integra todos de los elementos sintácticos mostrados en secciones anteriores, y otros nuevos, explicados en este extenso y detallado tutorial avanzado de bioćomputo en sistemas Linux.
Cuando estamos aprendiendo a programar, es muy importante leer mucho código bien escrito. No se aprende a programar bien sólo conociendo la sintaxis y generalidades, sino viendo cómo se pueden escribir programas completos de manera limpia, fácil de leer, mantener y modificar. Estos atributos, junto con el buen juicio y principios de programación, incluyendo la modularización en funciones, documentación adecuada de cada componente, incluyendo la interfaz de usuario, son los que debemos desarrollar. El estudio detallado e imitación de buenos programas es lo que más rápidamente nos lleva a escribir mejores programas.
El \(script\) run_phylip.sh, además de útil para hacer análisis filogenéticos basados en matrices de distancias, te presenta código que implementa buenas prácticas de programación \(Bash\).
8 Consideraciones finales y referencias recomendadas para continuar aprendiendo programación \(Shell\)
Llegamos al final de este tutorial, ¡¡¡enhorabuena!!! Has aprendido mucho, pero lo bueno es que todavía queda un largo camino por andar.
\(Bash\) sin duda es muy útil si vas a usarlo como lenguaje pegamento para conectar múltiples utilerías o programas como demuestra el \(script\) run_phylp.sh que presentamos en la sección anterior, pero no es el apropiado para escribir programas complejos que tienen que hacer procesamiento numérico, gráfico, estadístico de datos, interactuar con bases de datos, parsear archivos de estructura compleja, etc. Para ello deberás aprender otros lenguajes más poderosos, capaces de manejar estructuras de datos complejas (multidimensionales como hashes de hashes, arreglos de hashes, etc.) y que disponen de repositorios gigantescos de librerías o paquetes especializados de código para todo lo que puedas imaginar. Pero aprender primero \(Shell\) y \(AWK\) es sin duda lo mejor que puedes hacer para iniciarte en la programación en sistemas UNIX o GNU/Linux. Te facilitará el camino posterior y te permitirá dominar el sistema operativo, ya que el \(Shell\) es su interfaz programática.
Sabiendo \(Bash\) y \(AWK\) te será muy fácil entender \(Perl\), que con sus 42,000 paquetes y 197,000 módulos en el repositorio CPAN, sería una excelente elección como siguiente lenguaje a aprender. Sin duda Python y R hay que añadirlos también a la lista de lenguajes interpretados a aprender.
Recomiendo las siguientes guías para apoyarte en tu proceso de aprendizaje del \(Shell\). Disfruta el camino, saludos!