Tutorial for MassArray Primer Design

Targeted DNA methylation analysis by MassArray

The MassARRAY® System by Agena Bioscience allows a mass-spectrometry based analysis of DNA methylation patterns with single CpG resolution. The system offers a low-cost alternative and allows users to analyze DNA methylation of even formalin-fixed paraffin-embedded (FFPE) tissue. In short, the procedure begins with a bisulfite conversion of isolated genomic DNA from samples of interest. Thereby, unmethylated cytosines are converted to uracil residues while methylated cytosines are not affected. This ensures a difference in mass which can later on be analyzed by mass-spectrometry. A PCR amplification of the region of interest with a T7-promoter-tagged reverse primer is performed and the resulting fragments are in vitro transcribed. The latter process leads to an incorporation of uracil residues which generates a RNase cleavage pattern in the respective fragments. The fragments can be analyzed by the MassArray MALDI-TOF spectrometer and differ in mass due to the bisulfite treatment.

Source: “http://agenabio.com/wp-content/uploads/2015/06/51-20055R1.0-EpiTYPER-Brochure_WEB.pdf”

MassArray primer design

The MassArray primer design process requires several prerequesits compared to normal DNA primer design. First, it has to be assessed whether the analyzed CpGs are located within amplicons which show a mass difference in the MassArray workflow. CpGs which are located in fragments with an overlapping molecular weight (MW) cannot be analyzed as users will not be able to differentiate between CpGs in fragments with a mass-overlap. Second, it has to be determined which strand (Watson or Crick) has more CpG containing fragments which can be analyzed by MassArray. Bisulfite treatment destroys the complementary of both strands so that primers have to be specifically designed for either the coding (Watson) or non-coding (Crick) strand. And third, primers have to be designed for bisulfite converted DNA which requires adjusted primer design parameters.
Considering those prerequesits, MassArray primer design can be considered a tedious process. Users have to extract the genomic DNA sequence for the amplicons of interest, bisulfite convert the conding and non-coding strand in-silico, assess wether CpGs can be detected by MassArray and then design primer pairs. AmpliconDesign offers an all-in-one solution which extracts genomic DNA sequences from common reference genomes, analyzes fragment cleavage patterns, bisulfite converts DNA in-silico and returns a list of suggested primer pairs. Users solely have to provide the genomic coordinates of interest. A graphical output is generated which allows users to select the best primer pairs for MassArray.

Tutorial: Primer Design for MassArray

1.) Select the required genome build

AmpliconDesign currently supports three common genome builds: GRCh38/hg38, GRCh27/hg19 and GRCm38/mm10.
Select the required genome in the drop-down menu.

2.) Enter the coordinates of the genomic region of interest

The genomic regions of interest have to be entered in the form of genome coordinates (e.g. chr19:43203328-43203389)
Genomic regions can be further extended within the primer design workflow.
Do neither comma separate the numbers nor add a white space before the coordinates

3.) Inspect your Input Region

AmpliconDesign extract the region of interest from the reference genome and provides the user with the DNA sequence of the Watson (forward sequence) and Crick strand (reverse complement sequence). Users can furthermore inspect their input region in a tabular format. A slider allows the extension of the input region. The shown DNA sequences provide the user furthermore with information about repeat regions (green), SNP positions (violet) and CpG dinucleotides (red).

4.) Plot with Genome Annotations

Users can inspect the analyzed genomic region in a plot which includes transcript, SNP, CpG and repeat positions. The page furthermore includes a tabluar output which lists all the CpGs in the analyzed fragment.

5.) Analyze the in-silico Fragmentation

The amplicon prediction plots have been adapted from the MassArray R package (Thompson RF, Greally JM (2018). MassArray: Analytical Tools for MassArray Data. R package version 1.34.0.). The initial table shows in a binary format which of the covered CpGs can be analyzed by MassArray. The “summary” column inidicates whether a CpG is covered by putative fragments (1 = detectable; 0 = not detectable). The other columns show in which cleavage (C or T cleavage) and on which strand (+ or -) the CpG can be analyzed. It is crucial that users furthermore check the amplicon prediction plots below. Those show the putative fragments and the location of CpGs (filled circle). CpGs within a blue fragment can be uniquely analyzed in the MassArray workflow. Red fragments show overlaping molecular weights (MWs) with another fragment. Grey fragments are outside of the testable mass window. Black fragments do not contain a CpG but are uniquely assayable. Green fragments might contain a conversion control. An overlapping MW between CpG containing fragments is shown by linked arrowheads.

6.) Manual Primer Design

AmpliconDesign offers the possibility to manually design MassArray primers. Therefore users can select between the Watson and Crick strand, depending on the number of CpGs which are covered. In general it is recommended to select the strand with the highest number of fragments that can be analyzed by MassArray. A selector for the coding strand and one for the non-coding strand allows to mark the putative primer sequence in the genomic region. Please make sure that primer sequences do not overlap with CpGs or SNP (if unavoidable place them in the 5' position). The first sequence always shows the 5'-3' direction of the coding strand. So if the Watson strand is selected, this corresponds to the bisulfite converted coding strand. If the Crick strand is selected, this corresponds to the reverse complement of the Crick strand. This visualization is required to ensure the integrity of the genome annotations. The second sequence always shows the 3'-5' orientation. In case of the Watson selection, the complement of the bisulfite converted coding sequence and in case of the Crick, the bisulfite converted reverse Crick strand (3'-5').

The selected primer pair is shown in a table in a ready to order fashion. The table shows the forward primer and reverse primer in the correct orientation. Please ensure to add a T7-promoter tag to the primer of choice before ordering the primers for MassArray. The table furthermore shows primer sizes and the final amplicon size.

7.) Automated Primer Design

AmpliconDesign furthermore offers the possibility to automatically design primers by calling primer3 (Untergasser et al., NAR, 2012) with bisulfite adjusted parameters. Several input parameters have to be specified to allow a user-specific primer design:

Primer Size

Choose the optimal, minimal and maximal primer size.

Melting Temperature

Choose the optimal, minimal and maximal primer melting temperature.

Amplicon Size

Choose the minimal and maximal amplicon size.

Exclude CG from primer

Recommended: Choose this option to exclude CG dinucleotides from the primer sequence.

Include CpGs in the region

Users can choose if the covered CpGs should be included in the amplicon.

8.) Table with automatically designed primers

AmpliconDesign generates a table with automatically designed primers. This table includes the primer sequence for the plus strand (Amplicon 1) and the minus strand (Amplicon 2). The primer ID is given as the specified input region. Primer sequences for the forward and reverse primer are predicted (whereby t corresponds to bisulfite converted C). The primer begin is given with respect to the amplicon start site, followed by the GC content of the primer, melting temperatures and the total amplicon size. The Region column repeats the input region, and the UCSC coordinates for the forward and reverse primer (left and right) are given, as well as the final amplicon coordinates and the number of CpGs in each amplicon (#CpGs).

An overview of the putative fragment patterns, including MW, number of covered CpGs and collision (MW overlap with another fragment) are shown for each primer pair when clicking on a row in the primer table. Thereby users can check if their CpG of interest can be analyzed by MassArray.

9.) Download of the data

AmpliconDesign allows a download of the data and provides the user with all previously specified information, including selected primer sequences and a table with the automatically designed primer pairs.

Interactive Data

Data can be loaded and interactively explored.

Tutorial for AmpBS-Seq Primer Design

Amplicon Bisulfite Sequencing

Amplcion Bisulfite Sequencing (ampBS-Seq) allows a targeted in-depth analysis of genomic regions by next-generation sequencing. It is currently considered as state-of-the-art for DNA methylation biomarker development and is applied in many clinical and academic research laboratories worldwide (Bock et al., Nature Biotechnology, 2016). In short, the workflow starts with the bisulfite conversion of genomic DNA, followed by PCR amplification of the regions of interest. PCR primers must be specifically designed for bisulfite converted DNA and should not contain CpG dinucleotides. The PCR amplicons of each sample are pooled and barcoded which allows the analysis of several sample on one sequencing lane. A bioinformatic analysis can then distinguish between methylated and unmethylated CpGs due to the bisulfite conversion of unmethylated CpG to uracil.

Tutorial: Primer Design for ampBS-Seq

1.) Select the required genome build

AmpliconDesign currently supports three common genome builds: GRCh38/hg38, GRCh27/hg19 and GRCm38/mm10.
Select the required genome in the drop-down menu.

2.) Select the input

AmpliconDesign supports CpG-IDs, Genome Coordinates or FASTA files as single or batch input. CpG-IDs or Genome Coordinates have to be slash ('/') separated for batch input.

FASTA files: FASTA file input accepts standard FATSA files with one or more sequences. The genomic location of the sequence in the fasta file must be specified in the header (Example: >chr19:43203368-43203369)

Genome Coordinate: Genomic Coordinates of the CpG site of interest can be inserted. NOTE: No comma separation! (Example: chr19:43203368-43203369)

CpG ID: The genome assembly hg19 allows the input of CpG-IDs. (Example: cg04882394)

3.) Bisulfite Convert

Choose this option to bisulfite convert FASTA files first. For already bisulfite converted sequences: Capital Letter T indicates non-conversion derived thymins and lower case t for bisulfite conversion derived tymins.

Warning: Including CpG sites in the primer sequence might introduce a methylation bias in the data and should only be chosen by expert users with great care. When including CpG sites in the primer, the CpG location should be at the 5' end of the primer sequence as this minimizes the preferential amplification of an unmethylated allel.

4.) Primer Size

Choose the optimal, minimal and maximal primer size.

5.) Melting Temperature

Choose the optimal, minimal and maximal primer melting temperature.

6.) Amplicon Size

Choose the minimal and maximal amplicon size.

7.) Number of Primer to design

Choose how many primers to design per inserted genomic region.

7.) Exclude SNPs from the primer sequence

Common SNPs (hg38: dbSNP 151; hg19: dbSNP 151; mm10: dbSNP 142) can be excluded from the primer sequence to prevent a sample specific amplification bias.

9.) Advanced Setting:

Extend the input genomic region

The selected genomic regions can be extended by the chosen number of base pairs in both directions (not for FASTA).

Add Sequencing Adapter

Add sequencing adapter 5' of the primer sequence. Adapter sequences can be specified as Adapter Sequence Fwd and Adapter Sequence Rev if this option is selected.

Specify a target region (not for batch input)

For a single genomic region users can specify the position (either genome coordinate or relative base pair position of the amplicon) which should be covered in the amplicon.

9.) Overview of used Settings

After submitting the primer design job to the AmpliconDesign web server, users can re-evaluate their chosen settings in the Used Settings section.

10.) Evaluating the Results

AmpliconDesign generates a table with automatically designed primers. This table includes the primer sequence for the plus strand (Amplicon 1) and the minus strand (Amplicon 2) of each amplicon. The primer ID is given as the specified input region. Primer sequences for the forward and reverse primer are predicted (whereby t corresponds to bisulfite converted C). The primer begin is given with respect to the amplicon start site, followed by the GC content of the primer, melting temperatures and the total amplicon size. The Region column repeats the input region, and the UCSC coordinates for the forward and reverse primer (left and right) are given, as well as the final amplicon coordinates and the number of CpGs in each amplicon (#CpGs). Users can furthermore order the results based on the number of CpG sites in the final amplicon by clicking on the CpG column header in the results table.

Users can choose a primer pair by clicking on a row of the table. A pop-up window opens with a genome region plot, indicating primer binding sites, SNPs, transcript and CpG locations within the specified region.

Furthermore, two AmpliconDesign offers two addtional quality control mechanisms for each primer par: Users can inspect the underlying bisulfite converted amplicon sequence together with the SNPs, repeat or CpG locations and run an ePCR search to identify possible off-target amplicons (BisAlign primer against BS genome).

Interactive Data

Data can be loaded and interactively explored.

BisAlign Tutorial

Introduction to BisAlign

The alignment of designed primers against a reference genome is an important task to ensure the specificity and binding efficacy. This quality control step becomes even more important when working with bisulfite treated DNA as the conversion of cytosins to uracil (and later thymin after PCR amplification) lowers the complexity of the underling sequence from a 4- to a 3-letter code. Therefore, we have implemented BisAlign, a Bowtie based short sequence alignment tool which allows the in-silico alignment of primer sequences against a bisulfite converted genome.

Tutorial: Primer alignment using BisAlign

1.) Select the required alignment mode

Users can choose between a single primer alignment (BisAlign) or exploration of potential off-target amplicons resulting from a PCR with the respective forward and reverse primer (ePCR).
Select the required mode in the drop-down menu.

2.) Select the required genome build

AmpliconDesign currently supports three common genome builds: GRCh38/hg38, GRCh27/hg19 and GRCm38/mm10.
Select the required genome in the drop-down menu.

2.) Enter the primer sequence

Please enter the primer sequence which you want to align agains the bisulfite converted reference genome.
For ePCR mode you need to enter the respective forward and reverse sequence.

Note: You can just enter sequences with standard bases ACGT. For unknwon bases use N.

3.) Results section

The alignment results are shown as a table sorted based on the “top alignment score”. The strand column indicates on which strand a primer alignment was detected. The alignment position is given as chromosome and start of the binding sequence. The primer sequence or the reverse complement which was used for alignment is repeated in the sequence column. The mismatch to the primer sequence at a potential off-target binding sites is shown in the last column.

For the ePCR mode an additional sequence query has been implemented. When clicking on the respective row in the results table, a pop-up window appears which includes the sequence of the off-target amplicon together with a coloring index marking the mismatched bases and primer sequence.

Analysis Pipeline Tutorial

Introduction to the Analysis Pipeline

The analysis of ampBS-seq data requires adapter trimming, read alignment and methylation calling. Subsequent output files contain the number of methylated and unmethylated reads for each CpG sites. AmpliconDesign offers an out-of-the-box pipeline supporting each of these steps. Moreover, Bismark coverage files which are produced by the pipeline can be uploaded to the AmpliconDesign webserver for an interactive quality control and analysis. The following tutorial guides you through all the necessary steps of a successful ampBS-seq analysis workflow.

Tutorial: Analyzing ampBS-seq results using AmpliconDesign

1.) Install and configure the Snakemake pipeline

The snakemake pipeline can be simpely downloaded from github using a single command in the unix shell:

$ git clone https://github.com/HeyLifeHD/amplicon_bisulfite_seq_pipeline/

Open the config.yaml file in the “amplicon_bisulfite_seq_pipeline” directory using your favorite editor.

Set the paths for your working directory (make sure you have rights to write), Bismark reference genome and files for the R analysis pipeline which can be used optional.

Last create an folder called 01_fastq_raw in the working directory in which you move the .fastq files to be processed.

$ mkdir 01_fastq_raw

2.) Run the Snakemake pipeline

All you need to execute this pipeline is to install Snakemake via the Conda package manager (for more information please refer to the Snakemake documentation). Software needed by this workflow is automatically deployed into isolated environments by Snakemake.

$ conda install -c conda-forge -c bioconda snakemake

Test your configuration by performing a dry-run via:

$ snakemake --use-conda -n

Execute the workflow locally via

$ snakemake --use-conda all

3.) Upload files to the AmpliconDesign web server

Sample meta data, analyzed regions and Bismark coverage (methylation information) files can be uploaded to the webserver. Please make sure that your files fulfill the following criteria:
Sample Sheet: The sample sheet must contain the exact file names to which the meta data is assigned in a first column. Furthermore a column “UPN” specifying a unique identifier for each sample is required.

Regions: Genomic coordinates tab-separated in a common bed-file format (Chromosome, Start, End).

Bismark coverage files (.cov): Exact Bismark coverage-file (.cov) output (Chromosome, Start, End, Methylation Percentage, Methylated Reads, Unmethylated Reads). Multiple files can be uploaded. The file names have to match to the first column in the meta data file.

Coverage Cut-Off: CpG sites with less reads will be excluded from the downstream analysis.

4.) Overview of the uploaded data

The Overview tab shows information about the uploaded data. The coverage information for each sample is calculated in given in a table. “Covered CpG” counts the number of CpG sites retrieved from the Bismark .cov-files.

5.) Quality Control

The Quality Control and QC Filtered tabs allow information about the data quality. In the QC Filtered data, the coverage cut-off has been applied whereas Quality Control gives a general overview without filtering.

6.) PCA

A principal component analysis (PCA) can be computed. Users can choose which meta data column should be used for coloring of the samples in the two-dimensional space. Common CpG sites after coverage filtering are used for PCA.

7.) Heatmap

Common CpG sites after coverage filtering can be visualized in a heatmap. Users can choose whether to cluster rows or column, show column and row names, annotate meta data as column labels or change the colore scheme. Clustering is performed using hierarchical clustering (Method: Complete, Distance: Euclidean).

8.) Region

The methylation levels for all CpG sites in an uploaded region can be inspected. Not covered CpG sites are marked in grey. For this view, CpG sites are ordered based on their genomic location. The coverage filter can be applied optionally.