15  Ejercicios en equipo

Reúnanse en equipos de 3 personas para completar uno de los ejercicios propuestos en esta sección, algunos de los ejercicios provienen de la actividad de 6. Extracción de información de archivos fastq del curso de GNU/Linux para Bioinformática organizado por RSG Ecuador.

15.1 Extracción de información de archivos fastq

En este ejercicio se obtendrá información de archivos fastq, que cotienen secuencias de nucleótidos y la calidad de su proceso de secuenciación. Para esto usaremos los archivos secuencias1.fastq, secuencias2.fastq y secuencias3.fastq que se encuentra en la carpeta _files dentro del directorio raíz del repositorio de GitHub unix.bioinfo.rsgecuador .

Ejercicio 10.1 Contar el número de secuencias con diez “N” seguidas por cada archivo FASTQ

La letra N representa un nucleótido que no pudo ser leído correctamente, y se reporta como ninguno, o missing data. Imprime el número de secuencias que tenga diez N seguidas por cada uno de los archivos fastq de la carpeta _files.

Note

Es posible realizar esto con un comando de una sola línea, intenta resolverlo de esta forma.

#!/bin/bash
# File: countseq.sh

# Ruta a la carpeta que contiene los archivos .fastq
FOLDER="_files"

# Iterar sobre cada archivo .fastq en la carpeta
for file in "$FOLDER"/*[1-3].fastq; do
  # Imprimir el nombre del archivo
  echo "Archivo: $file"
  
  # Usar grep para buscar secuencias que contengan 10 'N' seguidas y contar cuántas veces aparecen
  count=$(grep -c 'NNNNNNNNNN' $file )
  
  # Imprimir el número de secuencias encontradas
  echo "Número de secuencias con 10 N seguidas: $count"
  echo "---------------------------------------"
done

Recuerda para ejecutar este codigo deberas encontrarte dentro de la carpeta descomprimida unix.bioinfo.rsgecuador-gh-pages/y ejecutar: ./countseq.sh

Salida:

  • secuencias1.fastq - 1
  • secuencias2.fastq - 58
  • secuencias3.fastq - 75

Ejercicio 10.2 Contar el número de secuencias con un string en cada archivo FASTQ

Ahora, determina el número de secuencias (lineas) de cada archivo .fastq por separado. Se conoce que los títulos de las corridas de las secuencias en cada archivo .fastq empiezan con el string @SRR098026.

#!/bin/bash
# File: countNseq.sh

# Ruta a la carpeta que contiene los archivos .fastq
FOLDER="_files"

# Iterar sobre cada archivo .fastq en la carpeta
for file in "$FOLDER"/*[1-3].fastq; do
  echo "Archivo: $file"
  
  # Usar grep para contar cuántas líneas empiezan con '@SRR098026'
  count=$(grep -c '^@SRR098026' $file )
  
  echo "Número de secuencias: $count"
  echo "---------------------------------------"
done

Salida:

  • secuencias1.fastq - 99
  • secuencias2.fastq - 75
  • secuencias3.fastq - 75

Ejercicio 10.3 Identifica el número de secuencias con calidad baja

En los archivos FASTQ, la tercera línea de cada entrada contiene el símbolo +, mientras que la cuarta línea representa la calidad de la secuencia en formato ASCII. Determina cuántas secuencias tienen una calidad promedio por debajo de un umbral específico (por ejemplo, una calidad promedio menor a 20).

Pista

Usa herramientas como grep y un script auxiliar para convertir los símbolos ASCII a valores numéricos y calcular su promedio.

Recuerda que los archivos se encuentran en unix.bioinfo.rsgecuador-gh-pages/_files/

# Umbral de calidad
umbral=20

# Extraer las líneas de calidad (cuarta línea de cada entrada)
grep -A 3 "^@" *.fastq | grep -v "^@" | grep -v "+" > quality_lines.txt

# Convertir las líneas de calidad ASCII a números, calcular el promedio, y contar las que sean menores al umbral
cat quality_lines.txt | while read line; do
  # Convertir caracteres ASCII a valores numéricos y calcular la calidad promedio
  sum=0
  count=0
  for i in $(echo "$line" | od -An -t u1); do
    sum=$((sum + i - 33))  # Restar 33 para convertir de ASCII a valor de calidad
    count=$((count + 1))
  done

  # Calcular la calidad promedio
  if [ "$count" -gt 0 ]; then
    calidad=$((sum / count))  # Calcular el promedio sin necesidad de bc
  else
    calidad=0
  fi

  # Comprobar si la calidad es menor que el umbral y almacenar
  if [ "$calidad" -lt "$umbral" ]; then
    echo "$calidad" >> low_quality_sequences.txt
  fi
done

# Contar secuencias con calidad promedio menor al umbral
echo "Número de secuencias con calidad promedio < $umbral: $(wc -l < low_quality_sequences.txt)"

Número de secuencias con calidad promedio < 20: 763

Ejercicio 10.4 Cuenta cuántas secuencias únicas hay en cada archivo

Determina cuántas secuencias únicas (sin repeticiones) existen en cada archivo. Asegúrate de ignorar las líneas que comienzan con @ o +, ya que no forman parte de las secuencias.

Pista

Usa grep para extraer solo las líneas correspondientes a las secuencias y después sort y uniq para contarlas.

Recuerda que los archivos se encuentran en unix.bioinfo.rsgecuador-gh-pages/_files/

# Extraer solo las secuencias (segunda línea de cada entrada)
grep -A 3 "^@" *.fastq | grep -v "^@" | grep -v "+" | grep -v "^--$" > sequences_only.txt

# Contar las secuencias únicas
echo "Número de secuencias únicas: $(sort sequences_only.txt | uniq | wc -l)"

Número de secuencias únicas: 620

Ejercicio 10.5 Encuentra la longitud promedio de las secuencias

Calcula la longitud promedio de las secuencias en cada archivo. Esto requiere sumar las longitudes de todas las secuencias y dividirlas entre el número total de secuencias.

Pista

Extrae las líneas correspondientes a las secuencias y utiliza wc -c para sumar los caracteres, luego divide por el número de secuencias.

# calcular_longitud_promedio.sh
# Extraer solo las secuencias (segunda línea de cada entrada)
grep -A 3 "^@" *.fastq | grep -v "^@" | grep -v "+" | grep -v "^--$" > sequences_only.txt

# Calcular la longitud total de las secuencias
total_length=0
while read -r sequence; do
  total_length=$((total_length + ${#sequence}))
done < sequences_only.txt

# Contar el número total de secuencias
num_sequences=$(wc -l < sequences_only.txt)

# Calcular la longitud promedio
if [ "$num_sequences" -gt 0 ]; then
  average_length=$((total_length / num_sequences))
else
  average_length=0
fi

echo "Longitud promedio de las secuencias: $average_length"

Longitud promedio de las secuencias: 58

Ejercicio 10.6 Calcular el porcentaje de GC con una lista de archivos

Tomando en cuenta el ejercicio anterior, calcula el porcentaje de GC de los tres archivos fastq secuencias1.fastq, secuencias2.fastq, secuencias3.fastq utilizando ciclos y/o bucles.

#!/bin/bash
# GC_per.sh

# Directorio donde están los archivos FASTQ
directorio="$HOME/unix.bioinfo.rsgecuador-gh-pages/_files"

# Loop para procesar cada archivo FASTQ
for archivo in "$directorio"/*.fastq; do
  if [ -f "$archivo" ]; then
    echo "Procesando archivo: $archivo"

    # Inicializar contadores
    total_bases=0
    gc_bases=0

    # Leer el archivo línea por línea
    contador=0
    while IFS= read -r linea; do
      ((contador++))
      
      # Solo procesamos las líneas de secuencias (segunda línea de cada bloque)
      if ((contador % 4 == 2)); then
        # Contamos las bases G y C
        gc_bases=$((gc_bases + $(echo "$linea" | grep -o 'G' | wc -l) + $(echo "$linea" | grep -o 'C' | wc -l)))
        total_bases=$((total_bases + ${#linea}))
      fi
    done < "$archivo"

    # Calculamos el porcentaje de GC
    if [ "$total_bases" -gt 0 ]; then
      porcentaje_gc=$((gc_bases * 100 / total_bases))
      echo "Porcentaje de GC en $archivo: $porcentaje_gc%"
    else
      echo "No se encontraron secuencias válidas en $archivo."
    fi
  else
    echo "El archivo $archivo no existe."
  fi
done

Procesando archivo: secuencias1.fastq: 27% Procesando archivo: secuencias2.fastq: 11% Procesando archivo: secuencias3.fastq: 2% Procesando archivo: secuencias_bash.fastq: 2%

Ejercicio 10.7 Calcular el contenido de bases por secuencias

Calcula el contenido de cada una de las bases en las secuencias del archivo secuencias1.fastq, A,T,G,C y N.

#!/bin/bash

# Archivo de entrada
archivo="secuencias1.fastq"

# Inicializar contadores
conteo_A=0
conteo_T=0
conteo_G=0
conteo_C=0
conteo_N=0

# Leer el archivo línea por línea
contador=0
while IFS= read -r linea; do
  ((contador++))

  # Procesar solo líneas de secuencias (segunda línea de cada bloque de 4)
  if ((contador % 4 == 2)); then
    conteo_A=$((conteo_A + $(echo -n "$linea" | grep -o 'A' | wc -l)))
    conteo_T=$((conteo_T + $(echo -n "$linea" | grep -o 'T' | wc -l)))
    conteo_G=$((conteo_G + $(echo -n "$linea" | grep -o 'G' | wc -l)))
    conteo_C=$((conteo_C + $(echo -n "$linea" | grep -o 'C' | wc -l)))
    conteo_N=$((conteo_N + $(echo -n "$linea" | grep -o 'N' | wc -l)))
  fi
done < "$archivo"

# Mostrar resultados
echo "Contenido de bases en $archivo:"
echo "A: $conteo_A"
echo "T: $conteo_T"
echo "G: $conteo_G"
echo "C: $conteo_C"
echo "N: $conteo_N"

Contenido de bases en secuencias1.fastq:

A: 643 T: 450 G: 494 C: 470 N: 1408

Ejercicio 10.8 Calcular la distribución de longitudes en una lista de archivos

Calcula la distribución de las longitudes de las secuencias en una lista de archivos fastq. En este ejercicio, deberán discutir qué hace cada uno de los pasos dentro del script y discutirlos con el resto de la clase.

#!/bin/bash
# LongitudArchivos.sh

# Directorio donde se encuentran los archivos FASTQ
directorio="$HOME/unix.bioinfo.rsgecuador-gh-pages/_files"

# Loop para procesar cada archivo FASTQ
for archivo in "$directorio"/*.fastq; do
  if [ -f "$archivo" ]; then
    echo "Procesando archivo: $archivo"

    # Inicializar contador de longitudes
    longitudes=()

    # Leer el archivo línea por línea
    contador=0
    while IFS= read -r linea; do
      ((contador++))
      
      # Solo procesamos las líneas de secuencias (segunda línea de cada bloque)
      if ((contador % 4 == 2)); then
        # Agregar la longitud de la secuencia al array
        longitudes+=(${#linea})
      fi
    done < "$archivo"

    # Calcular la distribución de longitudes
    echo "Distribución de longitudes para $archivo:"
    for longitud in "${longitudes[@]}"; do
      echo -n "$longitud "  # Imprimir todas las longitudes
    done

    echo -e "\nNúmero total de secuencias: ${#longitudes[@]}"

    # Opcional: generar un histograma (puedes ajustarlo según el rango de longitudes)
    echo "Histograma de longitudes de secuencias:"
    for longitud in "${longitudes[@]}"; do
      echo "$longitud"
    done | sort | uniq -c | sort -n  # Ordena y cuenta las frecuencias

  else
    echo "El archivo $archivo no existe."
  fi
done

Ejercicio 10.9 Encontrar la longitud de la secuencia más larga

Determina la longitud de la secuencia más larga entre los archivos FASTQ.

grep -A 1 "^@" *.fastq | grep -v "^@" | grep -v "^--$" | wc -L

La respuesta es: 72

Ejercicio 10.10 Encontrar un motivo específico en un archivo FASTQ

El objetivo de este ejercicio es buscar un motivo específico (un patrón de nucleótidos) dentro de las secuencias de un archivo FASTQ y contar cuántas secuencias contienen ese motivo.

Ejemplo de motivo: Supongamos que el motivo a buscar es AGCT

# Motivo a buscar
motivo="AGCT"

# Buscar el motivo en las secuencias del archivo FASTQ
grep -A 1 "^@" *.fastq | grep -v "^@" | grep -v "^--$" | grep -c "$motivo"

La respuesta es: 2

Y la historia continua ….

15.2 Material suplementario