El control de la calidad de las lecturas normalmente incluye dos pasos:
- Eliminación del extremo 3’ que está por debajo del umbral de calidad
- Detección y eliminación del adaptador
Datos de prueba
Protocolo de análisis (para completar el protocolo debéis entrar en el servidor de docencia)
0) descargar el fichero
wget -c http://bioinfo2.ugr.es/biocomputacion/wp-content/uploads/2017/11/DRR001650_part.fastq_.gz
1) descomprimir el fichero
gunzip DRR001650_part.fastq_.gz
el resultado es un fichero fastq_
2) renombrar el fichero (opcional)
mv nombre_original.fastq_ nombre_nuevo.fastq
3) visualizar parte del fichero
more nombre_nuevo.fastq
4) lanzar fastqc
Para determinar la calidad del conjunto de datos, vamos a usar el programa fastQC que se encuentra instalado en /opt/FastQC
Para lanzar el programa, solo tenemos que proporcionar el fichero
Página de fastQC: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Ayuda de fastqc: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/
fastqc nombre_nuevo.fastq &
El resultado de este comando es un fichero html (con el mismo nombre seguido de _fastqc.html). Para visualizarlo podéis descargarlo a vuestro ordenador y abrirlo con un navegador web (Chrome o Firefox). Explorad el report para ver los distintos parámetros de calidad.
5) detectar el adaptador y eliminar el extremo 3’ con baja calidad
trim_galore nombre_nuevo.fastq
Manual de trim_galore:
http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/trim_galore_User_Guide_v0.3.7.pdf
y del cutadapt (que usa el trim_galore)
https://cutadapt.readthedocs.io/en/stable
Trimgalore generará un fastq con las secuencias “recortadas”, es decir sin el adaptardor y eliminando los reads (o fragmentos de reads) de baja calidad
Cuestiones
- Comparar el resumen de calidad producido por fastqc antes y después de usar trim_galore (presta especial atención a los parámetros que se pueden ver afectados por el trimming: per base quality, adapter content).
- Explorar el parámetro –q del trim_galore (el umbral del PhredScore)
- ¿Que haríamos con los datos de prueba 1?
- Discutir el método empleado por trim_galore (cutadapt)
- Analizar el segundo conjunto de datos (test data 2). ¿Como evaluamos este segundo conjunto?
Preprocesar II
Datos: miRNA test data
Tareas y cuestiones:
- Eliminar el adaptador:
- ¿Cual es la secuencia del adaptador?
- ¿Que porcentage de lecturas empiezan con el 3′ adaptador?
- Buscar secuencias (pista, usa grep. Aunque también puedes hacer un script en python):
- ¿Cuantas veces existe la secuencia TCCTGTACTGAGCTGCCCCGAG en el fichero sin adaptadores?
- ¿Cuantas lecturas contienen TCCTGTACTGAGCTGCCCCGAG como subsecuencia (incluyendo la secuencia exacta)?
- Analizar la distribución de la longitud de las lecturas (frecuencia para longitud 1nt, 2nt, etc)
- ¿Cuantos maximos (incluyendo locales) se pueden distinguir?
- Interpretar los maximos (qué los pueden causar, tipos de RNA, etc)
- Lecturas únicas (lecturas diferentes)
- ¿Cuantas lecturas únicas hay en el fichero?
- ¿A que molécula corresponde la lectura mas frecuente? – ¿y la segunda mas frecuente?
- ¿Cuál es la diferencia entre la más frecuente y la segunda más frecuente?
- ¿Cuántas lecturas únicas tienen una logitud mayor o igual a 18 nt?
- Formato fasta de lecturas únicas
- Generar un multifasta con una entrada con cada lectura
- lanzar un blast local frente a la base de datos miRBase (ver comando abajo)
- Contar el número de lecturas mapeadas a cada secuencia de microRNA
Para hacer un blast local:
blastn -query 'ruta del fasta file' -db hsa_mature -outfmt 6 -out 'fichero de salida' -word_size 9
Explicación de la salida de Blastn
Más cuestiones (a resolver programáticamente)
- ¿Cuantas lecutras unicas se eliminan aplicando un minimo de 2?
- Aplicando un minimo de 5 lecturas, ¿que valor toma el conteo total?
- ¿Cuantas lecturas totales y unicas hay entre 20 y 23 nucleótidos?
- Parseando la salida del Blast, ¿cual es el microRNA más expresado?
Soluciones
Ayuda
- Explorar la opción –length y -o del trim_galore
- Explorar los programas grep y wc disponibles en linea de comando
- la función len() de python