Programación en shell

Scripts de shell

Los scripts (archivos con un conjunto de órdenes) permiten automatizar y acelerar el trabajo.

Vamos a crear nuestro primer script. Para ello en un editor de texto escribiremos lo siguiente y lo guardaremos con el nombre hola.sh.

#!/bin/bash
# Este es nuestro primer programa
echo Hola Mundo

A continuación iremos a la terminal y lo ejecutaremos:

~$ ./hola.sh

La primera línea (la shebang, #!/bin/bash) es la que le dice al sistema con qué programa ejecutarlo. La segunda es un comentario para consumo humano ya que va precedido de #. La tercera línea es la que indica qué queremos imprimir en pantalla (echo)

Cuestión

¿Podrías crear un script que tomara un archivo FASTA y contara el número de secuencias?

Queremos generar un script que funcione así:

contar_seqs.sh sequences.fasta

Pista: para que un script tome el argumento que pasamos luego (sequences.fasta en este caso) se hace llamando $1

Nota

#!/bin/bash
grep "^>" "$1" | wc -l

Las comillas dobles alrededor de "$1" no son decoración: si el argumento tuviera espacios (como los nombres de archivo mal escritos que viste en el módulo 2, del tipo "Muestra 01.fastq"), sin comillas el script lo trocearía en varios argumentos. Poner siempre "$1" entre comillas es un hábito que te evitará errores difíciles de depurar.

Variables

Podemos guardar valores en variables. Estas se asignan con nombre=valor (sin espacios alrededor del =) y se recuperan anteponiendo $:

#!/bin/bash
to_print='Hola caracola'
echo to_print ahora es $to_print
to_print=5.5
echo to_print ahora es $to_print

Los nombres de variable solo pueden contener letras, números y _, y por convención se escriben en minúsculas — así se evita chocar con variables del propio sistema, que suelen ir en mayúsculas (PATH, HOME…).

Bucles `for`

for repite un bloque de comandos una vez por cada elemento de una lista. Esa lista puede venir de tres formas:

Un rango numérico: for variable in {1..10}
Una lista de valores escritos a mano: for variable in SiteA SiteB SiteC
La salida de un comando: for variable in $(comando)

Cuestión

¿Qué valores iría tomando var1 en los siguientes casos?

for var1 in $(ls *.txt)
for var1 in $(grep '>' sequences.fasta)
for var1 in $(grep '>' sequences.fasta | sed 's/>//')
for var1 in $(grep '>' sequences.fasta | sed 's/>//' | cut -f 1 -d'|')

Tras un for los comandos se ejecutan dentro de un bloque do-done:

#!/bin/bash
for numero in {1..20};
do
  echo Este es el número: $numero
done

Cuestión

¿Podrías contar cuántas observaciones hay de cada site en species_observations.csv sin escribir los nombres de sitio a mano?

Nota

#!/bin/bash
for site in $(cut -d',' -f3 species_observations.csv | tail -n +2 | sort -u);
do
  n=$(grep "$site" species_observations.csv | wc -l)
  echo "$site: $n observaciones"
done

Renombrar archivos en bloque

Un uso muy habitual de for es renombrar muchos archivos a la vez — por ejemplo, los que descargas de un secuenciador o de una base de datos suelen traer nombres largos con información redundante, de los que solo te interesa una parte. Combinando for con cut (que ya conoces del módulo 3) y mv puedes quedarte con la parte significativa del nombre sin tocarlos uno a uno.

Imagina estos tres archivos, con el run del secuenciador, la lane en la que se corrió, el nombre de la muestra, la especie y la dirección de las lecturas (R1 o R2) todo junto en el nombre:

RUN2024_LANE3_S001_Quercus_robur_R1.fastq.gz
RUN2024_LANE3_S002_Pinus_sylvestris_R2.fastq.gz
RUN2024_LANE3_S003_Fragaria_vesca_R1.fastq.gz

Y quieres renombrarlos solo con el nombre de la muestra y la dirección de la lectura (S001_R1.fastq.gz, S002_R2.fastq.gz, …)

Cuestión

¿Podrías con las herramientas que conoces extraer del nombre del archivo la información que te interesa?
¿Podrías guardarlo en una variable en un script?
¿Podrías incluirlo en un for? ¿Sobre qué hacemos el for?

Nota

#!/bin/bash

for file in $(ls *.fastq.gz)
do
   sample=$(echo "$file" | cut -d'_' -f 3)
   direct=$(echo "$file" | cut -d'_' -f 6 | cut -d'.' -f1)
   
   new_file="${sample}_${direct}.fastq.gz"
   mv "$file" "$new_file"
done

Siempre, siempre, SIEMPRE, antes de lanzar un for (especialmente si hace uso de rm) asegúrate con echo de que los nombres son los que esperas.

`awk`: procesar texto por columnas

awk es un lenguaje pensado para archivos con estructura de columnas: lee el archivo línea a línea y, en cada línea, separa automáticamente los campos — por defecto por espacios o tabuladores — a los que te refieres como $1, $2, $3… ($0 es la línea completa, sin dividir). Su forma básica es:

awk '{acción}' archivo

-F indica el separador de campo, igual que -d en cut.

Imprimir columnas

print muestra lo que le pases, separado por comas en el propio print (que awk traduce a un espacio en la salida). Vamos a usar blast_tabular.tsv, un resultado de BLAST en formato tabular (-outfmt 6, separado por tabuladores) donde cada fila es el mejor alineamiento de un gen contra una base de datos de referencia:

~$ awk -F'\t' '{print $1, $2}' blast_tabular.tsv
pka ref_pka_Athaliana
fos ref_fos_Athaliana
cox ref_cox_Athaliana
lac ref_lac_Athaliana
p53 ref_p53_Athaliana
rna ref_rna_Athaliana

$1 es el identificador del gen consultado (query) y $2 el de la secuencia de referencia con la que ha alineado — las dos primeras columnas de cualquier salida de BLAST tabular. A diferencia de cut -f1,2, que solo puede extraer columnas tal cual, awk te deja además calcular con ellas.

Calcular una división: porcentaje de cobertura

Cada fila de blast_tabular.tsv trae, entre otras, la columna 4 (length, la longitud del alineamiento) y la columna 13 (qlen, la longitud total del gen consultado). La cobertura del gen es la proporción de su longitud total que ha quedado cubierta por el alineamiento: length / qlen. awk hace aritmética con + - * / directamente sobre los campos:

~$ awk -F'\t' '{printf "%s\t%.1f%%\n", $1, ($4/$13)*100}' blast_tabular.tsv
pka 90.0%
fos 93.8%
cox 100.0%
lac 45.0%
p53 100.0%
rna 100.0%

printf (en vez de print) permite controlar el formato: %s inserta el nombre del gen como texto, %.1f%% inserta el número con un decimal seguido de un % literal (por eso va doblado, %%), y \t/\n separan con un tabulador y saltan de línea al final, igual que en los echo/sed que ya conoces.

Nota

Fíjate en lac: solo un 45% de cobertura. Un pident (columna 3, porcentaje de identidad) alto no basta para confiar en un hit de BLAST — un alineamiento muy parecido pero muy corto respecto al gen completo puede ser un dominio compartido o una coincidencia parcial, no el gen entero.

Ejercicios

Escribe un script contar_observaciones.sh que reciba un CSV como argumento ($1) y muestre cuántas filas de datos tiene, sin contar las líneas de cabecera. Pruébalo con species_observations.csv (que tiene 3 líneas de cabecera antes de los datos).
Usando la forma de for con una lista de valores escritos a mano (for site in SiteA SiteB SiteC SiteD SiteE), escribe un script que cree un directorio por cada sitio y guarde dentro un archivo observaciones.csv solo con las filas de species_observations.csv de ese sitio.
Usando gene_expression.csv, escribe un awk que muestre, para cada gen, su identificador y el cambio entre sample_A y sample_B (sample_A/sample_B), con dos decimales.
Tienes fotos de cámara trampa con nombres como IMG_2024-05-12_CameraA_0001.jpg. Escribe un script que recorra con for todos los archivos IMG_*.jpg de un directorio y los renombre conservando solo la cámara y el número (CameraA_0001.jpg).
(Reto) Convierte el awk de la cobertura en un script cobertura.sh que reciba el archivo BLAST tabular como argumento ($1), de modo que puedas ejecutarlo como ./cobertura.sh blast_tabular.tsv.

Soluciones

#!/bin/bash
total=$(tail -n +4 "$1" | wc -l)
echo "$1 tiene $total observaciones"

~$ chmod +x contar_observaciones.sh
~$ ./contar_observaciones.sh species_observations.csv
species_observations.csv tiene 117 observaciones

#!/bin/bash
for site in SiteA SiteB SiteC SiteD SiteE
do
  mkdir -p "$site"
  grep "$site" species_observations.csv > "$site/observaciones.csv"
done

~$ tail -n +4 gene_expression.csv | awk '{printf "%s\t%.2f\n", $1, $2/$3}'
P12345   0.82
P67890   0.00
Q12345   0.90
P23456   0.92
P34567   1.11
P45678   1.14

tail -n +4 se salta la cabecera antes de pasar el archivo a awk; como gene_expression.csv ya usa espacios como separador, no hace falta -F.

#!/bin/bash
for file in IMG_*.jpg
do
  nuevo=$(echo "$file" | cut -d'_' -f3,4 --output-delimiter='_')
  mv "$file" "$nuevo"
done

#!/bin/bash
awk -F'\t' '{printf "%s\t%.1f%%\n", $1, ($4/$13)*100}' "$1"

~$ chmod +x cobertura.sh
~$ ./cobertura.sh blast_tabular.tsv
pka  90.0%
fos  93.8%
cox  100.0%
lac  45.0%
p53  100.0%
rna  100.0%

Scripts de shell

Variables

Bucles for

Renombrar archivos en bloque

awk: procesar texto por columnas

Imprimir columnas

Calcular una división: porcentaje de cobertura

Ejercicios

Soluciones

Bucles `for`

`awk`: procesar texto por columnas