Measuring methylation from arrays to sequencing Jovana Maksimovic

Measuring methylation: from arrays to sequencing Jovana Maksimovic, Ph. D jovana. maksimovic@mcri. edu. au @Jov. Maksimovic github. com/Jov. Maksimovic Bioinformatics Winter School, 3 July 2017

Talk outline • • Epigenetics DNA methylation Measuring DNA methylation Methylation arrays • How do they work? • What do they measure? • Example analysis • Methylation sequencing • What are the challenges? • How does it work? • Suggested analysis pipeline • Summary

I mostly work on human development & disease… I work at MCRI Me! …sometimes using mice or other models Ch. IPseq ATACseq BSseq Microarrays I analyse a lot of epigenetic data… I write software for analysing methylation array data miss. Methyl RNAseq Microarrays …and a lot of gene expression data

What is epigenetics? • Epigenetics refers to stable heritable traits not explained by changes in DNA sequence • Greek prefix “epi” means “on top of” genetics • Chromosome modifications that affect gene expression • Histones, DNA methylation • “Anything” that isn’t DNA! • Essential for normal development • Can be modified by environment • Can be disrupted in disease

Epigenetics brings DNA to life! identical DNA in every cell gland cell hormonesecreting cell B cell red blood cell haematopoietic stem cell fat cell embryogenesis sperm intestine cell kidney cell zygote muscle cell blastocyst egg lung germ cell embryonic stem cells neuronal progenitor skin cell • Important in all species sperm cell different epigenetic patterns T cell neuron astrocyte Modified from https: //biology. mit. edu/research/stemcell_epigenetics

Epigenetics is CRAZY complicated! Roy et al. (2010), Science Me Me • New sequencing & microarray technologies are enabling us to learn A LOT more about epigenetics • Different data types need different analysis • Today I’m only focussing on DNA methylation

What is DNA methylation? C C T C G A DNA methylation primarily occurs at Cp. G dinucleotides

DNA methylation in the genome • The human genome contains ~30, 000 Cp. Gs (~1%) • VERY different between different species • Cp. Gs are not evenly spaced across the genome Patterson et al. 2011, J Vis Exp Methylation correlation with distance • Tend to be present in clusters called Cp. G islands • Cp. G methylation is spatially correlated ~500 bp Eckhardt et al. 2007, Nature Genetics

Methylation can regulate gene expression Methylation at a single Cp. G vs. gene expression http: //meeting. dxy. cn/oemethylation 2012/article/i 18782. html Each point is one sample Plot from Peter Hickey

Methylation changes coat colour of Agouti mice Dolinoy 2008, Nutr Rev. This gene controls coat colour in Agouti mice These Cp. G sites in the promoter change PS 1 A expression depending on methylation These mice are genetically identical Hypomethylated Hypermethylated Coat colour different due to different maternal diet i. e. environment!

Methylation makes worker bees! These larvae are genetically identical Cridge et al. 2015, Nutrients Hypomethylated Hypermethylated

Methylation is cool What do we usually want to know about it?

Finding methylation differences can tell us a lot • Methylation is critical in determining cell type • Regulatory T-cell vs. Naïve T-cell • Methylation can be disrupted in disease • Cancer vs. Normal • Methylation is affected by the environment • Smokers vs. Non-smokers Collect appropriate samples Normal Cancer Extract DNA and measure methylation Statistical analysis

Epigenome-wide association studies (EWAS) • Similar to GWAS • Compare lots of cases to lots of controls • Often looking for small effects e. g. complex disease or environmental effects • Need lots of samples • 100 s or 1000 s of cases & controls https: //en. wikipedia. org/wiki/Epigenome-wide_association_study_(EWAS)

How do we measure methylation? • Bisulphite conversion • Create “SNPs” • Single nucleotide resolution • Array • Sequencing • Enrichment of methylated DNA • Restriction enzymes • Affinity • Regional resolution • Array • Sequencing

What is bisulphite conversion? PCR • Chemical process • Unmethylated Cs get converted to Ts • Methylated Cs are protected • Creates “SNP” • Used to call methylation

Methylation arrays What are they and how do they work?

Illumina Infinium Human. Methylation Bead. Chips 27 k array (2009) 1 chip = 12 samples >27, 000 unique Cp. G sites measured in each sample 450 k array (2011) 1 chip = 12 samples >450, 000 unique Cp. G sites measured in each sample 850 k array (2015) 1 chip = 8 samples >850, 000 unique Cp. G sites measured in each sample • Human only • Gene biased; selected to be relevant to human development & disease • eg. TSS, promoters, Cp. G islands, enhancers, . . . Modified slide from Belinda Phipson

Methylation arrays are based on SNP array technology What is this base? Measure fluorescence intensity • Methylation array “SNPs” (C/T) are created by bisulphite conversion • Comparing the intensity of C/T gives the proportion of methylation at single Cp. G

What methylation values can we get? A sample 1 0. 5 0 CH 3 Many cells in single sample • On an array, we measure methylation in a population of cells • Individual cell can be either 0, 0. 5 or 1 at one Cp. G • Across a population we get a continuous measurement between [0 -1]

Measures of methylation • Arrays measure both methylated (C) and unmethylated (T) signal to get proportion of methylation at a Cp. G Beta value Du et al. 2011, BMC Bioinformatics Intuitive, easy to interpret, great for visualisation Better statistical properties, recommended for statistical testing M value Can convert between them via a logit transformation

What does the data look like? Table of M-values Cp. G sites Sample A 1 Sample A 2 Sample A 3 Sample B 1 Sample B 2 Sample B 3 0. 213 0. 221 0. 311 0. 123 0. 216 0. 198 -0. 011 0. 001 -0. 016 2. 011 2. 002 2. 702 2. 213 2. 256 2. 698 0. 052 0. 101 0. 238 4. 567 5. 231 4. 982 4. 152 6. 216 4. 698 -4. 723 -3. 459 -5. 36 -5. 763 -5. 122 -4. 998 -5. 567 -4. 666 -4. 845 -4. 522 -4. 111 -3. 245 3. 421 5. 467 5. 554 5. 445 5. 298 4. 514 2. 981 3. 345 3. 512 -3. 534 -4. 311 -3. 889 3. 792 2. 987 3. 324 -0. 231 -0. 066 -0. 001 …. . .

Array analysis pipeline QC: b density plots, control probes, MDS/clustering plots, … Normalization: within and between arrays Software Remove bad samples and poor performing probes (Cp. Gs) minfi, methylumi, limma Transform data to remove unwanted variation minfi, miss. Methyl, wate. Rmelon Statistical testing for differential methylation, Cp. Gs & regions Estimate means and variances and limma, bumphunter, borrow information DMRcate across probes Annotation to genes, gene set testing, visualization, … Think about biological interpretation miss. Methyl, Gviz e. g. gene expression Genomic. Ranges Combine with other data types

M 28 M 29 M 30 r. Treg naive activated r. Treg activated naive

After QC, data exploration is your friend! Dimension 4 Dimension 2 MDS plots showing largest sources of variation in the data Dimension 1 Dimension 3 Clustering by individual and cell type

Statistical testing: Look for differences at single Cp. Gs Differential methylation Linear model : Can take into account any other covariates Adjust the p-values using Benjamini and Hochberg’s FDR Phipson & Oshlack 2015, Genome Biology Lots of differences between immune cell types! One test per Cp. G! Smyth, 2004 Modified slide from Belinda Phipson & Alicia Oshlack

Statistical testing: Differences across Cp. G dense region • Recall: Cp. G methylation is spatially correlated • Can we find consistent group-average level differences between Cp. Gs that are close together? • More functionally relevant than differences at individual Cp. Gs? Lots of DMRs between immune cell types! Aryee et al. 2014, Bioinformatics

You can do other cool stuff! • Unmethylated regions in r. Treg compared to naïve cells enriched for FOXP 3 binding motifs! DMR consensus motif matches Forkhead-binding motif Forkhead- binding motif Consensus motif from DMR seqs. Differences in cell types controlled by FOXP 3! Modified slide from Alicia Oshlack

Methylation array analysis is very mature: lots of methods! https: //www. bioconductor. org/ https: //f 1000 research. com/articles/5 -1281/v 3

Methylation sequencing AKA bisulphite sequencing: the good, the bad and the ugly

Two main types of bisulphite sequencing • Whole-genome bisulphite sequencing (BS-seq) • Gold standard • Genome-wide (~30, 000 Cp. Gs in human) • Expensive but covers almost everything • Need high (10 -30 x) coverage to reliably call methylation • Targeted BS-seq • Only sequence regions of interest • Reduced representation BS-seq (restriction enzyme) • Capture BS-seq (similar principal to exome) • Cheaper but can miss a lot of stuff • Can usually do higher (20 -60 x) coverage

What was bisulphite conversion again? DNA fragment All four of these can be sequenced!

What are the challenges? • Like calling SNPs, methylation in BS-seq inferred by comparison to unconverted reference sequence • Correct alignment is critical • More challenging than usual! • Aligned sequences do not exactly match reference • Complexity of libraries is reduced • Many Cs become Ts, so less info for mapping! • Methylation is not symmetrical • Two strands of DNA in the reference genome must be considered separately

DNA fragment Mapping (Bismark) BS conversion & PCR

DNA fragment Mapping (Bismark) BS conversion & PCR TCGGTATGTTTAAACGTT

DNA fragment Mapping (Bismark) BS conversion & PCR TCGGTATGTTTAAACGTT In silico read conversion C-to-T G-to-A TTGGTATGTTTAAATGTT TCAATATATTTAAACATT

DNA fragment Mapping (Bismark) BS conversion & PCR TCGGTATGTTTAAACGTT In silico read conversion Align to in silico bisulphite converted genome C-to-T G-to-A TTGGTATGTTTAAATGTT TCAATATATTTAAACATT Fwd strand C-to-T converted genome Fwd strand G-to-A converted genome …TTGGTATGTTTAAATGTT… …AACCATACAAATTTACAA… …CCAACATATTTAAACACT… …GGTTGTATAAATTTGTGA… Reverse complement

DNA fragment Mapping (Bismark) BS conversion & PCR TCGGTATGTTTAAACGTT In silico read conversion Align to in silico bisulphite converted genome C-to-T G-to-A TTGGTATGTTTAAATGTT TCAATATATTTAAACATT Fwd strand C-to-T converted genome Fwd strand G-to-A converted genome …TTGGTATGTTTAAATGTT… …AACCATACAAATTTACAA… …CCAACATATTTAAACACT… …GGTTGTATAAATTTGTGA… Reverse complement Read all alignment outputs simultaneously to determine if sequence can be mapped uniquely Reverse complement …TTGGTATGTTTAAATGTT… TCAATATATTTAAACATT …CCAACATATTTAAACACT… TCAATATATTTAAACATT x x x …AACCATACAAATTTACAA… TCAATATATTTAAACATT …GGTTGTATAAATTTGTGA… TCAATATATTTAAACATT x x x x x x x x x x x

DNA fragment Mapping (Bismark) BS conversion & PCR TCGGTATGTTTAAACGTT In silico read conversion G-to-A TTGGTATGTTTAAATGTT TCAATATATTTAAACATT Fwd strand C-to-T converted genome Fwd strand G-to-A converted genome …TTGGTATGTTTAAATGTT… …AACCATACAAATTTACAA… …CCAACATATTTAAACACT… …GGTTGTATAAATTTGTGA… Reverse complement Read all alignment outputs simultaneously to determine if sequence can be mapped uniquely Reverse complement …TTGGTATGTTTAAATGTT… TTGGTATGTTTAAATGTT TCAATATATTTAAACATT …CCAACATATTTAAACACT… TTGGTATGTTTAAATGTT TCAATATATTTAAACATT x xx xx x x xx …AACCATACAAATTTACAA… TTGGTATGTTTAAATGTT TCAATATATTTAAACATT …GGTTGTATAAATTTGTGA… TTGGTATGTTTAAATGTT TCAATATATTTAAACATT x x x xx xx xx xx xx x Align to in silico bisulphite converted genome C-to-T xx xx xx xx xx

Calling methylation …CCGGCATGTTTAAACGCT… TCGGTATGTTTAAATGTT …TTG TATGTTTAAATGTT …TCGGTATGTTTAAAT GTT… Genome reference …TCGGTATGTT AAACGTT… …TCGGTATGTTTAAATGTT …TCGGTATGTTT ATGTT… …TCGGTATGTTTAAAT TT… …TTGGTATGTTTA ATGTT… …TCGGTATGTTTAAACGT

Calling methylation …CCGGCATGTTTAAACGCT… TCGGTATGTTTAAATGTT …TTG TATGTTTAAATGTT Genome reference …TCGGTATGTTTAAATGTT… Good coverage is very important for reliable methylation calls!

Some real BS-seq mapping results https: //software. broadinstitute. org/software/igv/interpreting_bisulfite_mode

Methylation calling output Position of C in genome No. methylated reads chr 1 chr 1 chr 1 chr 1 chr 1 753479 753492 753540 753541 753667 753724 753763 753785 759932 760913 761299 761371 761377 761446 761460 762005 762114 762176 762180 This is what we work with! % methylation 50 66. 67 100 50 25 66. 67 0 0 100 80 100 92. 86 53. 85 100 0 1 2 1 1 1 2 0 0 1 0 2 8 10 13 7 1 0 0 0 1 1 0 1 3 1 2 1 0 2 0 1 6 0 5 7 8 Sum for total coverage No. unmethylated reads

Analysis pipeline Krueger et al. 2012, Nature Methods Thorough QC is VERY important for BS-seq …and adapters Need to be brutal with trimming off poor quality bases… As with SNP calling, removing PCR duplicates is a good idea for better methylation calling Other stuff to find cool biology!

Summary • Methylation arrays very popular • Only for human • Great for EWAS • Analysis very mature • Bioconductor is the place to go! • BS-seq best option for genome-wide single nucleotide resolution • Only option for species other than human • Pre-processing, mapping, etc. pretty good • Statistical analysis still developing • Bioconductor is a valuable resource • Downstream analysis dependent on biological question • Methylation is interesting & we know how to measure it • Best technology for the job depends on what you want to know!

Acknowledgments Murdoch Childrens Research Institute https: //f 1000 research. com/articles/5 -1281/v 3 • Alicia Oshlack • Belinda Phipson • MCRI Bioinformatics group! Johns Hopkins University • Peter Hickey jovana. maksimovic@mcri. edu. au @Jov. Maksimovic github. com/Jov. Maksimovic miss. Methyl https: //www. bioconductor. org/packages/ release/bioc/html/miss. Methyl. html