Sesion 2. Análisis de Calidad de Secuencias

Control de calidad, procesamiento y alineamiento de lecturas de secuenciación de alto rendimiento

Los avances en la tecnología de secuenciación están permitiendo a los investigadores secuenciar genomas con una profundidad mayor que nunca. Estos experimentos suelen generar millones de lecturas (reads), las cuales deben ser procesadas, verificadas en su calidad y alineadas antes de que podamos cuantificar la señal genómica de interés y aplicar métodos estadísticos o de aprendizaje automático (machine learning).

Por ejemplo, es posible que desees contar cuántas lecturas se solapan con un conjunto de promotores de tu interés, o cuantificar lecturas de RNA-seq que coinciden con exones.

En este capítulo, introduciremos los fundamentos del procesamiento y control de calidad de las lecturas, y mostraremos cómo realizar estas tareas en R. El control de calidad (QC) y el procesamiento son pasos fundamentales en todos los análisis de secuenciación de alto rendimiento. Por ejemplo, los análisis de RNA-seq, ChIP-seq y BS-seq que se muestran en los Capítulos 8, 9 y 10 requieren estos pasos de QC y procesamiento previos a cualquier análisis posterior. Durante mucho tiempo, las tareas de control de calidad y mapeo estuvieron fuera del dominio de R; sin embargo, hoy en día existen paquetes en R/Bioconductor capaces de realizar estas tareas de manera eficiente.

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Bloque 1: Paquetes generales y de gráficos (menos pesados)
BiocManager::install(c("ggplot2", "plot3D", "pheatmap", "cowplot", "cluster"))

# Bloque 2: Paquetes estadísticos y de reportes
BiocManager::install(c("qvalue", "fastqcr"))

# Bloque 3: Los pesos pesados de bioinformática (instálalos solos)
BiocManager::install("ShortRead")
BiocManager::install("Rqc")
BiocManager::install("QuasR")

library(Rqc)
folder = system.file(package="ShortRead", "extdata/E-MTAB-1147")

# feeds fastq.qz files in "folder" to quality check function
qcRes=rqc(path = folder, pattern = ".fastq.gz", openBrowser=FALSE)

Formatos FASTA y FASTQ

Las lecturas de secuenciación de alto rendimiento generalmente son entregadas por los centros de secuenciación como archivos de texto en un formato llamado “FASTQ” o “fastq”. Este formato se deriva de un formato anterior llamado FASTA.

El formato FASTA

El formato FASTA fue desarrollado como un estándar basado en texto para representar secuencias de nucleótidos o proteínas. Es el formato más común para almacenar genomas de referencia.

Un archivo FASTA se compone de dos elementos principales: 1. Línea de encabezado: Comienza con el símbolo > seguido de un identificador único para la secuencia (como el nombre de un cromosoma o un gen). 2. Línea de secuencia: El contenido de la secuencia biológica en sí.

El formato FASTQ

El formato FASTQ es la extensión natural del FASTA para la secuenciación de alto rendimiento, ya que no solo almacena la secuencia de nucleótidos, sino también la calidad de cada base llamada por el secuenciador.

Cada lectura en un archivo FASTQ se representa en cuatro líneas:

1. Línea 1: Comienza con @ y contiene el identificador de la lectura (información sobre el secuenciador, la celda de flujo y las coordenadas).

2. Línea 2: La secuencia de bases (A, C, G, T, N).

3. Línea 3: Comienza con un símbolo +, y opcionalmente repite el identificador de la línea 1.

4. Línea 4: Una cadena de caracteres que representan la calidad de cada base (puntuación Phred), codificada en formato ASCII.

Detalles del formato FASTA

La primera línea en un archivo FASTA generalmente comienza con el símbolo “>” (mayor que). Esta primera línea se denomina “línea de descripción” y puede contener información descriptiva sobre la secuencia que aparece en las líneas posteriores. Esta descripción puede ser el ID o el nombre de la secuencia, como por ejemplo, nombres de genes.

Aunque es muy poco frecuente, es posible encontrar líneas que comienzan con un “;” (punto y coma). Estas líneas se interpretan como comentarios y pueden contener información descriptiva adicional sobre la secuencia en las líneas siguientes.

El formato FASTQ: Secuencia y Calidad

Una extensión del formato FASTA es el formato FASTQ. Este formato está diseñado para manejar las métricas de calidad de base generadas por las máquinas de secuenciación. En este formato, tanto la secuencia como las puntuaciones de calidad se representan mediante caracteres ASCII únicos.

Estructura detallada del formato FASTQ

Cada lectura en un archivo FASTQ se define por cuatro líneas específicas:

• Línea 1: Comienza con el carácter @ seguido de un identificador de secuencia y una descripción opcional. Esta línea es utilizada por la tecnología de secuenciación y suele contener información específica del equipo, como el ID de la celda de flujo (flow cell), números de carril (lane numbers) e información sobre pares de lecturas (read pairs).
• Línea 2: Contiene las letras de la secuencia biológica (A, C, G, T, N).
• Línea 3: Comienza con el carácter +; marca el final de la secuencia y, opcionalmente, puede repetir el identificador de la secuencia de la Línea 1.
• Línea 4: Codifica los valores de calidad para la secuencia de la Línea 2. Debe contener exactamente el mismo número de símbolos que letras tiene la secuencia. Cada símbolo corresponde a una puntuación de calidad específica.

Puntuaciones de calidad Phred (Phred Quality Scores)

Aunque pueden existir diferentes definiciones para las puntuaciones de calidad, el estándar de facto en el campo es utilizar las “puntuaciones de calidad Phred”. Estas puntuaciones representan la probabilidad de que una base haya sido identificada incorrectamente por el secuenciador.

Formalmente, la puntuación se define mediante la fórmula:

\[Q_{phred} = -10 \log_{10}(e)\]

Donde \(e\) es la probabilidad de que la base sea errónea.

Debido a que la puntuación está en una escala logarítmica negativa: * Una puntuación más alta significa que es menos probable que la base sea errónea. * Por ejemplo, una puntuación \(Q = 30\) (\(10^{-3}\)) significa que hay una probabilidad del 0.1% de que la base sea incorrecta (99.9% de precisión).

1.2 Control de calidad de las lecturas de secuenciación

Las tecnologías de secuenciación suelen producir identificaciones de bases (basecalls) con calidades variables. Además, pueden surgir problemas específicos de la muestra durante la carrera de secuenciación, como la contaminación por adaptadores. El procedimiento estándar dicta que se debe verificar la calidad de las lecturas e identificar posibles problemas antes de proceder con análisis posteriores. Realizar un control de calidad adecuado y tomar decisiones informadas para el análisis downstream puede influir significativamente en los resultados finales de tu proyecto.

A continuación, te guiaremos a través de los pasos del control de calidad utilizando el paquete Rqc.

Uso del paquete Rqc

Primero, necesitamos proporcionar los archivos FASTQ a la función rqc() para obtener un objeto que contenga los resultados relacionados con la calidad de la secuencia. En este ejemplo, utilizaremos archivos FASTQ de muestra incluidos en el paquete ShortRead.

Calidad de secuencia por base/ciclo

Una vez que hemos obtenido el objeto qcRes (o qc_results en nuestro código previo), podemos graficar diversas métricas de calidad para nuestros archivos FASTQ.

La primera métrica que graficaremos es la “calidad de secuencia por base/ciclo”. Este gráfico muestra las puntuaciones de calidad de todas las bases en cada posición (ciclo) de las lecturas.

rqcCycleQualityBoxPlot(qcRes)

¿Por qué hablamos de “ciclos”?

En nuestro gráfico, el eje X está etiquetado como “ciclo” (cycle). Esto se debe a que, en cada “ciclo” de secuenciación, se añade un nucleótido marcado fluorescentemente para complementar la secuencia molde, y la máquina de secuenciación identifica qué nucleótido se ha incorporado. Por lo tanto, los ciclos corresponden a las bases/nucleótidos a lo largo de la lectura, y el número total de ciclos es equivalente a la longitud de la lectura.

Toma de decisiones: ¿Recortar o no recortar?

Las secuencias largas pueden presentar una degradación en la calidad hacia los extremos de las lecturas. Observar la distribución de la calidad sobre las posiciones de las bases nos ayuda a decidir si debemos realizar un recorte (trimming) hacia el final de las mismas o no.

• Muestra de buena calidad: Tendrá puntuaciones medianas de calidad por base superiores a 28.
• Muestra problemática: Si las puntuaciones caen por debajo de 20 hacia los extremos, deberías considerar seriamente recortar las lecturas.

El recorte de estas bases de baja calidad es un paso crítico, ya que los errores en las llamadas de bases al final de la lectura pueden impedir que los alineadores (como Rsubread) encuentren la posición correcta en el genoma de referencia, disminuyendo la tasa de mapeo de tu experimento de metagenómica.

Contenido de secuencia por base/ciclo

El contenido de secuencia por base muestra las proporciones de cada nucleótido (A, C, G, T) en cada posición de la lectura. En una librería de secuenciación aleatoria ideal, no debería existir un sesgo de nucleótidos; por lo tanto, las líneas de cada base deberían ser casi paralelas entre sí a lo largo de todo el gráfico.

El siguiente código muestra cómo generar este gráfico utilizando el objeto de resultados de Rqc:

rqcCycleBaseCallsLinePlot(qcRes)

Sesgos esperados en la composición de secuencias

Es importante notar que ciertos tipos de librerías de secuenciación pueden producir una composición de secuencias inherentemente sesgada.

1. Sesgo en RNA-Seq (Inicio de la lectura)

En experimentos de RNA-Seq, es muy común observar un sesgo en los primeros nucleótidos de las lecturas. Esto ocurre debido a la unión de cebadores aleatorios (random primers) al inicio de las lecturas durante la preparación de la librería de cDNA. Estos cebadores no son verdaderamente aleatorios, lo que genera una variación en la composición de bases al comienzo de la lectura.

A pesar de que los experimentos de RNA-seq suelen presentar estos sesgos, esto generalmente no afecta la capacidad de medir la expresión génica con precisión.

2. Sesgo en Secuenciación con Bisulfito (BS-Seq)

Por otro lado, algunas librerías son inherentemente sesgadas debido al tratamiento químico de la muestra. Por ejemplo, en los experimentos de secuenciación con bisulfito (utilizados para estudiar la metilación del ADN), la mayoría de las citosinas (C) se convierten en timinas (T) durante el proceso.

Esto creará una diferencia marcada y constante en las proporciones de las bases C y T a lo largo de toda la lectura. Sin embargo, este tipo de diferencia es normal y esperada en los experimentos de secuenciación con bisulfito, y no debe confundirse con una falla en la calidad de la secuenciación.

Gráfico de frecuencia de lecturas (Read frequency plot)

Este gráfico muestra el grado de duplicación de cada lectura en la librería. A continuación, mostramos cómo generar este gráfico y el resultado se puede observar en la siguiente figura.

Un nivel alto de duplicación (lecturas no únicas) suele indicar un sesgo de enriquecimiento. Esto puede ser causado por duplicados técnicos derivados de artefactos de la PCR. La PCR es un paso común en la preparación de librerías que crea muchas copias de un mismo fragmento de secuencia.

Interpretación según el tipo de experimento:

• En datos de RNA-seq: La proporción de lecturas no únicas puede superar el 20%. Sin embargo, estas duplicaciones pueden deberse simplemente a que ciertos genes se expresan a niveles muy altos. Esto significa que habrá muchas copias de los transcritos y, por ende, muchas copias del mismo fragmento. Dado que no podemos estar seguros de si estas lecturas duplicadas se deben a un sesgo de la PCR o a una alta tasa de transcripción, no se deben eliminar las lecturas duplicadas en los análisis de RNA-seq.
• En experimentos de ChIP-seq: Las lecturas duplicadas tienen una mayor probabilidad de ser el resultado de un sesgo de la PCR.

rqcReadFrequencyPlot(qcRes)

Otras métricas de calidad y herramientas de QC

La presencia de k-mers sobre-representados a lo largo de las lecturas puede servir como un control adicional. Si aparecen este tipo de secuencias, podría ser un indicativo de contaminación por adaptadores, los cuales deberían ser recortados. Los adaptadores son secuencias conocidas que se añaden a los extremos de las lecturas durante la preparación de la librería.

Este tipo de contaminación también suele ser visible en los gráficos de “contenido de secuencia por base”. Si conoces las secuencias de los adaptadores, puedes buscarlas en los extremos de las lecturas y eliminarlas mediante un proceso de recorte (trimming).

FastQC: El estándar de la industria

La herramienta más popular para el control de calidad de secuenciación es FastQC (Andrews 2010), la cual está escrita en Java. Esta herramienta genera los gráficos que hemos descrito anteriormente, además de gráficos específicos de sobre-representación de k-mers y de adaptadores.

El paquete de R fastqcr puede ejecutar esta herramienta de Java y producir informes y gráficos basados en R. Básicamente, el paquete llama a la herramienta Java y procesa sus resultados automáticamente. A continuación, mostramos cómo realizar este proceso:

library(fastqcr)

# install the FASTQC java tool
fastqc_install()

# call FASTQC and record the resulting statistics
# in fastqc_results folder
fastqc(fq.dir = folder,qc.dir = "fastqc_results")

Generación de reportes con fastqcr

Una vez que hemos ejecutado FastQC sobre nuestros archivos FASTQ, podemos importar los resultados a R para construir gráficos personalizados o informes completos.

La función qc_report() es especialmente útil, ya que permite crear un informe basado en R Markdown a partir de los resultados generados por FastQC. Esto facilita la documentación de tu análisis de metagenómica y permite compartir los resultados de calidad de forma clara y organizada.

A continuación, se muestra cómo generar este reporte:

# view the report rendered by R functions
qc_report(qc.path="fastqc_results", 
          result.file="reportFile", preview = TRUE)

Visualización personalizada de resultados de FastQC

Como alternativa a la creación de un informe automático, podemos importar los resultados directamente a R utilizando la función qc_read(). Esto nos da la flexibilidad de seleccionar y visualizar únicamente las métricas específicas que nos interesen mediante la función qc_plot().

Este enfoque es ideal si deseas integrar gráficos específicos de calidad dentro de un análisis más amplio o si quieres comparar métricas particulares entre tus muestras de metagenómica.

A continuación, mostramos cómo cargar los datos y generar gráficos individuales:

# read QC results to R for one fastq file
qc <- qc_read("fastqc_results/ERR127302_1_subset_fastqc.zip")

# make plots, example "Per base sequence quality plot"
qc_plot(qc, "Per base sequence quality")

Otras alternativas en Bioconductor para reportes de calidad

Además de fastqcr, existen otros paquetes potentes dentro del ecosistema de Bioconductor que pueden generar informes de calidad de manera similar a FastQC, aunque presentan algunas diferencias en el contenido y la cantidad de los gráficos generados:

• Rqc: Ofrece la función Rqc::rqcReport, que genera reportes interactivos y es muy eficiente para procesar múltiples muestras en paralelo.
• QuasR: Incluye la función QuasR::qQCReport. Este paquete es una solución “todo en uno” que no solo hace control de calidad, sino también alineamiento.
• systemPipeR: A través de la función systemPipeR::seeFastq, permite una integración fluida en flujos de trabajo automatizados (pipelines) de análisis genómico.
• ShortRead: Es uno de los paquetes base para el manejo de archivos FASTQ en R y cuenta con la función ShortRead::report para generar resúmenes de calidad.

Cada una de estas herramientas tiene sus propias ventajas. Por ejemplo, algunas son más rápidas para archivos extremadamente grandes, mientras que otras ofrecen gráficos más estéticos o interactivos. La elección dependerá de las necesidades específicas de tu proyecto de metagenómica y de qué tan integrado quieras que esté el control de calidad con el resto de tu código.

7.3 Filtrado y recorte de lecturas (Filtering and trimming)

Basado en los resultados del control de calidad, es posible que desees recortar o filtrar las lecturas. El control de calidad puede haber mostrado una cantidad considerable de lecturas con puntuaciones de calidad bajas. Es probable que estas lecturas no se alineen correctamente debido a posibles errores en la identificación de bases (base calling), o podrían alinearse en lugares incorrectos del genoma. Por lo tanto, es recomendable eliminar estas lecturas de tu archivo FASTQ.

Otro escenario común es que ciertas partes de las lecturas necesiten ser recortadas para poder alinearse. Existen dos situaciones principales: 1. Presencia de adaptadores: Secuencias técnicas que pueden aparecer en cualquiera de los lados de la lectura. 2. Errores técnicos: Problemas que provocan una disminución de la calidad de las bases hacia los extremos de las lecturas.

En ambos casos, esa porción de la lectura debe ser recortada para que la lectura pueda alinearse (o alinearse mejor) al genoma. A continuación, mostraremos cómo usar el paquete QuasR para recortar las lecturas. Otros paquetes como ShortRead también tienen capacidades para filtrar y recortar; sin embargo, la función QuasR::preprocessReads() proporciona una interfaz única para múltiples posibilidades de preprocesamiento.

Con esta función podemos: * Identificar secuencias de adaptadores y eliminarlas. * Eliminar lecturas de baja complejidad (que contienen secuencias repetitivas). * Recortar el inicio o el final de las lecturas por una longitud predefinida.

A continuación, configuraremos las rutas de los archivos FASTQ y los filtraremos basándonos en su longitud y en si contienen o no el carácter “N” (base no identificada). Con la misma función, también recortaremos 3 bases del final de las lecturas y eliminaremos segmentos del inicio si coinciden con la secuencia “ACCCGGGA”.

library(QuasR)

# obtain a list of fastq file paths
fastqFiles <- system.file(package="ShortRead",
                          "extdata/E-MTAB-1147",
                          c("ERR127302_1_subset.fastq.gz",
                            "ERR127302_2_subset.fastq.gz")
)

# defined processed fastq file names
outfiles <- paste(tempfile(pattern=c("processed_1_",
                              "processed_2_")),".fastq",sep="")

# process fastq files
# remove reads that have more than 1 N, (nBases)
# trim 3 bases from the end of the reads (truncateEndBases)
# Remove ACCCGGGA patern if it occurs at the start (Lpattern)
# remove reads shorter than 40 base-pairs (minLength)
preprocessReads(fastqFiles, outfiles, 
                nBases=1,
                truncateEndBases=3,
                Lpattern="ACCCGGGA",
                minLength=40)

##                  ERR127302_1_subset.fastq.gz ERR127302_2_subset.fastq.gz
## totalSequences                         20000                       20000
## matchTo5pAdapter                        3804                        3398
## matchTo3pAdapter                           0                           0
## tooShort                                 355                         426
## tooManyN                                   4                           1
## lowComplexity                              0                           0
## totalPassed                            19641                       19573

Filtrado personalizado con ShortRead

Como hemos mencionado, el paquete ShortRead cuenta con funciones de bajo nivel de las cuales también depende QuasR::preprocessReads(). Podemos utilizar estas funciones de bajo nivel para filtrar lecturas de formas que no son posibles mediante la función estándar de QuasR.

A continuación, mostraremos cómo leer un archivo FASTQ y filtrar las lecturas de modo que eliminemos aquellas donde cada una de sus puntuaciones de calidad esté por debajo de 20. Este es un filtro mucho más estricto y específico.

library(ShortRead)

# obtain a list of fastq file paths
fastqFile <- system.file(package="ShortRead",
                          "extdata/E-MTAB-1147",
                          "ERR127302_1_subset.fastq.gz")

# read fastq file
fq = readFastq(fastqFile)

# get quality scores per base as a matrix
qPerBase = as(quality(fq), "matrix")

# get number of bases per read that have quality score below 20
# we use this
qcount = rowSums( qPerBase <= 20) 

# Number of reads where all Phred scores >= 20
fq[qcount == 0]

## class: ShortReadQ
## length: 10699 reads; width: 72 cycles

Guardar las lecturas filtradas

Una vez que hemos aplicado nuestros criterios de filtrado personalizados y estamos satisfechos con la calidad de nuestras lecturas en el objeto rfq_filtrado, el último paso es exportar estos datos a un nuevo archivo físico.

Para ello, utilizamos la función ShortRead::writeFastq(). Esto creará un nuevo archivo FASTQ en tu servidor que contendrá únicamente las secuencias que pasaron tus filtros de calidad.

# Definimos el nombre del archivo
archivo_salida <- paste(fastqFile, "Qfiltered", sep="_")

# Si el archivo ya existe de un intento anterior, lo borramos
if(file.exists(archivo_salida)) { 
  file.remove(archivo_salida) 
}

## [1] TRUE

# write out fastq file with only reads where all 
# quality scores per base are above 20
writeFastq(fq[qcount == 0], archivo_salida)

Procesamiento eficiente de archivos grandes (Streaming)

Dado que los archivos FASTQ pueden ser extremadamente grandes, a menudo no es factible cargar un archivo de 30 Gigabytes directamente en la memoria RAM. Una forma más eficiente de gestionar la memoria consiste en leer el archivo pieza por pieza (por bloques).

Con este enfoque, podemos realizar nuestras operaciones de filtrado en cada bloque, guardar la parte filtrada en el disco y luego cargar un nuevo bloque. Afortunadamente, esto es posible en R utilizando la función ShortRead::FastqStreamer(). Esta función permite la transmisión (streaming) del archivo FASTQ en fragmentos, que son bloques del archivo con un número predefinido de lecturas.

Podemos acceder a los bloques sucesivos mediante la función yield(). Cada vez que llamamos a la función yield() después de haber abierto el archivo con FastqStreamer(), se cargará una nueva parte del archivo en la memoria.

A continuación, mostramos cómo implementar este flujo de trabajo para procesar archivos masivos sin agotar la RAM de tu servidor:

# Usamos un nombre ligeramente distinto para no chocar con el bloque anterior
archivo_stream <- paste(fastqFile, "Qfiltered_stream", sep="_")

# Si el archivo ya existe de un intento anterior, lo borramos
if(file.exists(archivo_stream)) { 
  file.remove(archivo_stream) 
}

## [1] TRUE

# set up streaming with block size 1000
f <- FastqStreamer(fastqFile, readerBlockSize=1000) 

# we set up a while loop to call yield() function to go through the file
while(length(fq <- yield(f))) {
  
    # remove reads where all quality scores are < 20 
    qPerBase = as(quality(fq), "matrix")
    qcount = rowSums( qPerBase <= 20) 
 
    # write fastq file with mode="a" (append)
    writeFastq(fq[qcount == 0], archivo_stream, mode="a")
}

Mapeo y alineamiento de lecturas al genoma

Tras el control de calidad y el preprocesamiento opcional, las lecturas están listas para ser mapeadas o alineadas al genoma de referencia. Este proceso consiste simplemente en encontrar el origen más probable de cada lectura en el genoma. Dado que pueden existir errores en la secuenciación y mutaciones en los genomas, es posible que no encontremos coincidencias exactas.

Una característica crucial de los algoritmos de alineamiento es su capacidad para tolerar posibles desajustes (mismatches) entre las lecturas y el genoma de referencia. Además, se requieren algoritmos y estructuras de datos eficientes para que el alineamiento se complete en un tiempo razonable.

Índices genómicos

Los métodos de alineamiento suelen crear estructuras de datos para almacenar y buscar eficientemente las lecturas coincidentes en el genoma. Estas estructuras se denominan índices genómicos, y su creación es el primer paso para el alineamiento.

Existen dos tipos principales de métodos basados en cómo se crean estos índices:

1. Tablas de Hash (Hash Tables): Funcionan como tablas de búsqueda simples donde todos los k-mers posibles apuntan a ubicaciones en el genoma. La idea general es que los k-mers superpuestos construidos a partir de una lectura pasan por esta tabla. Cada k-mer apunta a ubicaciones potenciales y, finalmente, la ubicación de la lectura se obtiene optimizando la cadena de k-mers según sus distancias en el genoma y en la lectura.
2. Transformación de Burrows-Wheeler (BWT): Estos algoritmos construyen índices creando una transformación compacta y buscable de todo el genoma. Aunque los detalles técnicos son complejos, estas herramientas (como BWA, Bowtie1/2 y SOAP) proporcionan un alineamiento más rápido y consumen menos memoria.

Alineamiento en R

El mapeo de lecturas en R puede realizarse con los paquetes gmapR, QuasR, Rsubread y systemPipeR. En esta sección, demostraremos el mapeo con QuasR, el cual utiliza el paquete Rbowtie (un adaptador para el alineador Bowtie).

Uso de qAlign()

Para realizar el alineamiento, utilizaremos la función qAlign(), la cual requiere dos argumentos obligatorios: 1. Archivo del genoma: En formato FASTA o como un paquete de tipo BSgenome. 2. Archivo de muestras (Sample file): Un archivo de texto que contiene las rutas a los archivos FASTQ y los nombres de las muestras.

Ejemplo del archivo de muestras: El archivo debe tener un formato de tabla como este:

FileName	SampleName
chip_1_1.fq.bz2	Sample1
chip_2_1.fq.bz2	Sample2

library(QuasR)

# copy example data to current working directory
file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

# genome file in fasta format
genomeFile <- "extdata/hg19sub.fa"

# text file containing sample names and fastq file paths
sampleFile <- "extdata/samples_chip_single.txt"

# create alignments 
proj <- qAlign(sampleFile, genomeFile)

¿Qué sucede dentro de qAlign()?

Es importante explicar qué ocurre internamente, ya que la función qAlign() hace que todo parezca muy sencillo. Esta función ha sido diseñada para ser intuitiva: por ejemplo, crea automáticamente un índice del genoma si este no existe, y buscará índices previos antes de intentar generar uno nuevo.

En el ejemplo anterior, solo proporcionamos dos argumentos: un archivo de texto con los nombres de las muestras y las rutas de los archivos FASTQ, y un archivo del genoma de referencia. Sin embargo, esta función posee muchos “ajustes” (knobs) y puedes cambiar su comportamiento suministrando diferentes argumentos para modificar la forma en que actúa Bowtie.

Características clave de qAlign():

• Parámetros de alineamiento: Puedes pasar parámetros específicos a Bowtie utilizando el argumento alignmentParameter.
• Optimización automática: La función está optimizada para diferentes tipos de problemas de secuenciación y selecciona los parámetros de alineamiento de forma automática. Está diseñada para funcionar en tareas de alineamiento y cuantificación para:
  • RNA-seq (Expresión génica).
  • ChIP-seq (Unión de proteínas al ADN).
  • Secuenciación de ARN pequeños (small-RNA).
  • Bisulfite sequencing (Metilación del ADN).
  • Análisis de alelo-específicos.

Nota importante: Si decides cambiar los parámetros por defecto de Bowtie, hazlo únicamente para problemas de alineamiento sencillos como ChIP-seq o RNA-seq. Para experimentos más complejos, los valores predeterminados suelen ser la opción más segura y eficiente.

Procesamiento posterior de las lecturas alineadas

Una vez completado el alineamiento, es posible que se requiera un procesamiento adicional. Sin embargo, estos pasos suelen depender específicamente del protocolo de secuenciación utilizado.

Por ejemplo: * Metilación (Bisulfite sequencing): Es necesario contar los desajustes (mismatches) de C a T para determinar los niveles de metilación. * Expresión génica (RNA-seq): Se deben contar las lecturas que se solapan con los transcritos conocidos para cuantificar la abundancia de cada gen.

Estas tareas de procesamiento posterior pueden realizarse mediante software especializado relacionado con el alineamiento o, en algunos casos, directamente dentro de R.