Proyecto de Técnicas y Herramientas de Bioinformática: Un pequeño mapeador de secuencias

MEI - FIB

Requerimientos

Este proyecto se ha desarrollado utilizando:

• Lenguaje: Python 3.x

• Librerías:

  • numpy: Usado principalmente para implementar de manera sencilla la matriz del algoritmo de programación dinámica.
  • Librerías estándar de Python: argparse, time, etc.

• Archivos necesarios:

  • Referencia: Archivo FASTA que contiene el genoma (por ejemplo, CP011330.1.fasta).
  • Lecturas: Archivo FASTQ que contiene lecturas simuladas o reales a mapear.

• Archivo generado:

  • SAM: Archivo de salida en formato estándar SAM con los resultados del mapeo.

Este script se ha creado para cumplir con los requisitos del proyecto final de la asignatura, implementando conceptos fundamentales vistos en clase como lectura de archivos FASTA/FASTQ, índices de k-mers, y alineamiento mediante programación dinámica.

---

Quick Init

Para ejecutar este script de forma sencilla, ubica los archivos necesarios en la carpeta actual y utiliza el siguiente comando:

python mapper.py -r CP011330.1.fasta -f lecturas.fastq -k 8 -o output.sam

• -r CP011330.1.fasta: Archivo FASTA del genoma usado como referencia.
• -f lecturas.fastq: Archivo FASTQ que contiene las lecturas a mapear.
• -k 8: Tamaño del k-mer utilizado para el índice (puede ajustarse según necesidad).
• -o output.sam: Archivo SAM resultante con los alineamientos.

---

Descripción general

Este proyecto realiza las siguientes acciones:

• Lectura del archivo FASTA con el genoma.
• Lectura de las secuencias de un archivo FASTQ.
• Creación de un índice de k-mers para optimizar la búsqueda.
• Identificación de posiciones candidatas en el genoma mediante el índice.
• Alineamiento aproximado de cada lectura usando programación dinámica, calculando distancias de edición y generando cadenas CIGAR.
• Generación del archivo SAM final con los resultados.

---

Funciones

read_fasta(file_path)

• Propósito:
  Lee un archivo FASTA y retorna la secuencia completa ignorando el encabezado.
• Argumento:
  • file_path (str): Ruta al archivo FASTA.
• Retorno:
  Una cadena (str) que contiene la secuencia concatenada de todas las líneas (sin los encabezados).
• Funcionamiento general:
  Se abre el archivo, se ignoran las líneas que comienzan con ">", y se concatenan las demás líneas para formar la secuencia de referencia.

---

read_fastq(file_path)

• Propósito:
  Lee un archivo FASTQ y retorna una lista de lecturas.
• Argumento:
  • file_path (str): Ruta al archivo FASTQ.
• Retorno:
  Una lista de tuplas, donde cada tupla es (título, secuencia, calidad).
• Funcionamiento general:
  Se lee el archivo en bloques de 4 líneas (título, secuencia, línea separadora y calidad) y se agregan a una lista hasta llegar al final del archivo.

---

build_kmer_index(sequence, k)

• Propósito:
  Construye un índice de k-mers a partir de la secuencia de referencia.
• Argumentos:
  • sequence (str): Secuencia completa (genoma de referencia).
  • k (int): Tamaño del k-mer a utilizar para la indexación.
• Retorno:
  Un diccionario donde cada clave es un k-mer y el valor es una lista de posiciones (enteros) donde ese k-mer aparece en la secuencia.
• Funcionamiento general:
  Se recorre la secuencia con una ventana de tamaño k y se almacena cada k-mer junto con su posición inicial. Esto permite acceder rápidamente a las posiciones donde aparece cada fragmento.

---

get_candidate_positions(kmer_index, read, k)

• Propósito:
  A partir de una lectura, extrae k-mers y utiliza el índice para obtener posiciones candidatas en la referencia.
• Argumentos:
  • kmer_index (dict): Índice de k-mers generado a partir de la referencia.
  • read (str): Secuencia de la lectura que se desea mapear.
  • k (int): Tamaño del k-mer.
• Retorno:
  Un conjunto de posiciones candidatas (enteros) ajustadas según el offset del k-mer en la lectura.
• Funcionamiento general:
  Para cada k-mer de la lectura se consulta el índice; las posiciones encontradas se ajustan (se resta el índice del k-mer en la lectura) y se acumulan en un conjunto para evitar duplicados.

---

dp_align(query, ref_segment)

• Propósito:
  Realiza un alineamiento por programación dinámica entre la secuencia de consulta y un segmento de la referencia.
• Argumentos:
  • query (str): La secuencia de la lectura (consulta).
  • ref_segment (str): Segmento de la referencia de igual longitud que la lectura.
• Retorno:
  Una tupla (distancia, cigar):
  • distancia (int): Distancia de edición (Levenshtein) entre la lectura y el segmento.
  • cigar (str): Cadena CIGAR simple que indica el tipo de operación en cada posición:
    • "M" para coincidencias,
    • "X" para sustituciones,
    • "I" para inserciones y
    • "D" para deleciones.
• Funcionamiento general:
  Se construye una matriz DP de tamaño (len(query)+1) x (len(ref_segment)+1), inicializando la primera fila y columna con penalizaciones de inserción o deleción. Luego se llena la matriz comparando cada base, tomando el mínimo entre coincidencia/sustitución, inserción o deleción. Finalmente, se realiza el retroceso (traceback) para generar la cadena CIGAR correspondiente.

---

write_sam(sam_path, alignments, ref_name, ref_length)

• Propósito:
  Escribe un archivo SAM con los alineamientos obtenidos.
• Argumentos:
  • sam_path (str): Ruta y nombre del archivo SAM de salida.
  • alignments (list): Lista de alineamientos. Cada alineamiento es una tupla con:
    (QNAME, FLAG, POS, MAPQ, CIGAR, SEQ, QUAL).
  • ref_name (str): Nombre de la referencia (usualmente se toma el nombre del archivo FASTA).
  • ref_length (int): Longitud de la secuencia de referencia.
• Retorno:
  No retorna valor; escribe el archivo SAM.
• Funcionamiento general:
  Se escribe la cabecera SAM básica (HD y SQ), y para cada alineamiento se escribe una línea con los campos requeridos según la especificación SAM. Se usa un FLAG de 0 para lecturas mapeadas y 4 para no mapeadas.

---

main()

• Propósito:
  Función principal del script, que integra todas las funcionalidades:
  1. Procesa los argumentos de la línea de comandos.
  2. Lee el genoma de referencia (FASTA) y las lecturas (FASTQ).
  3. Construye el índice de k-mers de la referencia.
  4. Para cada lectura, obtiene posiciones candidatas y realiza alineamiento mediante DP.
  5. Selecciona la mejor posición (con menor distancia de edición y dentro del umbral máximo).
  6. Genera y escribe un archivo SAM con los resultados.
• Argumentos (vía argparse):
  • -r / --reference: Ruta al archivo FASTA del genoma de referencia.
  • -f / --fastq: Ruta al archivo FASTQ con las lecturas.
  • -k / --kmer: Tamaño del k-mer para el índice (valor por defecto: 8).
  • --max_edit: Distancia de edición máxima permitida para considerar un mapeo (valor por defecto: 3).
  • -o / --output: Nombre del archivo SAM de salida (valor por defecto: output.sam).
• Funcionamiento general:
  Se miden los tiempos de ejecución para evaluar el rendimiento. Para cada lectura, se buscan candidatos usando el índice de k-mers y se realiza un alineamiento DP para determinar la similitud. Si se encuentra un mapeo dentro del umbral (max_edit), se registra como lectura mapeada; de lo contrario, se marca como no mapeada. Finalmente, se escribe el archivo SAM con los resultados.

---

Conclusiones

Desarrollar este proyecto como parte de la asignatura "Tècniques i Eines Bioinformàtiques" me ha permitido reforzar conceptos esenciales de bioinformática y programación, especialmente sobre el uso práctico de índices de k-mers y algoritmos de alineamiento dinámico. Además, ha sido una buena introducción a formatos estándar como FASTA, FASTQ y SAM, fundamentales para análisis posteriores. Considero que este mapeador proporciona una base útil en técnicas de bioinformática.