Limpieza archivos Illumina

Limpieza de secuencias pair-end de Illumina

Los equipos de Illumina, si se corrieron con pair-end, generan dos archivos que llevan en algún lugar del nombre R1 y R2, éstos son las secuencias en forwards y en reverse. Puede suceder que el equipo no limpie completamente adaptadores, barcodes (indexes) y bases de baja calidad, por lo que es conveniente analizar ambos archivos fastq con fastQC.

Una buena revisión de los problemas al usar el kit Nextera de Illumina se encuentra aquí.

FastQC

En el directorio donde están ambos archivos fastq podemos correr el comando:

$ fastqc file_R1.fastq file_R2.fastq

FastQC generará archivos con los resultados y podemos consultar el .html correspondiente a cada uno de los dos fastq.

Limpieza

fastx_trimmer

Con esta información podemos programar la limpieza, como casi todas las posiciones tienen una alta calidad (>Q30) solo tendríamos que eliminar las primeras 15 bases y a partir de la 190. Para lo cual podemos correr el siguiente comando para cada uno de los dos archivos:

$ fastx_trimmer -f 15 -l 190 -i file.fastq -o file_trim2.fastq

-f Número de bases a recortar a la izquierda (3').

-l Recortar a partir de que número de base a la derecha (5').

-i Archivo a recortar, usar el que se generó con el script anterior (file_trim.fastq).

-o Nombre del nuevo archivo recortado de salida.

cutadapt

Otra opción mas completa es usar cutadapt para la búsqueda de adaptadores y a la vez limpieza por calidad y ambigüedades:

$ cutadapt -u 15 -U 15 -a CAAGCAGAAGACGGCATACGAGAT -a CTGTCTCTTATACACATCT -A AATGATACGGCGACCACCGAGATCTACAC -A CTGTCTCTTATACACATCT --times 2 -q 30,30 --trim-n -o filename.R1.fastq -p filename.R2.fastq forward.fastq reverse.fastq

Opciones que realiza ANTES de la búsqueda de adaptadores:

-u Recortar bases a la derecha (3') si el valor es positivo, si negativo, entonces se cortan desde la izquierda (5') para la secuencias R1 (forward).

-U igual a -u pero para la segunda secuencia (R2, reverse). Estas bases se cortan ANTES de los adaptadores por lo que hay que ver si eso es lo que se quiere, podría, por lo tanto, no reconocer adaptadores completos.

Búsqueda y eliminación de adaptadores:

-a Secuencia del adaptador en la región 3' de la primer secuencia pareada (R1), puede ponerse varias secuencias de adaptadores repitiendo -a.

-A Secuencia del adaptador en la región 3' de la segunda secuencia pareada (R2).

--times 2 Hace dos pases para eliminar bien los adaptadores.

Opciones que realiza DESPUÉS de la eliminación de adaptadores:

-m Longitud mínima, lecturas menores a este valor serán eliminadas, no necesario para secuencias Illumina.

--trim-n Recorta Ns al final de las secuencias, si existieran.

-q Bases con calidad Phred (Q) a recortar en ambas secuencias si se usa doble (-q 30,30).

Unión de lecturas pair end

PEAR

Si requerimos unir (ensamblar) secuencias pair end, podemos usar PEAR; este programa genera tres tipos de archivos:

1. Assembled, las lecturas que pudo ensamblar (un archivo).

2. Discarded, lecturas descartadas (un archivo).

3. Unassambled, lecturas que no ensambló por que no se traslaparon (dos archivos).

$ pear -f forward.fastq -r reverse.fastq -o sample -j 4

-f Archivo con las lecturas en sentido forward.

-r Archivo con las lecturas en sentido reverse.

-o Un nombre que identifique a la muestra y que será antepuesto a los archivos de salida.

-j Número de núcleos a usar.

FLASh

Otra opción para ensamblar secuencias pair end es flash:

$ flash forward.fastq reverse.fastq -o sample -t 2

-o Un nombre que identifique a la muestra y que será antepuesto (prefix) a los archivos de salida.

-t Número de núcleos a usar.

Arreglo de lecturas pair end

Suele suceder que cuando se limpian las secuencias pair end, y especialmente cuando se hace un recorte o trimming, se elimina toda una secuencia de solo uno de los dos archivos, entonces los archivos quedan "desbalanceados", uno tiene unas secuencias de mas o de menos. Cuando esto sucede, algunos de lo programas en análisis posteriores como SPAdes, no pueden continuar. Por lo tanto debemos corregir este desbalanceo volviendo a "parear" ambos archivos. Esta corrección la podemos hacer con Pairfq:

$ pairfq makepairs -f R1.fastq -r R2.fastq -fp R1_pair.fastq -rp R2_pair.fastq -fs R1_unpair.fastq -rs R2_unpair.fastq --stats

Archivos de entrada:

-f for the file of forward reads

-r for the reverse reads

Archivos de salida (pueden ser cualquier nombre pero que no se repitan, ver el ejemplo arriba):

-fp for the file of paired forward reads

-rp for the file of paired reverse reads

-fs for the file of forward singleton/unpaired reads

-rs for the singleton/unpaired reverse reads

Convertir fastq a fasta

En caso que se necesite convertir los archivos fastq a fasta para continuar con ciertos análisis, se puede hacerlo con un sencillo comando:

$ sed -n '1~4s/^@/>/p;2~4p' file.fastq > file.fasta

Varias muestras

A veces queremos limpiar los archivos fastq de varias muestras que están en subdirectorios, para ésto podemos hacer un script relativamente sencillo que vaya entrando a cada uno de los subdirectorios y ejecute una limpieza con cutadapt, por ejemplo. Esto es útil si ya sabemos que todas las secuencias son similares en calidad y posibles adaptadores contaminantes.

Para estos casos tenemos también un script llamado multiple_fastq_cleaner, solo tenemos que ejecutarlo en el directorio que contenga subdirectorios con los dos archivos fastq (_R1_ y _R2_) que pueden estar comprimidos con terminación .gz

La estructura de este directorio sería similar a:

|_Sample01

| |_Sample01_S001_R1_001.fastq.gz

| |_Sample01_S001_R2_001.fastq.gz

|_Sample02

|_Sample02_S002_R2_001.fastq.gz

|_Sample03_S002_R2_001.fastq.gz

Esta estructura es como sale de los secuenciadores Illumina, donde el nombre de la muestra sólo es la primera parte antes del primer guión bajo (Sample01 y Sample02), el _S001_ lo pone el equipo y corresponde a un numero de muestra de la corrida; el _R1_ o _R2_, son el archivo pair en forward y el reverse, respectivamente.

El script genera una nuevo directorio arriba llamado clean_fastq con carpetas con las secuencia limpias adentro de cada una, ahora el nombre de los archivos sera solo la muestra y el número de pair-end, por ejemplo Sample01.R1.fastq y Sample01.R2.fastq.

Tiene parámetros por default que corte las primeras 15 bases de cada archivo, que usualmente tienen adaptadores, y elimine las secuencias inferiores a Q30 (un probable error en cada 1,000 bases), además elimina terminaciones con varias Ns. Estos parámetros son mas o menos estándar para nuestro equipo Illumina Miniseq.