Articles

Análisis de datos de CHIP-seq

Este tutorial se inspiró en los esfuerzos de Mo Heydarian y Mallory Freeberg. Björn Grüning, Marius van den Beek y otros miembros de la IUC han envuelto las herramientas aquí iluminadas. Dave Bouvier y Martin Cech ayudaron a ajustar e implementar herramientas en el servidor público de Galaxy.

En este tutorial vamos a:

  • lecturas de secuenciación previa al proceso
  • lecturas de mapa
  • datos mapeados posteriores al proceso
  • evaluar la calidad y la potencia de la señal de ChIP
  • mostrar gráficas de cobertura en un navegador de genoma
  • picos de chip de llamada con MACS2
  • inspeccionar las llamadas obtenidas
  • buscar motivos de secuencia dentro de li>
  • observe la distribución de las regiones enriquecidas entre los genes.

Datos

Los conjuntos de datos para este tutorial fueron proporcionados por Shaun Mahony y se generaron en el laboratorio de Frank Pugh.

Reb1 ChIP-exo

Para este análisis utilizaremos conjuntos de datos ChIP-exo. Para este experimento se realizó inmunoprecipitación con antobodios contra Reb1. Reb1 reconoce una secuencia específica (TTACCCG) y participa en muchos aspectos de la regulación transcripcional de las tres ARN polimerasas de levadura y promueve la formación de regiones libres de nucleosomas (NFR) (Hartley & Madhani:2009; Raisner:2005).

Aunque se trata de datos ChIP-exo, en este tutorial lo analizaremos como si fuera ChIP-seq estándar. Explicaremos las peculiaridades del análisis ChIP-exo en un tutorial dedicado.

Descripción de los datos

Hay cuatro conjuntos de datos:

Dataset Description
Reb1_R1 ChIP experiment, Replicate 1
Input_R1 Input DNA, Replicate 1
Reb1_R2 ChIP experiment, Replicate 2
Input_R2 Input DNA, Replicate 2

Data location

×

Importing from History

Estos conjuntos de datos se depositan en una biblioteca Galaxy (vea el video sobre cómo importar datos de una biblioteca):

Biblioteca de datos Galaxy que contiene las lecturas. Aquí puedes ver dos repeticiones (R1 y R2). Estos son datos de un solo extremo. Cargue conjuntos de datos en un nuevo historial seleccionando todos los conjuntos de datos y haciendo clic en el botón to History. Asigne un nombre al nuevo historial y haga clic en Import (vea este video).

×

Crear una colección de conjuntos de datos | datos de un solo extremo

Cargar

Después de cargar conjuntos de datos en el historial de Galaxy, combinaremos todos los conjuntos de datos en una única colección de conjuntos de datos. Esto simplificará el procesamiento posterior de los datos. El proceso para crear una colección para este tutorial se muestra aquí.

Asignación y postprocesamiento

Asignación

En este caso particular, los datos son de muy alta calidad y no es necesario recortarlos ni procesarlos de ninguna manera antes de la asignación. Procederemos mapeando todos los datos con la última versión del genoma de levadura sacCer3:

Mapeando todos los datos a la vez. Tenga en cuenta que Seleccionar tipo de entrada se establece en Single fastq y, al seleccionar el botón folder (), puede seleccionar como colección completa de conjuntos de datos fastq. Importante: aquí también podemos establecer readgroups automáticamente activando Conjunto readgroups desplegable información a Set readgroups (SAM/BAM specification) y configuración de todos los Auto-asignar botón Yes.

Ejecutar BWA en una colección generará otra colección de BAM archivos. Nombre de esta colección mapped data(para obtener ayuda sobre cómo cambiar el nombre de una colección, consulte este video).
×

cambiar el nombre de una colección

Post-procesamiento

Para el post-procesamiento vamos a eliminar todos los no-exclusiva de lecturas asignadas. Esto se puede hacer simplemente filtrando todas las lecturas con una calidad de asignación inferior a 20 utilizando NGS: SAMtools → Filtrar SAM o BAM:

Filtrar lecturas de asignaciones múltiples restringiendo los datos a lecturas con una calidad de asignación superior a 20. Tenga en cuenta que al seleccionar el botón folder () puede seleccionar como colección completa de conjuntos de datos BAM para filtrar a la vez.

Ejecutar Filter SAM or BAM en una colección generará otra colección de BAM archivos. Nombre de esta colección filtered data(para obtener ayuda sobre cómo cambiar el nombre de una colección, consulte este video).

Evaluación de la calidad del chip

Después de mapear y filtrar las lecturas, es hora de hacer algunas inferencias sobre cuán buenos son los datos subyacentes.

Correlación entre muestras

En el experimento out hay dos réplicas, cada una de las cuales contiene conjuntos de datos de tratamiento y de entrada (control). Lo primero que podemos comprobar es si las muestras están correlacionadas (en otras palabras, si las muestras de tratamiento y control de las dos réplicas contienen este mismo tipo de señal). Para hacer esto, primero generamos la matriz de conteo de lectura usando NGS: DeepTools → multiBamSummary.

Ejecutar multiBAMsummary en una colección de conjuntos de datos BAM (como antes, puede seleccionar la colección pulsando el botón folder ()).

Esta herramienta divide el genoma en contenedores de tamaño fijo (10,000 bp en nuestro ejemplo) y calcula el número de lecturas que caen dentro de cada contenedor. Aquí hay un fragmento de su salida:

#'chr' 'start' 'end' 'Reb1_R1' 'Input_R1' 'Input_R1' 'Reb1_R2'chrVI 0 1000 19.0 41.0 3.0 6.0chrVI 1000 2000 29.0 30.0 13.0 5.0chrVI 2000 3000 0.0 0.0 0.0 0.0chrVI 3000 4000 0.0 2.0 0.0 0.0chrVI 4000 5000 7447.0 139.0 7.0 2645.0

luego podemos alimentar esta matriz en NGS: DeepTools → plotCorrelation para generar un mapa de calor como este:

A. Ejecuta plotCorrelation en la salida de multiBamSummary.

B. Mapa de calor de cuatro muestras: Tratamientos (Rab1) y controles (Entrada) están bien correlacionados entre sí.

Aquí podemos ver que hay buenas correlaciones entre réplicas (entre Reb1_R1 y Reb1_R2, y entre input_R1 y input_R2), mientras que las correlaciones entre tratamientos (Reb1) y controles (input) son débiles. Esta es una buena señal que implica que hay alguna señal en nuestros datos.

Evaluar la intensidad de la señal

¿Cómo sabemos si tenemos señal procedente del enriquecimiento de chips? Una forma de hacerlo es el Escalado de Extracción de señales (SES) propuesto por Díaz:2012. SES funciona de la siguiente manera. Supongamos que tenemos dos conjuntos de datos: ChIP y ADN de entrada. Dividimos el genoma en N ventanas no superpuestas (N = 10 en el ejemplo siguiente) y para cada ventana calculamos el número de lecturas. De esta manera, terminamos con dos listas: una lista de lectura cuenta para ChIP (Lista de chips) y la otra para Entrada (Lista de entrada):

Window ChIP-count Input-count------------------------------- 1 3 3 2 4 3 3 2 1 4 1 3 5 3 3 6 27 2 7 18 3 8 2 2 9 45 310 8 3

Luego ordenamos la lista de chips en orden ascendente y movemos elementos de la lista de entrada para que coincidan con este orden:

 Window ChIP-count Input-count-------------------------------4 1 33 2 18 2 21 3 35 3 32 4 310 8 37 18 36 27 29 45 3

Ahora agreguemos otras dos columnas a este conjunto de datos. Estas columnas mostrarán el porcentaje de lecturas sumando a cada fila para datos de ChIP y de entrada. Por ejemplo, 0.044 en la fila 3 es (1 + 2 + 2)/113 = 0.044.

 1 2 3 4 5------------------------------- 4 1 3 0.008 0.115 3 2 1 0.026 0.153 8 2 2 0.044 0.230 1 3 3 0.070 0.346 5 3 3 0.097 0.230 2 4 3 0.132 0.57610 8 3 0.203 0.692 7 18 3 0.362 0.807 6 27 2 0.601 0.884 9 45 3 1.000 1.000------------------------ 113 26 Where:1 = Window, 2 = read count in ChIP, 3 = read count in Input4 = % or read to this point in ChIP 5 = % of read to this point 

En la matriz de arriba, una gran parte de las lecturas de chips (columna 4) se concentra en los pocos contenedores cercanos a la parte inferior. Este no es el caso de las lecturas de entrada (columna 5). Si trazamos dos últimas columnas de esta matriz, obtendremos una curva como esta:

parcela SES para nuestro ejemplo de juguete. La mayoría de las» lecturas » en el experimento de chips se concentran en los últimos tres contenedores.

DeepTools proporcionar una buena explicación de cómo el éxito de un ChIP experimento puede ser juzgado en función de SES (también llamado de huellas dactilares) parcelas:

DeepTools explicación de SES parcelas.

Así que apliquemos esto a nuestros propios datos usando NGS: DeepTools → plotFingerprint:

A. Ejecutando plotFingerprint en los datos filtrados (15.).

B. SES huella dactilar de cuatro muestras: Los tratamientos (Rab1) muestran una forma característica que indica la señal de ChIP. Aproximadamente el 30% de las lecturas están contenidas en varios % del genoma.

Generar conjuntos de datos bigWig para visualización

En esta sección convertiremos los archivos BAM generados con bwa en formato bigWig que nos permitirá ver la distribución de cobertura de lectura en todo el genoma. También «precalentaremos» un navegador genómico para mostrar los picos que llamaremos en la siguiente sección.

Generar conjuntos de datos bigWig

Usaremos NGS: DeepTools → bamCoverage:

Ejecutar bamCoverage en una colección de conjuntos de datos BAM filtrados (como antes, puede seleccionar la colección pulsando el botón folder ()). Aquí establecemos el tamaño de la bandeja en 25. A continuación, establecemos el tamaño efectivo del genoma en user specified e introducimos 12000000 (tamaño aproximado del genoma de Saccharomyces cerevisiae). Debido a que esta herramienta tiene una interfaz particularmente larga, recortamos secciones importantes para hacer esta imagen (vea los paneles a continuación).

Finalmente establecemos lecturas extendidas al tamaño promedio de fragmento dado a 150. Esto se debe a que en este experimento en particular se seleccionó el tamaño del ADN entre 120 y 170 pb para la preparación de la biblioteca.

Ejecutar bamCoverage en una colección de conjuntos de datos BAM generará una colección de conjuntos de datos bigWig. Nombre de esta colección coverage(para obtener ayuda sobre cómo cambiar el nombre de una colección, consulte este video).

Mostrar pistas de cobertura en un navegador

×

Mostrar múltiples conjuntos de datos en IGV

Ahora podemos mostrar conjuntos de datos bigWig generados en la sección anterior en un navegador de genoma. Hay una variedad de navegadores disponibles. En este tutorial usaremos el navegador IGV (este video muestra cómo mostrar múltiples conjuntos de datos en IGV).

la Colección de peces gordos producido por bamCoverage arriba. Tenga en cuenta que en el conjunto de datos expandido (Reb_R2) there is ase muestra en el enlace IGV.

Haciendo clic en este enlace en los cuatro conjuntos de datos (tendrá que expandir cada conjunto de datos haciendo clic en él. Esto expondrá el enlace IGV) y el navegador enfocado en el gen MPH1 (YIR002C) producirá la siguiente imagen:

Distribución de cobertura dentro de IGV. Aquí las réplicas de chips son de color naranja y los controles son azules. Las cuatro pistas se establecieron en el valor máximo de 70.

Llamar a picos

Mientras que los picos que se muestran en la captura de pantalla del navegador anterior son bastante claros y consistentes en las dos réplicas, observar todo el genoma en el navegador no es una forma sostenible de identificar todos los picos. Hay varias maneras de identificar eventos de unión en todo el genoma. Se resumen en la figura que figura a continuación:

Esquema de tres métodos de detección de eventos de unión ChIP-seq. Los métodos de búsqueda de picos normalmente desplazan las ubicaciones de la etiqueta ChIP-seq en una dirección de 3′ a la mitad de la longitud esperada del fragmento, o extienden la longitud de la etiqueta en una dirección de 3′ para que sea igual a la longitud esperada del fragmento. Las etiquetas de cadenas opuestas se fusionan para construir paisajes de densidad de etiquetas sin marca, y las ubicaciones de eventos de enlace se predicen a partir de las ubicaciones con la máxima cobertura de etiquetas dentro de cada región que contiene un enriquecimiento significativo de etiquetas ChIP-seq (es decir, la cumbre de pico). Análisis de métodos de emparejamiento de picos entre la asignación de lecturas a hebras + y -. Veamos estos resultados:

Replicate 1 Replicate 2
Peak model and lag between strands.

In the case of these data peaks are very sharp and have narrow gap between them: 27 and 33 bp for replicate 1 and 2, respectively. Usaremos un promedio de estos valores, 30, como --extsize parámetro para llamar a picos usando NGS: Llamada a picos → pico de MACS2:

Llamar a picos con MACS2 en datos agrupados. Aquí elegimos múltiples entradas presionando el botón y seleccionando ambos conjuntos de datos de ChIP en el Archivo de Tratamiento ChIP-Seq y ambos conjuntos de datos de ADN de entrada en el Archivo de Control ChIP-Seq. Luego seleccionamosSaccharomyces cerevisiae genoma como el tamaño efectivo del genoma. MACS2 la interfaz es larga y la dividimos en varias partes en esta figura. Vea también la sección inferior , ¡ es importante!

En esta parte inferior de MACS2 la interfaz establece el modelo de compilación en Do not build the shifting model (ya lo hemos hecho con preductd en el paso anterior) y *Establezca el tamaño de la extensión en 30 (el número que estimamos en el paso anterior). Finalmente, solo le pediremos a MACS2 que produzca dos salidas: Peak summits y el que produjo por defecto, que contiene coordenadas de pico.

Si establece los parámetros como se muestra anteriormente MACS2 producirá dos salidas (si produjo más, busque las llamadas narrow peaks y summits). Vamos a hacer clic en el icono de lápiz() adyacente a summits y narrow peak conjuntos de datos y renombrarlos como se muestra a continuación:

MACS2 output with summits and narrow peak datasets renamed. .

Next, we will run MACS2 on BAM datasets for Replicate 1 only:

Picos de llamada con MACS2 en R1 Con la excepción de seleccionar solo conjuntos de datos R1, todos los demás parámetros deben establecerse como en la figura anterior.

Ahora haga esto usted mismo:

  • cambie el nombre de los conjuntos de datos resultantes como R1 summits y R1 peaks
  • run MACS2 run on Replicate 2
  • rename resulting summits y narrow peak conjuntos de datos como R2 summits y R2 peaks.

al final, usted debe tener algo como esto:

Todos los MACS2 conjuntos de datos. Después de ejecutar MACS2 tres veces, deberíamos haber encuestado, R1 y R2 conjuntos de datos.

Inspección de picos

¿Qué hay en la salida?

Mirando los datos de MACS2, hemos obtenido los siguientes números de picos:

Pooled Replicate 1 Replicate 2
974 955 784

Peaks data is generated in the following format:

 1 2 3 4 5 6 7 8 9 10--------------------------------------------------------------------chrI 35 491 MACS2_peak_1 176 . 10.35332 21.51081 17.68957 101chrI 87135 87212 MACS2_peak_2 127 . 7.71763 15.89278 12.78060 26chrI 92612 92793 MACS2_peak_3 153 . 9.22373 18.72748 15.31966 49chrI 119739 119782 MACS2_peak_4 78 . 6.08885 10.52482 7.82302 25

where columns are:

  1. Cromosoma
  2. Inicio
  3. Fin
  4. Id iterativo dado porMACS2
  5. Puntuación entera para visualización
  6. Strand (irrelevante en este caso)
  7. Cambio de plegado (enriquecimiento de plegado para esta cumbre de pico contra distribución aleatoria de Poisson con lambda local)
  8. -log10P-valor (por ejemplo, 17.68 es 1 x 10-17)
  9. -log10Q-valor del procedimiento Benjamini-Hochberg–Yekutieli
  10. Posición relativa de la cumbre al inicio del pico

¿Cuántos picos son comunes entre las réplicas?

Para ver cuántos picos son comunes entre los conjuntos de datos agrupados y las dos réplicas, usaremos la herramienta Operar en Intervalos Genómicos → Unir dos veces.

Galaxy main tiene dos herramientas llamadas Unir. ¡No los confundas! Aquí estamos usando el de la sección Operar en Intervalos Genómicos.

en Primer lugar vamos a unir Peaks pooledPeaks R1:

Unirse en común, y R1 resultados con la etiqueta Join herramienta. Tenga en cuenta que, debido a que cambiamos el nombre de los conjuntos de datos, ahora se pueden seleccionar fácilmente.

A continuación uniremos el resultado de la operación anterior con Peaks R2:

Uniendo los resultados de Pooled/R1 con R2 con Join herramienta. Tenga en cuenta que, debido a que cambiamos el nombre de los conjuntos de datos, ahora se pueden seleccionar fácilmente.

Estos resultados en 723 regiones se comparten entre los picos encuestados, R1 y R2. Llamemos a este conjunto de Alta confianza. Sin embargo, antes de que podamos usarlo, recortemos solo las columnas relevantes. Dado que hemos producido este conjunto de datos uniendo otros tres conjuntos de datos, es tres veces más ancho (30 columnas). Para cortar estas tres primeras columnas podemos usar la herramienta Manipulación de texto → Cortar columnas:

Cortar columnas de salida Join.

Cambie el nombre del último conjunto de datos como High confidence set. Esto hará que sea fácil de encontrar a medida que continuamos.
Usando Cut columns la herramienta produce un conjunto de datos de tipo tabular. Sin embargo, al cortar las primeras diez columnas, hemos creado un conjunto de datos en formato de CAMA. Por lo tanto, debemos informar a Galaxy de eso restableciendo los metadatos como se muestra a continuación.

A continuación, debemos asegurarnos de que la salida de Cut columns tool tenga el tipo BED. Para hacer esto, editaremos sus metadatos como se muestra a continuación:

Estableciendo los metadatos en el tipo de datos BED. Haga clic en el icono de lápiz() junto al conjunto de datos y elija la pestaña Tipo de datos. Allí podrá configurarlo en BED.

Veamos todo en el navegador

Visualicemos picos fusionados, así como picos y cumbres estrechos producidos por MACS2 en IGV haciendo clic en display with IGV local enlaces adyacentes a Peaks pooled y High confidence set conjuntos de datos (ya debería tener el navegador abierto):

Una descripción general en IGV. Aquí puede ver los conjuntos de datos originales de bigWig junto con los picos pronosticados.

Qué motivos de secuencia se encuentran dentro de los picos

En este experimento, se han utilizado anticuerpos contra la proteína Reb1 para la inmunoprecipitación. El sitio de reconocimiento para Reb1 es TTACCCG(Badis: 2008 y Harbison: 2004). Para averiguar qué motivos de secuencia se encuentran dentro de nuestros picos, primero necesitamos convertir las coordenadas en secuencias subyacentes. Esto se hace utilizando la herramienta Obtener alineaciones / Secuencias → Extraer ADN genómico:

Extrayendo ADN genómico correspondiente a picos ChIP-seq. Aquí usamos Merged peaks conjunto de datos generado unos pasos antes.

A continuación, debemos asegurarnos de que todas las secuencias sean lo suficientemente largas para encontrar patrones. MEME, las herramientas que usaremos para encontrar motivos, requerían secuencias de al menos 8 nucleótidos de longitud. Por lo tanto, eliminaremos secuencias cortas utilizando la herramienta Manipulación de FASTA → Filtrar secuencias por longitud:

Filtrar FASTA por longitud. Aquí estamos eliminando todas las secuencias con cortocircuito de 8 nucleótidos.

Ahora podemos ejecutar Herramientas de motivos → MEME:

Ejecutar MEME en secuencias FASTA filtradas de longitud del paso anterior. Tenga en cuenta que la configuración de opciones se establece en Advanced y Comprobar complemento inverso se establece en Yes.

MEME genera un número de salidas. El más interesante es el informe HTML. Muestra que 620 regiones contienen TTACCCG motivo:

Motivo de MEME encontrado en 620 secuencias correspondientes a regiones de pico comunes.

Resumiendo el enriquecimiento de la señal de ChIP en todos los genes

Cuántos genes contienen regiones aguas arriba enriquecidas en etiquetas de ChIP. Esto a menudo se representa como un mapa de calor:

Ejemplo de mapa de calor de la documentación de DeepTools.

Para generar el mapa de calor primero debemos producir conjuntos de datos normalizados para los dos replicados que tenemos. Esto se hace mediante NGS: DeepTools → bamCompare herramienta:

Ejecutar bamCompare en replicar 1. Aquí establecemos el método a usar para escalar la muestra más grande a la más pequeña a SES (aunque es posible que desee probar otros métodos también. El SES se examinó brevemente más arriba.

Realice el mismo análisis en Replicar 2 conjuntos de datos y cambie el nombre de los dos elementos resultantes a R1 normalized y R2 normalized.

Debido a que queremos trazar el enriquecimiento alrededor de los genes, necesitamos descargar la anotación de genes. Usaremos Get Data → UCSC Main para esto:

Obtener datos de UCSC. Aquí asegúrese de seleccionar ensamblado llamado sacCer3 y está eligiendo SGD Genes. Al hacer clic en obtener salida, se mostrará la siguiente pantalla que se muestra a continuación.

Aquí simplemente haga clic en Enviar consulta a Galaxy.

A continuación, para preparar los datos necesarios para dibujar el mapa de calor, usaremos NGS: DeepTools → Utilidad computeMatrix:

Matriz de computación: los datos a partir de los cuales se construirá el mapa de calor. Aquí se seleccionan ambos conjuntos de datos normalizados dentro del cuadro de archivos de puntuación, los genes de levadura que acabamos de descargar de UCSC se eligen como Regiones para trazar. «el punto de referencia se establece como opción principal de computeMatrix y, finalmente, las distancias aguas arriba y aguas abajo se establecen en 2.000 bp. Obviamente, le invitamos a jugar con estos parámetros.

Finalmente, podemos visualizar el mapa de calor utilizando NGS: DeepTools → Herramienta plotHeatmap:

Dibujar mapa de calor con plotHeatmap herramienta.

La imagen resultante muestra que una fracción significativa de 6.692 genes presentes en los datos de anotación que hemos utilizado contienen sitios de unión a Reb1 dentro de sus regiones ascendentes:

Mapa de calor que muestra la distribución de sitios de unión a Reb1 a través de regiones ascendentes de 6.692 genes de levadura.

No hemos falsificado este

Este análisis completo está disponible como historial de galaxias aquí. Impórtalo y juega con él.

Y así sigue…

Esperemos que este tutorial le haya dado el gusto de lo que es posible. Hay más herramientas por ahí, ¡así que experimenta! Si las cosas no funcionan, quejese usando el botón Open Chat a continuación o nuestro foro de soporte.