Articles

Analisi dei dati CHIP-seq

Questo tutorial è stato ispirato dagli sforzi di Mo Heydarian e Mallory Freeberg. Strumenti higlighted qui sono stati avvolti da Björn Grüning, Marius van den Beek e altri membri IUC. Dave Bouvier e Martin Cech hanno contribuito alla messa a punto e alla distribuzione di strumenti sul server pubblico di Galaxy.

In questo tutorial ci sarà:

  • pre-processo di sequenziamento di legge
  • mappa di legge
  • post-process mapping di dati
  • valutare la qualità e la forza del ChIP-segnale
  • display copertura trame in una genome browser
  • chiamata ChIP con picchi MACS2
  • ispezionare ottenuto chiamate
  • cercare motivi di sequenza all’interno di detta cime
  • guarda distribuzione, arricchito regioni attraverso i geni.

Data

I set di dati per questo tutorial sono stati forniti da Shaun Mahony e sono stati generati nel laboratorio di Frank Pugh.

Reb1 ChIP-exo

Per questa analisi useremo set di dati ChIP-exo. Per questo esperimento l’immunoprecipitazione è stata eseguita con antobodies contro Reb1. Reb1 riconosce una sequenza specifica (TTACCCG) ed è coinvolto in molti aspetti della regolazione trascrizionale da tutte e tre le RNA polimerasi di lievito e promuove la formazione di regioni libere da nucleosomi (NFRs) (Hartley& Madhani:2009; Raisner:2005).

Sebbene si tratti di dati ChIP-exo, in questo tutorial lo analizzeremo come se si trattasse di CHIP-seq standard. Spiegheremo le peculiarità dell’analisi ChIP-exo in un tutorial dedicato.

Descrizione dei dati

Ci sono quattro set di dati:

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

Questi set di dati sono depositati in una libreria Galaxy (guarda il video su come importare i dati da una libreria):

Galaxy data library contenente le letture. Qui puoi vedere due repliche (R1 e R2). Si tratta di dati single-end. Carica i set di dati in una nuova cronologia selezionando tutti i set di dati e facendo clic sul pulsante to History. Assegnare un nome alla nuova cronologia e fare clic su Import (guarda questo video).

×

Creazione di una raccolta di dati/Dati di fine singolo

Caricamento

Dopo aver caricato i set di dati nella cronologia di Galaxy, combineremo tutti i set di dati in un’unica raccolta di dati. Ciò semplificherà l’elaborazione a valle dei dati. Il processo per la creazione di una raccolta per questo tutorial è è mostrato qui.

Mappatura e post-elaborazione

Mappatura

In questo caso particolare i dati sono di altissima qualità e non hanno bisogno di essere tagliati o post-elaborati in alcun modo prima della mappatura. Procederemo mappando tutti i dati rispetto all’ultima versione del genoma del lievito sacCer3:

Mappando tutti i dati contemporaneamente. Si noti che Select input type è impostato su Single fastq e selezionando folder () pulsante è possibile selezionare come intera raccolta di set di dati fastq. Importante: qui abbiamo anche impostare readgroups automaticamente attivando Set readgroups informazioni discesa a Set readgroups (SAM/BAM specification) e impostando tutti i pulsante Auto-assign a Yes.

L’esecuzione diBWA su una raccolta genererà un’altra raccolta di file BAM. Assegna un nome a questa raccolta mapped data (per informazioni su come rinominare una raccolta vedere questo video).
×

Rinominare una raccolta

Post-elaborazione

Per la post-elaborazione rimuoveremo tutte le letture mappate non univocamente. Questo può essere fatto semplicemente filtrando tutte le letture con una qualità di mappatura inferiore a 20 usando NGS: SAMtools → Filter SAM o BAM:

Filtering multi-mapped reads by restricting the data to reads with mapping quality above 20. Si noti che selezionando folder () pulsante è possibile selezionare come intera raccolta di set di dati BAM per filtrare in una sola volta.

L’esecuzione diFilter SAM or BAM su una raccolta genererà un’altra raccolta di file BAM. Assegna un nome a questa raccolta filtered data (per informazioni su come rinominare una raccolta vedere questo video).

Valutazione della qualità del chip

Dopo aver mappato e filtrato le letture è il momento di fare alcune inferenze su quanto siano buoni i dati sottostanti.

Correlazione tra campioni

In out experiment ci sono due repliche, ciascuna contenente set di dati di trattamento e input (controllo). La prima cosa che possiamo verificare è se i campioni sono correlati (in altre parole se i campioni di trattamento e controllo tra le due repliche contengono lo stesso tipo di segnale). Per fare ciò generiamo prima la matrice di conteggio delle letture usando NGS: DeepTools → multiBamSummary.

Esecuzione di multiBAMsummary su una raccolta di set di dati BAM (come prima è possibile selezionare la raccolta premendo il pulsante folder ()).

Questo strumento rompe genoma in bidoni di dimensioni fisse (10.000 bp nel nostro esempio) e calcola il numero di letture che rientrano in ogni bin. Ecco un frammento del suo output:

#'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

possiamo quindi alimentare questa matrice in NGS: DeepTools → plotCorrelation per generare una mappa di calore come questa:

A. Esecuzione diplotCorrelation sull’uscita dimultiBamSummary.

B. Heatmap di quattro campioni: Trattamenti (Rab1) e controlli (Input) sono ben correlati tra loro.

Qui possiamo vedere che ci sono buone correlazioni tra repliche (tra Reb1_R1 e Reb1_R2 e tra input_R1 e input_R2), mentre le correlazioni tra trattamenti (Reb1) e controlli (input) sono deboli. Questo è un buon segno che implica che c’è qualche segnale sui nostri dati.

Valutare la potenza del segnale

Come facciamo a dire è che abbiamo segnale proveniente da chip di arricchimento? Un modo per farlo è Signal Extraction Scaling (SES) proposto da Diaz:2012. SES funziona come segue. Supponiamo di avere due set di dati: ChIP e Input DNA. Dividiamo genoma in N finestre non sovrapposte (N = 10 nell’esempio qui sotto) e per ogni finestra calcolare il numero di letture. In questo modo finiamo con due liste: una lista di lettura conta per ChIP (lista ChIP) e l’altra per Input (lista Input):

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

Ordiniamo quindi l’elenco dei CHIP in ordine crescente e spostiamo gli elementi dall’Input-list in modo che corrispondano a questo ordine:

 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

Ora aggiungiamo altre due colonne a questo set di dati. Queste colonne mostreranno la percentuale di letture che si sommano a ciascuna riga per i dati di CHIP e input. Ad esempio, 0.044 sulla riga 3 è (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 

Nella matrice sopra una grande porzione di letture di CHIP (colonna 4) è concentrata nei pochi bidoni vicino al fondo. Questo non è il caso per le letture di input (colonna 5). Se tracciamo due ultime colonne di questa matrice otterremo una curva come questa:

SES plot per il nostro esempio di giocattolo. La maggior parte delle “letture” nell’esperimento ChIP sono concentrate negli ultimi tre bidoni.

DeepTools fornisce una bella spiegazione di come il successo di un esperimento di chip può essere giudicato in base ai grafici SES (chiamati anche impronte digitali):

DeepTools spiegazione dei grafici SES.

Quindi applichiamo questo ai nostri dati usando NGS: DeepTools → plotFingerprint:

A. Esecuzione di plotFingerprint su dati filtrati (15.).

B. SES impronta digitale di quattro campioni: Trattamenti (Rab1) mostrano caratteristica forma che indica di ChIP-segnale. Circa il 30% delle letture sono contenute in diversi % del genoma.

Generazione di set di dati bigWig per la visualizzazione

In questa sezione convertiremo i file BAM generati conbwa in formato bigWig che ci permetterà di visualizzare la distribuzione di copertura di lettura in tutto il genoma. Ci sarà anche “pre-caldo” un browser genoma per la visualizzazione di picchi chiameremo nella prossima sezione.

Generazione di set di dati bigWig

Useremo NGS: DeepTools → bamCoverage:

Eseguendo bamCoverage su una raccolta di set di dati BAM filtrati (come prima puoi selezionare la raccolta premendo il pulsante folder ()). Qui impostiamo la dimensione del cestino su 25. Successivamente impostiamo la dimensione effettiva del genoma suuser specified e inseriamo12000000 (dimensione approssimativa del genoma di Saccharomyces cerevisiae). Poiché questo strumento ha un’interfaccia particolarmente lunga, abbiamo tagliato sezioni importanti per creare questa immagine (vedi i riquadri sotto).

Infine impostiamo Extend reads alla dimensione media del frammento data a 150. Questo perché in questo particolare esperimento il DNA è stato selezionato per essere tra 120 e 170 bp per la preparazione della biblioteca.

L’esecuzione dibamCoverage su una raccolta di set di dati BAM genererà una raccolta di set di dati bigWig. Assegna un nome a questa raccolta coverage (per informazioni su come rinominare una raccolta vedere questo video).

Visualizzazione di tracce di copertura in un browser

×

Visualizzazione di più set di dati in IGV

Ora possiamo visualizzare set di dati bigWig generati nella sezione precedente in un browser genoma. C’è una varietà di browser disponibili. In questo tutorial useremo IGV Browser (questo video mostra come visualizzare più set di dati in IGV).

Raccolta di pezzi grossi prodotti da bamCoverage sopra. Si noti che in un set di dati espanso (Reb_R2) there is avisualizza il collegamento IGV`.

Facendo clic su questo collegamento in tutti e quattro i set di dati (è necessario espandere ciascun set di dati facendo clic su di esso. Questo esporrà IGV link) e concentrandosi browser su MPH1 (YIR002C) gene produrrà la seguente immagine:

Copertura distribuzione all’interno IGV. Qui le repliche di chip sono colorate in arancione e i controlli sono blu. Tutte e quattro le tracce sono state impostate sul valore massimo di 70.

Chiamare picchi

Mentre i picchi mostrati nello screenshot del browser sopra sono abbastanza chiari e coerenti tra le due repliche, guardando l’intero genoma nel browser non è certo un modo sostenibile per identificare tutti i picchi. Esistono diversi modi per identificare gli eventi di legame a livello genomico. Sono riassunti nella figura seguente:

Schema di tre metodi di rilevamento degli eventi di binding CHIP-seq. I metodi di ricerca del picco in genere spostano le posizioni dei tag CHIP-seq in una direzione 3′ della metà della lunghezza prevista del frammento o estendono la lunghezza del tag in una direzione 3′ per essere uguale alla lunghezza prevista del frammento. I tag dei trefoli opposti vengono uniti per costruire paesaggi di densità dei tag non marchiati e le posizioni degli eventi vincolanti sono previste dalle posizioni con la massima copertura dei tag all’interno di ciascuna regione che contiene un significativo arricchimento dei tag chip-seq (cioè il vertice del picco). Analisi dei metodi di accoppiamento di picco tra la mappatura delle letture a + e-trefoli. Diamo un’occhiata a questi risultati:

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. Useremo una media di questi valori, 30, come --extsize parametro per la chiamata di picchi di utilizzo di NGS: Picco di chiamata → MACS2 callpeak:

Chiamare con picchi MACS2 dati raccolti. Qui scegliamo più ingressi premendo il pulsante e selezionando entrambi i set di dati di chip nel file di trattamento ChIP-Seq e entrambi i set di dati di input DNA nel file di controllo ChIP-Seq. Selezioniamo quindiSaccharomyces cerevisiae genoma come dimensione effettiva del genoma. MACS2s interfaccia è lunga e abbiamo diviso in più pezzi in questa figura. Vedi anche la sezione inferiore – è importante!

In questa parte inferiore del MACS2 interfaccia, impostare la costruzione di un modello di Do not build the shifting model (abbiamo già fatto con preductd nel passaggio precedente) e *Impostare la dimensione dell’estensione 30 (il numero si stima che nel passaggio precedente). Infine, chiederemo solo MACS2 per produrre due uscite: Peak summits e quello che ha prodotto per impostazione predefinita, che contiene le coordinate di picco.

Se si impostano i parametri come è stato mostrato sopra MACS2 produrrà due uscite (se ne ha prodotte di più basta trovare quelle chiamate narrow peaks e summits). Facciamo clic sull’icona della matita () adiacente ai set di datisummits enarrow peak e rinominiamo quindi come mostrato di seguito:

MACS2 output with summits and narrow peak datasets renamed. .

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

Chiamare picchi conMACS2 su R1 Ad eccezione di selezionare solo set di dati R1, tutti gli altri parametri devono essere impostati come nella figura precedente.

Ora fare questo da soli:

  • rinomina risultante set R1 summits e R1 peaks
  • esegui MACS2 esegui su Replicare 2
  • rinomina risultante summits e narrow peak set R2 summits e R2 peaks.

Alla fine dovresti avere qualcosa del genere:

Tutti i MACS2 dataset. Dopo aver eseguitoMACS2 tre volte avremmo dovuto interrogare i set di dati R1 e R2.

Ispezione dei picchi

Cosa c’è nell’output?

Guardando i dati MACS2 abbiamo ottenuto i seguenti numeri di picchi:

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. Start
  3. Fine
  4. Iterativo id offerto da MACS2
  5. Integer punteggio per display
  6. Strand (irrilevante in questo caso)
  7. Fold-change (fold arricchimento per questo picco di vertice contro casuali con distribuzione di Poisson con locale lambda)
  8. -log10P-valore (ad esempio, 17.68 1 x 10-17)
  9. -log10Q-valore da Benjamini–Hochberg–Yekutieli procedura
  10. Relativa vertice posizione a picco start

Come molte cime sono comuni tra repliche?

Per vedere quanti picchi sono comuni tra i set di dati in pool e le due repliche useremo Operate on Genomic Intervals → Join tool twice.

Galaxy main ha due strumenti chiamati Join. Non confonderli! Qui stiamo usando quello da Operare sulla sezione Intervalli genomici.

Prima di tutto ci uniremo Peaks pooled con Peaks R1:

Entrare in Pool e R1 risultati con Join strumento. Si noti che, poiché abbiamo rinominato i set di dati, ora sono facilmente selezionabili.

Successivamente uniremo il risultato dell’operazione precedente con Peaks R2:

Unire Pooled/R1 con i risultati R2 con Join strumento. Si noti che, poiché abbiamo rinominato i set di dati, ora sono facilmente selezionabili.

Questo si traduce in 723 regioni sono condivisi tra intervistati, R1, e picchi R2. Chiamiamo questo set di alta fiducia. Prima di poterlo utilizzare, tuttavia, tagliamo solo le colonne pertinenti. Poiché abbiamo prodotto questo set di dati unendo altri tre set di dati, è tre volte più ampio (30 colonne). Per tagliare queste prime tre colonne possiamo usare lo strumento Manipolazione testo → Taglia colonne:

Tagliare colonne da Join output.

Rinominare l’ultimo set di dati comeHigh confidence set. Questo lo renderà facile da trovare mentre continuiamo.
Utilizzando Cut columns lo strumento produce un set di dati di tipo tabulare. Tuttavia, tagliando le prime dieci colonne abbiamo creato un set di dati in formato LETTO. Quindi dobbiamo far sapere a Galaxy che ripristinando i metadati come mostrato di seguito.

Quindi dobbiamo assicurarci che l’output di Cut columns abbia il tipo BED. Per fare ciò modificheremo i metadati come mostrato di seguito:

Impostando i metadati sul tipo di dati BED. Fare clic sull’icona a forma di matita() adiacente al set di dati e scegliere la scheda Tipo di dati. Lì sarai in grado di impostarlo su BED.

vediamo tutto nel browser

visualizzare Unite sulle vette Stretti picchi e Vette prodotto da MACS2 in IGV cliccando su display with IGV local link adiacente al Peaks pooled e High confidence set dataset (si dovrebbe già avere il browser aperto):

Una panoramica IGV. Qui puoi vedere i set di dati originali bigWig insieme ai picchi previsti.

Quali motivi di sequenza si trovano all’interno dei picchi

In questo esperimento gli anticorpi contro la proteina Reb1 sono stati utilizzati per l’immunoprecipitaion. Il sito di riconoscimento per Reb1 è TTACCCG (Badis:2008 e Harbison:2004). Per scoprire quali motivi di sequenza si trovano all’interno dei nostri picchi, dobbiamo prima convertire le coordinate in sequenze sottostanti. Questo viene fatto usando Fetch Allineamenti / Sequenze → Estrai lo strumento DNA genomico:

Estraendo il DNA genomico corrispondente ai picchi di CHIP-seq. Qui usiamoMerged peaks dataset generato pochi passaggi prima.

Successivamente, dobbiamo assicurarci che tutte le sequenze siano sufficientemente lunghe per trovare modelli. MEME, gli strumenti che useremo per trovare motivi, sequenze necessarie per essere almeno 8 nucleotidi lunghi. Quindi rimuoveremo sequenze brevi usando FASTA manipulation → Filter sequences by length tool:

Filtering FASTA by length. Qui stiamo rimuovendo tutte le sequenze in corto di 8 nucleotidi.

Ora possiamo eseguire Motif Tools → MEME:

Esecuzione di MEME su sequenze FASTA filtrate in lunghezza dal passaggio precedente. Si noti che la configurazione delle opzioni è impostata su Advancede Check reverse complement è impostata su Yes.

MEME genera un numero di uscite. Il più interessante è il report HTML. Mostra che 620 regioni contengonoTTACCCG motif:

MEME Motif trovato in 620 sequenze corrispondenti alle regioni di picco comuni.

Riassumendo CHIP arricchimento del segnale in tutti i geni

Quanti geni contengono regioni a monte arricchiti in tag ChIP. Questo è spesso rappresentato come una heatmap:

Esempio di Heatmap dalla documentazione di DeepTools.

Per generare la heatmap dobbiamo prima produrre set di dati normalizzati per i due replicati che abbiamo. Questo viene fatto usando NGS: DeepTools → bamCompare tool:

In esecuzione bamCompare sulla replica 1. Qui impostiamo il metodo da utilizzare per ridimensionare il campione più grande al più piccolo a SES (anche se potresti voler provare anche altri metodi. SES è stato brevemente discusso sopra.

Eseguire la stessa analisi su Replicare 2 set di dati e rinominare i due elementi risultanti comeR1 normalized eR2 normalized.

Poiché vogliamo tracciare l’arricchimento attorno ai geni, dobbiamo scaricare l’annotazione genica. Useremo Get Data → UCSC Main per questo:

Ottenere dati da UCSC. Qui assicurati di selezionare assembly chiamatosacCer3 e stai scegliendoSGD Genes. Facendo clic su Ottieni output verrà visualizzata la schermata successiva mostrata di seguito.

Qui basta fare clic su Invia query a Galaxy.

Successivamente, per preparare i dati necessari per disegnare la heatmap useremo NGS: DeepTools → computeMatrix utility:

Computing matrix – i dati da cui verrà costruita la heatmap. Qui entrambi i set di dati normalizzati sono selezionati all’interno della casella file Score, i geni di lievito che abbiamo appena scaricato da UCSC sono scelti come regioni da tracciare. “il punto di riferimento è impostato come opzione principale computeMatrix e, infine, le distanze a monte e a valle sono impostate su 2.000 bp. Ovviamente siete invitati a giocare con questi parametri.

Infine, possiamo visualizzare la heatmap usando NGS: DeepTools → plotHeatmap tool:

Disegnando heatmap conplotHeatmap tool.

L’immagine risultante mostra che una frazione significativa di 6,692 geni presenti nell’annotazione dei dati abbiamo utilizzato contenere Reb1 siti di legame all’interno delle loro regioni a monte:

Heatmap che mostra la distribuzione di Reb1 siti di legame tra regioni a monte di 6,692 lievito geni.

Non abbiamo falsificato questo

Questa intera analisi è disponibile come cronologia della galassia qui. Importarlo e giocare con esso.

E così va…

Speriamo che questo tutorial ti abbia dato il gusto per ciò che è possibile. Ci sono più strumenti là fuori così esperimento! Se le cose non funzionano – si lamentano utilizzando Open Chat pulsante qui sotto o il nostro forum di supporto.