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.
En estas figuras podemos notar la calidad de las bases para ambos archivos fastq y como éstas disminuyen hacia el final (posición >190 bp) pero también en las primeras bases.
Si vemos el contenido de bases por posición (figuras), podemos notar mucho ruido al principio (1-15) y al final (a partir de la posición 190). Al parecer este ruido es causado por la fragmentación del kit Nextera, que no es tan aleatoria como se sugiere.
Limpieza
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.
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
Si requerimos unir (ensamblar) secuencias pair end, podemos usar PEAR; este programa genera tres tipos de archivos:
Assembled, las lecturas que pudo ensamblar (un archivo).
Discarded, lecturas descartadas (un archivo).
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.
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
SCRIPT nextera_cleaner
En biobacter tenemos ya un script que puede realizar todos estos pasos y otros de una manera mas eficiente: nextera_cleaner.sh como su nombre lo indica, esta diseñado para eliminar adaptadores generados con el kit Nextera de Illumina, eliminar bases con calidad menos a Q30 y secuencias de puras Ns. Además se le pueden indicar cuantas bases eliminar al principio de las secuencias y recortarlas a cierto tamaño eliminando bases al final.
Además permite:
Unir secuencias pair-end (merge) con FLASH.
convertir de fastq a fasta.
analizar con FastQC con los archivos finales.
$ nextera_cleaner file_R1.fastq file_R2.fastq
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.