Skip to main content Link Menu Expand (external link) Document Search Copy Copied
Curso de transcriptomica

Introducción a datos NGS

Que son los datos NGS? Son las secuencias tal cual salen de la plataforma de secuenciación (Illumina, IonTorrent, PacBio, entre otros). Es decir los reads (lecturas).

Los datos crudos (sin importar la plataforma y el método de laboratorio utilizado) conllevan cierto nivel de basura, es decir:

  • reads de baja calidad que deben ser descartados por completo o en parte (trimmed) antes de proceder a los análisis biológicos de los datos.

  • Secuencias sobrerepresentadas, por ejemplo dímeros de adaptadores.

La longitud y distribución de la calidad de los reads varían de plataforma a plataforma, pero también el tipo de error más común:

  • Illumina: errores puntuales en la asignación de un nucleótido
  • IonTorrent y símiles: se les complica determinar el número correcto de nucleótidos en cadenas como AAAAA.

Dado que los datos crudos son muy pesados (espacio de disco) y que buena parte de los datos crudos son basura, en general los datos crudos a este nivel no se guardan en repositorios públicos (e.g. SRA) hasta que hayan pasado por los filtros del pre-procesamiento.

Como veremos más adelante, los filtros de pre-procesamiento ayudan a identificar las buenas secuencias de la basura a partir de su calidad asociada.

Información en los archivos FASTQ

Representación de secuencias

En cómputo las secuencias de ADN son una string (cadena) de caracteres.

  • Secuencias genómicas

{A,C,G,T}+

  • Secuencias mRNA

{A,C,G,U}+

Secuencias simples: FASTA/Pearson Format

Línea 1: información de la secuencia

Línea 2: la secuencia.


>gi|365266830|gb|JF701598.1| Pinus densiflora var. densiflora voucher Pf0855 trnS-trnG intergenic spacer, partial sequence; chloroplast
GCAGAAAAAATCAGCAGTCATACAGTGCTTGACCTAATTTGATAGCTAGAATACCTGCTGTAAAAGCAAG
AAAAAAAGCTATCAAAAATTTAAGCTCTACCATATCTTCATTCCCTCCTCAATGAGTTTGATTAAATGCG
TTACATGGATTAGTCCATTTATTTCTCTCCAATATCAAATTTTATTATCTAGATATTGAAGGGTTCTCTA
TCTATTTTATTATTATTGTAACGCTATCAGTTGCTCAAGGCCATAGGTTCCTGATCGAAACTACACCAAT
GGGTAGGAGTCCGAAGAAGACAAAATAGAAGAAAAGTGATTGATCCCGACAACATTTTATTCATACATTC
AGTCGATGGAGGGTGAAAGAAAACCAAATGGATCTAGAAGTTATTGCGCAGCTCACTGTTCTGACTCTGA
TGGTTGTATCGGGCCCTTTAGTTATTGTTTTATCAGCAATTCGCAAAGGTAATCTATAATTACAATGAGC
CATCTCCGGAGATGGCTCATTGTAATGATGAAAACGAGGTAATGATTGATATAAACTTTCAATAGAGGTT
GATTGATAACTCCTCATCTTCCTATTGGTTGGACAAAAGATCGATCCA
Fastq: formato para secuencias de secuenciación de siguiente generación

Secuencia fasta + detalles calidad de la información (la Q es de Quality).

  • Línea 1: Encabezado (Header): comienza con @. Sigue el identificador (identifier). Si son datos crudos contiene info del secuenciador que identifica a esta secuencia y el read pair (/1 o /2), si son datos ya procesados en SRA contiene una descripción de la secuencia.

  • Línea 2: la secuencia.

  • Línea 3: Comienza con +. Puede ser sólo el símbolo + o repetir la info del Header.

  • Línea 4: Información de la calidad de secuenciación de cada base. Cada letra o símbolo representa a una base de la secuencia codificado en formato ASCII.

La info de calidad se codifica en ASCII porque esto permite en 1 sólo caracter codificar un valor de calidad. Por lo tanto la línea 2 y la 4 tienen el mismo número de caracteres.

Ordenados de menor a mayor estos son los caracteres ASCII usados para representar calidad en FASTQ:

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~

Dependiendo del tipo y versión de plataforma de secuenciación se toman diferentes caracteres (pero sin desordenarlos):

(Tomé la imágen de aquí)

Pero en todos los casos el valor máximo de calidad es = ~40 y los valores < 20 se consideran bajos.

Ojo: @ y + están dentro de los caracteres ASCII utilizados para codificar la calidad, lo que hace que usar grep en archivos FASTQ pueda complicarse, ojo.

Ejemplos:

Ejemplo de datos FASTQ recién salidos de Illumina:

@HWI-ST999:102:D1N6AACXX:1:1101:1235:1936 1:N:0:
ATGTCTCCTGGACCCCTCTGTGCCCAAGCTCCTCATGCATCCTCCTCAGCAACTTGTCCTGTAGCTGAGGCTCACTGACTACCAGCTGCAG
+
1:DAADDDF<B<AGF=FGIEHCCD9DG=1E9?D>CF@HHG??B<GEBGHCG;;CDB8==C@@>>GII@@5?A?@B>CEDCFCC:;?CCCAC

Y uno más:

@OBIWAN:24:D1KUMACXX:3:1112:9698:62774 1:N:0:
TAATATGGCTAATGCCCTAATCTTAGTGTGCCCAACCCACTTACTAACAAATAACTAACATTAAGATCGGAAGAGCACACGTCTGAACTCAGTCACTGACC
+
CCCFFFFFHHHHHIJJJJJJJJJJJJIIHHIJJJJJJJJJJJJJJJJJJJJIJJJJJJIJJJJIJJJJJJJHHHHFDFFEDEDDDDDDDDDDDDDDDDDDC

¿Quieres saber cuáles son las partes del Header? Clic aquí. Y sí, el último ejemplo es real y por lo tanto hay un Illumina HiSeq2000 que se llama Obiwan :)

Ejemplo de datos FASTQ del SRA:

@SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC
+SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9IC

Los datos FASTQ típicamente están comprimidos en formato gzip (.gz) o tar (.tar.gz o .tgz).

Análisis de datos NGS

CONTROL CALIDAD

Antes de saltar a filtrar tus datos con filtros de calidad que la terminal ejecute muy obediente, lo mejor es ver algunos gráficos básicos que nos dicen mucho más que una serie semi-eterna de caracteres ASCII.

FASTQC es un programa que hace una serie de análisis básicos y estándar de calidad. La mayoría de las empresas de secuenciación efectúan este análisis y te mandan los resultados junto con tus datos crudos.

Los análisis de FASTQC son útiles para identificar problemas que pudieron surgir durante el laboratorio o durante la secuenciación.

El análisis de FASTQC consiste en los siguientes campos:

Notas importantes:

  • FASTQ automáticamente dice si nuestra muestra “pasó” o “falló” la evaluación. Sin embargo debemos tomar esto dentro del contexto de lo que esperamos de nuestra librería, ya que FASTQ espera una distribución al diversa y al azar de nucleótidos, lo que puede no cumplirse en algunos protocolos.

Vamos a la página de dicho programa a ver ejemplos de:

¿Qué que son los dímeros de adaptadores?

Los adaptadores se ligan al ADN de nuestras muestras en un paso de ligación, sin embargo, también pueden ligarse entre sí y luego pegarse a la flow cell (que lo traduzca quién sepa cómo). Resultado: son secuenciados pero no proven datos útiles, simplemente la secuencia de los adaptadores repetida muchas veces. Adelante veremos cómo lidiar con ellos bioinformáticamente, pero se recomienda intentar deshacerse de ellos desde el laboratorio (con pequeños, pequeños imanes como Agencourt o símiles de otras marcas).

¿Qué tanto importa el análisis FASTQC?

Mucho, a partir del análisis FASTQC es que decidirás si tu secuenciación fue exitosa y qué parámetros de pre-procesamiento deberás utilizar para deshacerte del ruido y quedarte con datos limpios.

Escoger los parámetros adecuados de pre-procesamiento es vital ya que todas las corridas de secuenciación son diferentes. Además entender bien tu FASTQC puede permitirte rescatar datos usables incluso dentro de una mala corrida.

fastqc -t 2 ngsexample/data/*fastq.gz -o ngsexample/results/quality/

FILTRADO DE READS

Pre-procesamiento

El pre-procesamiento se refiere al filtrado y edición de los datos crudos para quedarnos con los datos limpios, que son los que se analizarán para resolver preguntas biológicas.

El input son archivos .fastq y el output son también archivos .fastq (o más posiblemente sus versiones comprimidas).

El pre-procesamiento por lo común incluye los siguientes pasos:

Recortar secuencias por calidad (Sequence Quality Trimming)

Recortar (quitar) las bases de baja calidad de la secuencia. En Illumina por lo general se trata de las últimas bases (-3’ end). Cuántas bases cortar puede determinarse tras la inspección visual de análisis FASTQC o automáticamente con un parámetro de calidad.

Factor a considerar: algunos programas de ensamblado requieren que las secuencias a ensamblar sean del mismo largo, por lo que si ocuparás uno de esos programas es necesario recortar todos tus datos al mismo largo (incluso si se secuenciaron en lanes distintas y una tiene buena calidad).

Recortar secuencias (Trimming)

Recortar (quitar) x bases de la secuencia porque no nos interesa conservarlas para el análisis (por ejemplo barcodes o sitios de restricción).

Filtrar secuencias por calidad

Remueve del conjunto de datos todas las secuencias que estén por debajo de un mínimo de calidad (número de bases con calidad <x, promedio de calidad <x y símiles).

Quitar adaptadores

Busca la secuencia de los adaptadores y los recorta de las secuencias finales. También es posible limitar las secuencias finales a sólo aquellas con un adaptador especificado (en vez de otro que pudiera ser contaminación).

Filtrar artefactos

Detecta primers de PCRs, quimeras y otros artefactos y los desecha de los datos finales.

Separar por barcodes “demultiplexear” (demultiplexing)

Identifica las secuencias que contienen uno o más barcodes (también llamado índices), que son secuencias cortas (4-8 bp por lo general) que se incluyen en los adaptadores y que son únicos por muestra. Esto permite identificar y separar todas las secuencias pertenecientes a una muestra de otra secuenciadas al mismo tiempo.

Requiere que le demos la lista de barcodes utilizados y en qué extremo se localizan. Muchos programas tendrán como output un archivo llamado algo como GATCATGC.fastq.gz, donde se encuentran todas las secuencias del barcode GATCATGC. El nombre de tu muestra deberás ponerlo en un paso subsecuente.

Ojo Tu lista barcodes-nombremuestra es de la info más valiosa de tu proyecto, no la pierdas.

Paired end merging

Si se realizó secuenciación Illumina a ambos lados (pair end) es posible unir las lecturas si se detecta que coinciden (aunque sea parcialmente). Esto permite corregir errores de secuenciación al tomar la base de la lectura de mayor calidad.

Remover otras secuencias no deseadas

Busca secuencias no deseadas, como genoma de E. coli, restos de PhiX o partes del genoma que no son de interés (e.g. cpDNA).

Trim_galore 

#trim_galore 
trim_galore --quality 20 --fastqc --length 25 --output_dir ngsexample/results/trimmed/ ngsexample/data/*.gz

GTF y GFF

“Genetic transfer format” y “Genomic feature format” Ref

Sirven para representar características genéticas con una función anotada (gen, mRNA, rRNA, tRNA, etc).

Contenido del formato GTF: seqname(#chr) program feature start end score strand frame attribute(gene_id; txpt_id, etc)

Cada intervalo toma una línea, con info en difernetes columnas. Columnas 1-9 separadas por tab y campos dentro de la col 9 separados por espacio. La columna 9 es compuesta y puede tener varios atributos: (mínimo el identificador id del gen), separados por ;. Las coordenadas son 1-based

Ejemplo:

1 transcribed_unprocessed_pseudogene  gene        11869 14409 . + . gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; 
1 processed_transcript                transcript  11869 14409 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; gene_name "DDX11L1"; gene_sourc e "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana";

ALINEAMIENTO

Formato SAM/BAM

“Sequence Alignment Map”, su versión binaria (comprimida) BAM. Ref

Formato para representar un alineamiento de NGS seq.

Programa asociado: Samtools

Alineamiento: Mapear las letras (bases) de dos o más secuencias, permitiendo algunos espaciadores (indels).

Variaciones dentro de un alineamiento: Indels Substituciones

El ADN se alinea de forma continua al genoma

Pero el mRNA puede tener un alineamiento dividido (spliced aligment) -> <- y forman un continuo

Contenido del formato SAM:

Header: líneas que empiezan con @ y dan información del alineamiento, la longitud de cada cromosoma, programa con el que se hizo, etc.

Alineamiento: una línea por cada alineamiento.

Contenido de las columnas:

Read id
FLAG
Chr
Start
Mapping quality
CIGAR (aligment)
MateChr
MateStart
MateDist
QuerySeq
QueryBaseQuals
AlignmentScore
Edot-distance-to-reference
Number-of-hits
Strand
Hit-index

Ejemplo:

Un alineamiento así:

Coor     12345678901234  5678901234567890123456789012345
ref      AGCATGTTAGATAA**GATAGCTGTGCTAGTAGGCAGTCAGCGCCAT
+r001/1        TTAGATAAAGGATA*CTG
+r002         aaaAGATAA*GGATA
+r003       gcctaAGCTAA
+r004                     ATAGCT..............TCAGC
-r003                            ttagctTAGGC
-r001/2                                        CAGCGGCAT

En SAM se codifica así:

@HD VN:1.5 SO:coordinate
@SQ SN:ref LN:45
r001   99 ref  7 30 8M2I4M1D3M = 37  39 TTAGATAAAGGATACTG *
r002    0 ref  9 30 3S6M1P1I4M *  0   0 AAAAGATAAGGATA    *
r003    0 ref  9 30 5S6M       *  0   0 GCCTAAGCTAA       * SA:Z:ref,29,-,6H5M,17,0;
r004    0 ref 16 30 6M14N5M    *  0   0 ATAGCTTCAGC       *
r003 2064 ref 29 17 6H5M       *  0   0 TAGGC             * SA:Z:ref,9,+,5S6M,30,1;
r001  147 ref 37 30 9M         =  7 -39 CAGCGGCAT         * NM:i:1

Preparacion del index 


## descarga de genoma de referencia

efetch -db nucleotide -id NC_001224.1 -format fasta > ngsexample/data/NC_001224.1.fasta

## creamos una carpeta para el indexado de genoma

mkdir -p ngsexample/data/index

## realizamos el indexado del genoma

bowtie2-build --threads 2 ngsexample/data/NC_001224.1.fasta ngsexample/data/index/mitocp

Alineamiento de secuencias 

lista="SRR21688982_mt_ SRR21688985_mt_ SRR21688981_mt_ SRR21688987_mt_"

for i in $lista
do
	bowtie2 --end-to-end -I 0 -X 1000 -p 3 \
	-x ngsexample/data/index/mitocp \
	-1 ngsexample/results/trimmed/${i}1_trimmed.fq.gz \
	-2 ngsexample/results/trimmed/${i}2_trimmed.fq.gz \
	-S ngsexample/results/map/${i}0.sam

	samtools view -u@ 2 ngsexample/results/map/${i}0.sam | \
	samtools sort -@ 2 -o ngsexample/results/map/${i}.sorted.bam -

	samtools index ngsexample/results/map/${i}.sorted.bam

done