Required external software tools

MAGAR depends on a funtional installation of PLINK for handling genotyping data. Follow the installation instructions available here and specify the path to the installation for the package qtlSetOption(plink.path="path_to_plink").

Additionally, depending on the type of the analysis, more external software tools are required (). More specifically, if you want to perform imputation bgzip and tabix from the htslib package are needed, as well as vcftools for sorting VCF files. Lastly, if methQTL calling is to be performed through fastQTL, the software tool needs to be installed and the location of the executable specified to MAGAR using the qtlSetOption function with the value bgzip.path, tabix.path, vcftools.path, fast.qtl.path.

DNA methylation data (microarrays)

MAGAR utilizes the widely used RnBeads software package for DNA methylation data import. Thus, MAGAR supports the various input options available in RnBeads, including a direct download from the Gene Expression Omnibus (GEO), IDAT, and BED files. For further options, we refer to the RnBeads vignette and documentation. In addition to the raw methylation data, a sample annotation sheet specifying the samples to be analyzed needs to be provided. The sheet contains a line for each sample and looks as follows:

SampleID,age,sex,barcode
Sample_1,14,f,209054857842_R01C01
Sample_2,42,f,209054857842_R02C01
Sample_3,45,m,209054857842_R03C01

For further details on the import process, we refer to the RnBeads vignette. Most importantly, analysis options need to be specified for the import and preprocessing modules of RnBeads. MAGAR provides a default setting, which is available in extdata/rnbeads_options.xml. You can use this file as a template for your own setting and then specify it to the package:

opts <- rnb.xml2options(system.file("extdata","rnbeads_options.xml",package="MAGAR"))
rnb.options(identifiers.column="geo_accession",
    import.idat.platform="probes450")
xml.fi <- file.path(getwd(),"rnbeads_options.xml")
cat(rnb.options2xml(),file=xml.fi)
qtlSetOption(rnbeads.options = xml.fi)

Genotyping data

geno.dir <- "plink_data"
rnb.set <- load.rnb.set("rnb_set_dir")
s.anno <- "sample_annotation.csv"
data.loc <- list(idat.dir=rnb.set,geno.dir=plink.dir)
qtlSetOption(geno.data.type="plink",
    meth.data.type="rnb.set")
meth.qtl <- doImport(data.location=data.loc,
    s.anno=s.anno,
    s.id.col="SampleID",
    out.folder=getwd())

SampleID,age,sex,barcode,GenoSentrixPosition
Sample_1,14,f,209054857842_R01C01,9701756058_R05C01
Sample_2,42,f,209054857842_R02C01,9701756058_R07C01
Sample_3,45,m,209054857842_R03C01,9742011016_R04C01

qtlSetOption(
  impute.geno.data=TRUE,
  imputation.reference.panel="[email protected]",
  imputation.phasing.method="shapeit",
  imputation.population="eur"
)

Perform data import

The doImport function requires the paths to the respective genotyping and DNA methylation data, as well as a sample annotation sheet as discussed earlier. In this vignette, we will describe the import of DNA methylation data in IDAT format and genotyping data as PLINK files. First, you’ll have to specify the paths to the corresponding IDAT and plink files. Additionally, you have to specify the sample identifier column in the sample annotation sheet that determines the samples in both the genotyping and DNA methylation data. For larger files, we recommend to activate the option to store large matrices on disk rather than in main memory (hdf5dump).

For imputed data, no further processing is performed on the genotyping data and the dosage values are used as they are:

idat.dir <- "idat_dir"
geno.dir <- "geno_dir"
anno.sheet <- "sample_annotation.tsv"
qtlSetOption(hdf5dump=TRUE)
imp.data <- doImport(data.location = c(idat.dir=idat.dir,geno.dir=geno.dir),
                    s.anno = anno.sheet,
                    s.id.col = "ID",
                    tab.sep = "\t",
                    out.folder = getwd())

Please note that the recode.allele.frequencies option specifies, if, according to the cohort analyzed, SNP reference and alternative allele are to be recoded according to the allele frequencies in the available samples. Alternatively, a path to a local version of dbSNP (Sherry et al. 2001) can be provided through db.snp.ref, and reference/alternative allele information will be automatically parsed from the database. This is especially crucial, if imputation is to be performed, since the Michigan Imputation Server is sensitive to reference mismatches.

Compute CpG correlation blocks

Since neighboring CpGs are often highly correlated, using each CpG independently as a potential methQTL candidate leads to many redundant results. We thus aimed to approximate DNA methylation haplotypes by determining highly correlated CpGs in close vicinity. The procedure itself is split into six steps, and is performed for each chromosome independently:

1. Compute the (Pearson) correlation matrix between all CpGs (futher correlation types available in option correlation.type)
2. Construct the distance matrix from the correlation matrix
3. Discard all interactions with a correlation lower than a given threshold (option: cluster.cor.threshold)
4. Weight the distance according to the genomic distance between the two CpGs with a Gaussian (option: standard.deviation.gauss). Higher values for the standard deviation lead to a lower penalty on distal CpGs, thus the clusters will become larger.
5. Discard all interactions ranging longer than the option absolute.distance.cutoff
6. Compute the Louvain clustering on the undirected, weighted graph induced by the distance matrix

This will return a clustering according to the correlation structure between neighboring CpGs that we will later use for methQTL calling. Note that we used simultation experiments to determine the parameters for each data type individually. They will be automatically loaded for the dataset that is used and are:

• 450k: cluster.cor.threshold=0.2, standard.deviation.gauss=5,000, absolute.distance.cutoff=500,000
• EPIC: cluster.cor.threshold=0.2, standard.deviation.gauss=3,000, absolute.distance.cutoff=500,000
• RRBS/WGBS: cluster.cor.threshold=0.2, standard.deviation.gauss=250, absolute.distance.cutoff=500,000

Call methQTL per correlation block

From the list of correlation blocks, MAGAR computes methQTL interactions with all SNPs on the same chromosome. The process is split into three steps:

1. Compute a representative CpG (tag-CpG) per correlation block, as specified with the option representative.cpg.computation (default: row.medians).
2. Discard all SNPs that are further than absolute.distance.cutoff (default: 1,000,000) away from the representative CpG
3. Call methQTL by using linear models. Multiple options of methQTL calling are available and can be selected via the option linear.model.type (default: classical.linear). Alternatively, fastQTL can be set as an option for meth.qtl.type. This will tell the package to use the fastQTL software (Ongen et al. 2016).

The meth.qtl.type tells, how a methQTL interaction is defined and provides three options, in addition to the already mentioned fastQTL:

1. oneVSall: A CpG can only be influenced by one SNP. We choose the one with the lowest p-value.
2. twoVSall: A CpG can both positively and negatively be influenced by two independent SNPs. The package will output those fulfilling the p-value cutoff.
3. allVSall: For each CpG, all SNPs showing a p-value lower than the p-value cutoff will be returned.

In the methQTL calling process, potential covariates can be specified using the option sel.covariates. We recommend to include at least age and sex as covariates, as they have a strong influence on the DNA methylation pattern.

How to use methQTLResult

The above procedure will create an object of class methQTLResult, which contains the methQTL that are called in the previous step. To get a table of all the methQTL, you need to extract the information from the object. In the majority of the function calls below, there is the option type, which takes on the values: * ‘SNP’: To characterize the SNPs that influence any DNA methylation state * ‘CpG’: To characterize the representative CpGs per correlation block that are influences by any SNP genotype * ‘cor.block’: To characterize all CpGs, which are part of a correlation block, whose representative CpG is influenced by any genotype

Furthermore, you can obtain genomic annotations for both the CpGs and the SNPs involved in the methQTL interactions:

result.table <- getResult(meth.qtl.res)
head(result.table)

##           CpG          SNP         Beta     SE.Beta     P.value Chromosome
## 18 cg00014996    rs9948869 -0.009183646 0.002665375 0.001087206      chr18
## 14 cg01750507 18:12842:G:A -0.009789160 0.004595578 0.037568569      chr18
## 20 cg01848352   rs62075713 -0.018601404 0.008481772 0.032469421      chr18
## 43 cg02931136   rs78288743  0.021477662 0.007926956 0.008926114      chr18
## 35 cg03122836   rs78288743  0.023898710 0.009456236 0.014344688      chr18
## 40 cg03484476   rs78288743 -0.047466160 0.019065425 0.015782950      chr18
##    Position.CpG Position.SNP Distance p.val.adj.fdr
## 18       158348        73323    85025             1
## 14       158229        12842   145387             1
## 20       159265        48740   110525             1
## 43       475647        65068   410579             1
## 35       332912        65068   267844             1
## 40       457239        65068   392171             1

anno.meth <- getAnno(meth.qtl.res,"meth")
head(anno.meth)

##            Chromosome Start   End Strand Strand.1 AddressA AddressB Design
## cg00703566      chr18 10862 10863      +        + 40648584       NA     II
## cg21164455      chr18 10870 10871      +        + 49733949       NA     II
## cg15244183      chr18 11143 11144      +        + 32645179 38797212      I
## cg20937418      chr18 11145 11146      +        + 12701867       NA     II
## cg21218768      chr18 11155 11156      +        + 43687103       NA     II
## cg18836174      chr18 11724 11725      -        - 54678938 28656821      I
##            Color Context Random HumanMethylation27 HumanMethylation450
## cg00703566  Both      CG  FALSE              FALSE                  NA
## cg21164455  Both      CG  FALSE              FALSE                  NA
## cg15244183   Grn      CG  FALSE              FALSE                  NA
## cg20937418  Both      CG  FALSE              FALSE                  NA
## cg21218768  Both      CG  FALSE              FALSE                  NA
## cg18836174   Grn      CG  FALSE              FALSE                TRUE
##            Mismatches.A Mismatches.B CGI.Relation CpG GC SNPs.3 SNPs.5
## cg00703566            0            0  North Shore   9 62      0      0
## cg21164455            0            0  North Shore   8 60      0      0
## cg15244183            0            0  North Shore  15 74      0      0
## cg20937418            0            0  North Shore  15 74      0      0
## cg21218768            0            0  North Shore  13 70      0      0
## cg18836174            0            0       Island   6 68      0      0
##            SNPs.Full Cross.reactive
## cg00703566         0              0
## cg21164455         0              0
## cg15244183         0              0
## cg20937418         0              0
## cg21218768         0              0
## cg18836174         0              0

anno.geno <- getAnno(meth.qtl.res,"geno")
head(anno.geno)

##              Chromosome Start cM Allele.1 Allele.2 Allele.1.Freq Allele.2.Freq
## 18:11149:C:T      chr18 11149 NA        C        T     0.9310345    0.06896552
## 18:11543:G:A      chr18 11543 NA        G        A     0.8965517    0.10344828
## 18:11860:G:A      chr18 11860 NA        G        A     0.6551724    0.34482759
## 18:12842:G:A      chr18 12842 NA        G        A     0.8965517    0.10344828
## 18:13034:G:A      chr18 13034 NA        G        A     0.5172414    0.48275862
## 18:13413:C:A      chr18 13413 NA        C        A     0.8793103    0.12068966

For more detailed information about the output, also see the function getResultsGWASMap.

Plots

To visualize methQTL, MAGAR provides some plotting functions. Most functions return an object of type ggplot, which can be subsequently stored or viewed. Either all methQTL can be simultaneously visualized in a single plot, or a specific methQTL can be visualized:

result.table <- result.table[order(result.table$P.value,decreasing=FALSE),]
qtlPlotSNPCpGInteraction(imp.data,result.table$CpG[1],result.table$SNP[1])

## `geom_smooth()` using formula = 'y ~ x'

qtlDistanceScatterplot(meth.qtl.res)

Interpretation functions

The package provides a bunch of interpretation functions to characterize the detected methQTLs. This includes LOLA enrichment analysis(Sheffield and Bock 2016) (qtlLOLAEnrichment), genomic annotation enrichment based on putative regulatory elements defined by the Ensembl Regulatory Build(Zerbino et al. 2015) (qtlAnnotationEnrichment), enrichment analysis of different base substitutions in SNPs (qtlBaseSubstitutionEnrichment), or TFBS motif enrichment using TFBSTools. Enrichment is compared for the methQTLs that are available in the provided methQTLResult (for a single input), or to the overlapping QTLs for a list of methQTLResult. The background of the enrichment is defined as all the SNPs/CpGs that have been used as input to the methQTL calling.

res <- qtlBaseSubstitutionEnrichment(meth.qtl.res)

Lists of methQTL results

Most of the functions discussed above either support a single methQTLResult as input, or a list of such objects. In case a list is specified, the functions will typically overlap the methQTLs found and compare those with all SNPs/CpGs that have been used for methQTL calling. Additionally, there are functions that particularly work on a list of methQTLResult objects and that perform overlapping, or determine the methQTLs specific to a dataset.

meth.qtl.res.2 <- loadMethQTLResult(system.file("extdata","MethQTLResult_chr18",package="MAGAR"))
meth.qtl.list <- list(First=meth.qtl.res,Second=meth.qtl.res.2)
qtlVennPlot(meth.qtl.list,out.folder=getwd())

## Loading required namespace: VennDiagram

qtlUpsetPlot(meth.qtl.list,type = "cor.block")

## 2024-12-05 21:51:15     2.4  STATUS STARTED Constructing universe
## 2024-12-05 21:51:15     2.4  STATUS     STARTED Obtaining correlation blocks for object 1
## 2024-12-05 21:51:15     2.4  STATUS     COMPLETED Obtaining correlation blocks for object 1
## 2024-12-05 21:51:15     2.4  STATUS     STARTED Obtaining correlation blocks for object 2
## 2024-12-05 21:51:15     2.5  STATUS     COMPLETED Obtaining correlation blocks for object 2
## 2024-12-05 21:51:15     2.5  STATUS     STARTED Overlapping
## 2024-12-05 21:51:15     2.5  STATUS         STARTED Obtaining correlation blocks for object 1
## 2024-12-05 21:51:15     2.5  STATUS         COMPLETED Obtaining correlation blocks for object 1
## 2024-12-05 21:51:15     2.5  STATUS         STARTED Obtaining correlation blocks for object 2
## 2024-12-05 21:51:15     2.5  STATUS         COMPLETED Obtaining correlation blocks for object 2
## 2024-12-05 21:51:15     2.5  STATUS     COMPLETED Overlapping

spec.first <- getSpecificQTL(meth.qtl.list$First,meth.qtl.list[-1])

MAGAR options

MAGAR is a flexible software packages that allows for multiple flavors of methQTL analyses. Here, we present the different options that describe the analysis in a table and also discuss possible problems with choosing non-default options. Most options have reasonable default values that have been determined using a simulation experiment. For further information, see the documentation of qtlSetOption

Employ MAGAR on a scientific compute cluster

MAGAR can automatically distribute jobs across a high performance compute cluster, which has been setup using the Sun Grid Engine (SGE) technology. You can pass the option cluster.submit to doMethQTL and thus activate the cluster submission. Note that you’ll also have to specify a path to an executable Rscript and potentially specify resource requirements using the option setting cluster.config.

qtlSetOption(cluster.config = c(h_vmem="60G",mem_free="20G"))
qtlSetOption(rscript.path = "/usr/bin/Rscript")
meth.qtl.res <- doMethQTL(meth.qtl = imp.data,
                        cluster.submit = T)