Analyses bioinformatiques et statistiques de données de ChIP-seq sous UNIX et R - Journée 2

Transcript

1. Analyses bioinformatiques et statistiques de données de ChIP-seq sous UNIX et R Agrocampus Ouest - Avril 2017 Pierre-François ROUX Département « Biologie cellulaire & infection » Unité « Organisation nucléaire & oncogenèse » Institut Pasteur – Paris pierre-francois.roux@pasteur.fr CELL BIOLOGY & INFECTION 1 Journée 2
2. Bioconductor, c’est quoi ? • Bibliothèque R regroupant plus de 800 packages pour l’analyse de données génomiques • Documentation très complète pour chaque package (vignette) • Mise à jour tous les 6 mois • Mise à disposition de matériel de cours • Exemples de workflows • Chaine F1000 pour la publication de workflow • Communauté d’utilisateurs / développeurs très vaste et active Bioconductor Analyse de données NGS sous R 2
3. Bioconductor Analyse de données de ChIP-seq sous R 3 Ressources Bioconductor pour l’analyse de données de ChIP-seq • GenomicRanges : Grammaire et fonctions pour la manipulation de données d’intervalles • Rsamtools : Grammaire et fonction pour la manipulation des BAM • DiffBind : Analyse différentielle de données de ChIP-seq • rtracklayer : Imports / exports de données de séquences et d’annotation • DESeq : Analyse différentielle pour données de comptage • edgeR : Analyse différentielle pour données de comptage • chipseq : Outils pour l’analyse de données de ChIP-seq • ChIPpeakAnno : Annotation et enrichissement fonctionnel pour le ChIP-seq • MotifDb : Collection de motifs et PWMs • motifStack : Outils graphique pour l’affichage des PWMs • PWMEnrich : Outils pour les tests d’enrichissement en motifs • …
4. Bioconductor Analyse de données de ChIP-seq sous R 4 Ressources Bioconductor pour l’analyse de données de ChIP-seq (peak-calling) • BayesPeak : Chaines de Markov cachées (HMM) et algorithme Bayésien • PICS : Inférence probabilistique • MOSAiCS : Model-based analysis of ChIP-Seq data • iSeq : Modèle bayésien Ising caché • ChIPseqR : Positionnement de nucléosome • CSAR : Tests basés sur la loi de Poisson • SPP : peak-caller développé dans le cardre du projet ENCODE • NarrowPeaks : Analyse en composante principale fonctionnelle (FPCA)
5. Aperçu d’un pipeline Bioconductor de traitement de données ChIP-seq Analyse de données de ChIP-seq sous R Données brutes .fastq Données trimmées .fastq QuasR QuasR Rsamtools .bam .bam Prétraitement Données alignées Données filtrées QuasR • Qualité des reads • Contenu en GC • Taux de duplication • Séquences sur-représentées QuasR • Taux d’alignement csaw • Estimation de la taille des fragments
6. Aperçu d’un pipeline Bioconductor de traitement de données ChIP-seq Analyse de données de ChIP-seq sous R Données brutes .fastq Données trimmées .fastq QuasR QuasR Rsamtools .bam .bam Prétraitement Données filtrées .bam Pics .bed csaw Analyses préliminaires csaw Analyse différentielle Régions différentielles .txt Données alignées Données filtrées R / edgeR • Clustering • PCA • Meta-profiles • Normalisation
7. Aperçu d’un pipeline Bioconductor de traitement de données ChIP-seq Analyse de données de ChIP-seq sous R 7 Données brutes .fastq Données trimmées .fastq QuasR QuasR Rsamtools .bam .bam Prétraitement Données filtrées .bam Pics .bed csaw Analyses préliminaires csaw Analyse différentielle Régions différentielles .txt Visualisation Analyse de motifs Analyses fonctionnelles QuasR/ Gviz Données de visualisation .wig PWMEnrich Motifs enrichis .logo chipenrich Termes enrichis .txt Données alignées Données filtrées
8. Présentation des données Analyse de données de ChIP-seq sous R 8 TP Cellule de la lignée Ishikawa • Cellules épithéliales dérivées d’un adénocarcinome de l’endomètre Anticorps anti-ESR1 polyclonalde lapin (Santa Cruz Biotechnology sc-543) • Récepteur aux estrogènes, facteur de transcription présentant plusieurs sites d’interaction protéine – protéine, hormone – protéine et protéine – ADN • Localisation nucléaire • Forme des homo et hétéro-dimers avec ESR2 • Essentiel pour le développement des caractères liés au genre • Impliqué dans le développement de certains cancers du sein, de l’endomètre ainsi que l’ostéoporose ChiP-seq Illumina Genome Analyzer SE 36 bp • 2 x 107 cellules • Sonication • Amplification par 15 cycles de PCR
9. Présentation des données Analyse de données de ChIP-seq sous R 9 TP Echantillons • 2 réplicats biologiques + 0.02 % DMSO durant une heure • 2 réplicats biologiques + 10 nM oestradiol durant une heure • Pas de librairie contrôle (Sono-seq) Replicat Nom Séquences Totales (M) Séquences Chr19 (M) Estradiol Rep 1 ESR1_ES_REP1 18 0.50 Rep 2 ESR1_ES_REP2 23 0.46 DMSO Rep 1 ESR1_DMSO_REP1 21 0.50 Rep 2 ESR1_DMSO_REP1 26 0.40
10. ETAPE 1 : Préparation des données - Import des données brutes - Analyse de données de ChIP-seq sous R 10 TP Utilisez le package ShortRead pour procéder à l’import des données sous R
11. ETAPE 1 : Préparation des données - Import des données brutes - Analyse de données de ChIP-seq sous R 11 TP Explorez les données importées
12. ETAPE 1 : Préparation des données - Import des données brutes - Analyse de données de ChIP-seq sous R 12 TP Explorez les données importées
13. ETAPE 1 : Préparation des données - Qualité des données brutes - Analyse de données de ChIP-seq sous R 13 TP Utilisez le package Rqc pour évaluer la qualité des 4 libraires
14. ETAPE 1 : Préparation des données - Qualité des données brutes - Analyse de données de ChIP-seq sous R 14 TP Utilisez le package Rqc pour évaluer la qualité des 4 libraires
15. ETAPE 1 : Préparation des données - Qualité des données brutes - Analyse de données de ChIP-seq sous R 15 TP Utilisez le package Rqc pour évaluer la qualité des 4 libraires
16. ETAPE 1 : Préparation des données - Filtrer et trimmer - Analyse de données de ChIP-seq sous R 16 TP Utilisez le package QuasR pour filtrer et trimmer les lectures
17. ETAPE 1 : Préparation des données - Filtrer et trimmer - Analyse de données de ChIP-seq sous R 17 TP Evaluez l’efficacité du trimming avec Rqc
18. ETAPE 1 : Préparation des données - Filtrer et trimmer - Analyse de données de ChIP-seq sous R 18 TP Evaluez l’efficacité du trimming avec Rqc
19. ETAPE 1 : Préparation des données - Filtrer et trimmer - Analyse de données de ChIP-seq sous R 19 TP Evaluez l’efficacité du trimming avec Rqc
20. ETAPE 2 : L’alignement - Génome de référence - Analyse de données de ChIP-seq sous R 20 Téléchargez les FASTA associé au chromosome 19 de l’assemblage hg38 du génome humain sur Ensembl (version soft mask et mask) Importer les deux FASTA dans R via le package Biostrings TP
21. ETAPE 2 : L’alignement - Génome de référence - Analyse de données de ChIP-seq sous R 21 TP Explorez les données importées
22. ETAPE 2 : L’alignement - Génome de référence - Analyse de données de ChIP-seq sous R 22 TP Explorez les données importées
23. ETAPE 2 : L’alignement - Alignement avec Bowtie - Analyse de données de ChIP-seq sous R 23 TP Importez le fichier décrivant les échantillons à analyser Procédez à l’alignement avec Bowtie via QuasaR
24. ETAPE 2 : L’alignement - Qualité des alignements - Analyse de données de ChIP-seq sous R 24 TP Evaluez la qualité des alignements via QuasR
25. ETAPE 3 : Traitement des alignements - Tri et indexation - Analyse de données de ChIP-seq sous R 25 TP Procédez au tri et à l’indexation des alignement via Rsamtools
26. ETAPE 3 : Traitement des alignements - Manipuler les SAM / BAM sous R - Analyse de données de ChIP-seq sous R 26 • Volume d’information trop important pour charger la totalité d’un BAM dans R • Package GenomicAligments which : régions génomiques du BAM à considérer what : informations à extraire du BAM flag : caractéristiques des alignements à considérer ScanBamParam () which = what = flag=
27. ETAPE 3 : Traitement des alignements - Manipuler les SAM / BAM sous R - Analyse de données de ChIP-seq sous R 27 which : régions génomiques du BAM à considérer what : informations à extraire du BAM flag : caractéristiques des alignements à considérer ScanBamParam () • Volume d’information trop important pour charger la totalité d’un BAM dans R • Package GenomicAligments
28. ETAPE 3 : Traitement des alignements - Manipuler les SAM / BAM sous R - Analyse de données de ChIP-seq sous R 28 Importez un des BAM générés précédemment en utilisant ScanBamParam et readGAlignments TP Explorez les informations contenues
29. ETAPE 4 : Le peak-calling - Les régions blacklistées - Analyse de données de ChIP-seq sous R 29 Téléchargez le fichier BED référençant le régions blacklistées du génome hg19 Convertissez les coordonnées hg19 vers hg38 via liftOver TP
30. ETAPE 4 : Le peak-calling - Estimer la taille des fragments - Analyse de données de ChIP-seq sous R 30 Etablissez le profil de cross-corrélation pour chacune des librairies avec csaw TP
31. ETAPE 4 : Le peak-calling - Estimer la taille des fragments - Analyse de données de ChIP-seq sous R 31 Etablissez le profil de cross-corrélation pour chacune des librairies avec csaw TP
32. ETAPE 4 : Le peak-calling - Peak-calling sans contrôle - Analyse de données de ChIP-seq sous R 32 Principe du peakcallingsans contrôle avec csaw Comptage Comptage Signal / Bruit Valeur seuil Alignements Evaluation du signal Evaluation du bruit Ratio Signal / Bruit median
33. ETAPE 4 : Le peak-calling - Peak-calling sans contrôle - Analyse de données de ChIP-seq sous R 33 Générez les matrices de comptages dans des fenêtres chevauchantes de 150 bp Générez les matrices de comptages dans des fenêtres non-chevauchantes de 2000 bp Evaluez l’enrichissement local TP
34. ETAPE 4 : Le peak-calling - Peak-calling sans contrôle - Analyse de données de ChIP-seq sous R 34 Définissez le seuil discriminant signal / background
35. ETAPE 4 : Le peak-calling - Peak-calling sans contrôle - Analyse de données de ChIP-seq sous R 35 Définissez le seuil discriminant signal / background Filtrez les régions pour lesquelles le ratio signal / background < log2 (2.5) TP
36. ETAPE 4 : Le peak-calling - Peak-calling sans contrôle - Analyse de données de ChIP-seq sous R 36 Comparez les 4 échantillons TP
37. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 37 Spécificité de la normalisation pour les données de ChIP-seq Pourquoi normaliser Différence de taille entre les librairies (composition) Différence d’efficacité de l’immuno-précipitation (efficacité) Variabilité technique classique Comment normaliser • Pour la composition : • Normalisation par le nombre de reads totaux (RPM) • Scaling factors (TMM & DESeq2) • Pour l’efficacité : • Scaling factors (TMM & DESeq2) • Normalisation non-linéaire (Lowess, quantiles) A prendre en compte Mutuellement exclusives
38. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 38 Spécificité de la normalisation pour les données de ChIP-seq • Est-ce que l’on considère que les différences systématiques dans les régions où les comptages sont élevés sont de vrais positifs ? • OUI : Normalisation pour la composition Permet de préserver ces différences observées • NON : Normalisation pour l’efficacité Permet de lisser le différentiel d’efficacité de l’IP • Dépend de la question biologique : • KO vs wild-type : nombreuses différences attendues - Normalisation pour la composition • Profilage d’histones : peu de différences attendues - Normalisation pour l’efficacité Hypothèse : La plupart des pics ne sont pas différentiels
39. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 39 Spécificité de la normalisation pour les données de ChIP-seq Evaluer l ’effort de normalisation Ratio des comptages Moyenne des comptage Composition Efficacité MA plot Rep 1 vs Rep 2 Signal Background MA plot Condition 1 vs Condition 2 Régions différentielles Régions non différentielles
40. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 40 Spécificité de la normalisation pour les données de ChIP-seq Ligne matérialisant la normalisation pour la composition à passe à travers le nuage « background » Ligne matérialisant la normalisation pour l’efficacité à passe à travers le nuage « signal » Efficacité Composition MA plot Rep 1 vs Rep 2 Signal Background MA plot Condition 1 vs Condition 2 Régions différentielles Régions non différentielles
41. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 41 Spécificité de la normalisation pour les données de ChIP-seq M ≠ 0 pour le background à biais de composition M ≠ 0 pour le signal à différence systématique pour les régions très couvertes, artefactuelle ou non biais d’efficacité possible Efficacité Composition MA plot Rep 1 vs Rep 2 Signal Background MA plot Condition 1 vs Condition 2 Régions différentielles Régions non différentielles
42. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 42 Evaluez les facteurs de normalisation (composition et efficacité) TP
43. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 43 Evaluez l’effort de normalisation TP
44. ETAPE 5 : Analyse différentielle - Normalisation - Analyse de données de ChIP-seq sous R 44 Evaluez l’effort de normalisation DMSO Rep1 vs Rep2 ES Rep1 vs Rep2 Composition Efficacité TP
45. ETAPE 5 : Analyse différentielle - edgeR - Analyse de données de ChIP-seq sous R 45 Procédez à l’analyse différentielle avec EdgeR TP Common Tagwise Trended
46. ETAPE 5 : Analyse différentielle - edgeR - Analyse de données de ChIP-seq sous R 46 Procédez à l’analyse différentielle avec EdgeR TP • BCV = racine carrée de la dispersion (loi NB) • Dispersion tagwise à variabilité importante pour les régions à faible comptage • Dispersion common à 20 %, les comptages sont très reproductibles, peu de variabilité • Dispersion trended à le nuage de point est réparti au dessus et en dessous, suggérant que le modèle a été correctement ajusté • Pour la plupart des jeux de données, on s’attend à une tendance qui décroit pour se stabiliser à un plateau à plus grande robustesse pour les comptages plus importants Common Tagwise Trended
47. ETAPE 5 : Analyse différentielle - Ajustement pour les tests multiples - Analyse de données de ChIP-seq sous R 47 Clusterisezles fenêtres adjacentes et évaluer le FDR Définissez, pour chaque cluster, la fenêtre la plus représentative TP
48. ETAPE 5 : Analyse différentielle - Ajustement pour les tests multiples - Analyse de données de ChIP-seq sous R 48 Préparez un GRanges synthétisant les résultats Sélectionnez les clusters différentiels TP
49. ETAPE 5 : Analyse différentielle - Sortie BED et WIG - Analyse de données de ChIP-seq sous R 49 Exportez les résultats de l’analyse différentielle et les données de couverture TP
50. ETAPE 5 : Analyse différentielle - Sortie BED et WIG - Analyse de données de ChIP-seq sous R 50 Qu’est ce que les formats WIG / BIGWIG / BEGRAPH ? • Format utiliser pour stocker des informations de visualisation denses et continues (e.g. contenu en GC, probabilité, couverture en reads) • BIGWIG : format binaire, plus léger et en général recommandé • WIG : format ASCII, données continues • BEDGRAPH : format ASCII, données sparses Qu’est ce que le format BED ? • « Browser Extensible Data » • Format tabulaire adapté à l’affichage de données d’intervalles sur les « genome browsers » (UCSC, Ensembl, IGV …)
51. ETAPE 6 : Visualisation & interprétation - Visualisation avec Gviz - Analyse de données de ChIP-seq sous R 51 Préparez les données d’annotation utilisées par Gviz TP
52. ETAPE 6 : Visualisation & interprétation - Visualisation avec Gviz - Analyse de données de ChIP-seq sous R 52 Préparez les données relatives aux échantillons TP
53. ETAPE 6 : Visualisation & interprétation - Visualisation avec Gviz - Analyse de données de ChIP-seq sous R 53 Générez la sortie graphique Ideo G_Axe HT G_Reg TP
54. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 54 Générez une matrice de comptage normalisée (centre des pics ± 1000 bp) TP
55. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 55 Générez une matrice de comptage normalisée (centre des pics ± 1000 bp) TP
56. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 56 Procédez à un lissage par fenêtre glissante TP
57. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 57 Calculez les vecteurs de couverture moyenne pour les pics « UP » et « DOWN » TP
58. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 58 Générez les méta-profils avec ggplot2 TP
59. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 59 Préparez les données d’annotation pour les heatmaps TP
60. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 60 Créez une fonction personnalisée basée sur ComplexHeatmap pour générez les heatmaps TP
61. ETAPE 6 : Visualisation & interprétation - Meta-profils et binding heatmaps - Analyse de données de ChIP-seq sous R 61 Générez les heatmaps TP
62. ETAPE 6 : Visualisation & interprétation - Enrichissement en motifs - Analyse de données de ChIP-seq sous R 62 Convertissez les coordonnées des pics différentiels de hg38 vers hg19 TP
63. ETAPE 6 : Visualisation & interprétation - Enrichissement en motifs - Analyse de données de ChIP-seq sous R 63 Identifiez les motifs les plus enrichis pour les pics « UP » TP
64. ETAPE 6 : Visualisation & interprétation - Enrichissement fonctionnel - Analyse de données de ChIP-seq sous R 64 Pourquoi utiliser des outils adaptés au ChIP-seq ? FET et dérivés DAVID, METASCAPE Chaque gène a une probabilité identique d’être « significatif » Proportion similaire de gènes « significatifs » dans chaque gene set testés • NGS : corrélation positive entre la longueur du gène et la probabilité de le définir « significatif » • ChIP-seq : la probabilité qu’un gène présente un pic est proportionnelle à sa longueur • Les gènes plus longs apportent une contribution disproportionnée au signal d’enrichissement à Hypothèses du FET invalidées • Par ailleurs, la distribution en taille n’est pas homogène dans les gene sets Nervous system development (GO BP) rRNA processing (GO BP) Hypothèses
65. ETAPE 6 : Visualisation & interprétation - Enrichissement fonctionnel - Analyse de données de ChIP-seq sous R 65 Pourquoi utiliser des outils adaptés au ChIP-seq ? Test binomial GREAT Le nombre de pics dans un locus est proportionnel à la taille du locus La variabilité en terme de nombre de pics par gène est proportionnelle à la taille du locus • 63 ENCODE ChIP-seq GM12878 Hypothèses Ryan P. Welch, Chee Lee, Paul M. Imbriano, Snehal Patil, Terry E. et al. “ChIP-Enrich: gene set enrichment testing for ChIP-seq data.” NAR 2014
66. ETAPE 6 : Visualisation & interprétation - Enrichissement fonctionnel - Analyse de données de ChIP-seq sous R 66 Pourquoi utiliser des outils adaptés au ChIP-seq ? Chip-Enrich Ajustement pour la relation entre taille du locus et nombre de pics Ajustement pour la mappabilité Hypothèses Ryan P. Welch, Chee Lee, Paul M. Imbriano, Snehal Patil, Terry E. et al. “ChIP-Enrich: gene set enrichment testing for ChIP-seq data.” NAR 2014 Modèle 1 – FET Modèle 2 – Binomial Modèle 3 – Chip-Enrich Observations Relationship between the probability of a gene being assigned a peak and locus length.
67. ETAPE 6 : Visualisation & interprétation - Enrichissement fonctionnel - Analyse de données de ChIP-seq sous R 67 Pourquoi utiliser des outils adaptés au ChIP-seq ? Chip-Enrich Ajustement pour la relation entre taille du locus et nombre de pics Ajustement pour la mappabilité Hypothèses Ryan P. Welch, Chee Lee, Paul M. Imbriano, Snehal Patil, Terry E. et al. “ChIP-Enrich: gene set enrichment testing for ChIP-seq data.” NAR 2014 Modèle 1 – FET Modèle 2 – Binomial Modèle 3 – Chip-Enrich Observations Relationship between the probability of a gene being assigned a peak and locus length.
68. ETAPE 6 : Visualisation & interprétation - Enrichissement fonctionnel - Analyse de données de ChIP-seq sous R 68 Procédez au test d’enrichissement fonctionnel en utilisant chipenrich TP