Targeted bisulfite sequencing for biomarker discovery

Abstract

Cytosine methylation is one of the best studied epigenetic modifications. In mammals, DNA methylation patterns vary among cells and is mainly found in the CpG context. DNA methylation is involved in important processes during development and differentiation and its dysregulation can lead to or is associated with diseases, such as cancer, loss-of-imprinting syndromes and neurological disorders. It has been also shown that DNA methylation at the cellular, tissue and organism level varies with age. To overcome the costs of Whole-Genome Bisulfite Sequencing, the gold standard method to detect 5-methylcytosines at a single base resolution, DNA methylation arrays have been developed and extensively used. This method allows one to assess the status of a fraction of the CpG sites present in the genome of an organism. In order to combine the relatively low cost of Methylation Arrays and digital signals of bisulfite sequencing, we developed a Targeted Bisulfite Sequencing method that can be applied to biomarker discovery for virtually any phenotype. Here we describe a comprehensive step-by-step protocol to build a DNA methylation-based epigenetic clock.

Keywords: Biomarker discovery; DNA methylation; Epigenetic clock; Next-generation sequencing; Target bisulfite-seq.

Figure 1:. Overview of the Targeted Bisulfite Sequencing Protocol.

Genomic DNA is extracted from collected blood samples (section 4), fragmented and subject to NGS library preparation (5.1 5.10). Adapter ligated libraries are then pooled (5.10.14), concentrated (6.1) and incubated with RNA biotinylated probes to enrich for target regions (6.2 6.33 and Figure 2). Captured DNA fragments are then bisulfite-treated (7.1 7.17), PCR amplified, Quality Controlled (7.18 7.22) and sequenced (7.23). Data is then analyzed (see Section 8 and Figure 4 for more details).

Figure 2:. Hybridization Capture setup.

The Hybridization mix (top, yellow boxes – 2.2.7) is prepared and incubated in a Thermomixer for 10 min at 60°C (6.4). After 5 min of Room Temperature (RT) incubation (6.5), 18.5 μl are transferred into PCR tubes (6.8). The HYB tubes are stored at RT until step 6.9. Dried library pools from step 6.1 are resuspended for 15 min at 60°C in a Thermomixer after the addition of 7 μl of H2O and 5 μl of Blockers Mix (2.2.6) (Step 6.2). The blockers + libraries mixture (12 μl) is then transferred into PCR tubes that are now called LIBs (6.6). LIB tubes are then transferred into a thermocycler for denaturation and hybridization temperature equilibration (step 6.7). After the thermocycler reaches the hybridization temperature (65°C), HYB tubes are added into the machine (6.9). After 5 minutes at 65°C, 18 μl of the HYB mix are transferred into each LIB tube (still in the thermocycler) (6.10). HYB tubes are discarded, and the LIB+HYB tubes are incubated for 16–20 hours at the hybridization temperature (6.11).

Figure 3:. Final Libraries Quality Control.

Agilent TapeStation 2200 High Sensitivity D1000 ScreenTape Assay (7.21) of the Final Libraries obtained with the Targeted Bisulfite Sequencing Approach (7.19.14). (a) Gel View of the Ladder and Final Libraries. (b) Electropherogram View with information about the Dimers Region (orange, 60–160 bp) and the Final Libraries Region (blue, 180–1000 bp).

Figure 4:. Data Analysis Pipeline described in Section 8.

FastQ files (8.1.1) are obtained from NGS-sequencing (7.23) and subject to quality control (FastQC) before and after Adapter Trimming (cutadapt) (8.2.3). Genome fasta files can be obtained from sequence databases (8.1.2) and then subject to Genome Indexing (BSBolt Index – 8.2.2). Both Genome Index and trimmed FastQ files are used as input for Reads alignment (BSBolt Align – 8.2.4). After Alignment, duplicated reads are removed (samtools markdup – 8.2.5) and methylation is called for every cytosine (BSBolt CallMethylation – 8.2.6), creating CGmap files. Several CGmap files can be combined into a single matrix using BSBolt AggregateMatrix (8.3). Methylation values for missing sites can be imputed using BSBolt Impute (OPTIONAL – 8.3.2). Model fitting is then performed in a Jupyter Notebook (8.4.2), where external libraries (8.4.3) and phenotypic data (8.4.4) are imported. After Samples QC to detect outliers (8.4.5), the model is trained using a leave-one-out elastic-net penalized regression (8.4.6). The model is then evaluated by fitting a trendline between the known and the predicted phenotypic value (8.4.7).

Figure 5:. Targeted Bisulfite Sequencing Epigenetic Clock.

Epigenetic age predictions for (n=48) samples made using penalized regression models compared to the chronological age of each sample with a line of best fit. The chronological (observed) age is represented on the x-axis, while the predicted epigenetic age is on the y-axis.