## ----setup, include=FALSE, cache=FALSE----------------------------------------
knitr::opts_chunk$set(echo = TRUE)
#knitr::opts_knit$set(root.dir="")

## ----installation, echo=TRUE, eval=FALSE--------------------------------------
#  if(!requireNamespace("BiocManager")){
#      install.packages("BiocManager")
#  }
#  if(!requireNamespace("MAGAR")){
#      BiocManager::install("MAGAR")
#  }

## ----loading, eval=TRUE-------------------------------------------------------
suppressPackageStartupMessages(library(MAGAR))  

## ----options, eval=TRUE-------------------------------------------------------
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)

## ----plink_folder, eval=FALSE-------------------------------------------------
#  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())

## ----imputation options, eval=TRUE--------------------------------------------
qtlSetOption(
  impute.geno.data=TRUE,
  imputation.reference.panel="apps@hrc-r1.1",
  imputation.phasing.method="shapeit",
  imputation.population="eur"
)

## ----imputed_import,eval=FALSE------------------------------------------------
#  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())

## ----call_methQTL, eval=TRUE--------------------------------------------------
imp.data <- loadMethQTLInput(system.file("extdata","reduced_methQTL",package="MAGAR"))
qtlSetOption(standard.deviation.gauss=100,
    cluster.cor.threshold=0.75)
meth.qtl.res <- doMethQTL(imp.data,default.options=FALSE,p.val.cutoff=0.05)

## ----show_results, eval=TRUE--------------------------------------------------
result.table <- getResult(meth.qtl.res)
head(result.table)
anno.meth <- getAnno(meth.qtl.res,"meth")
head(anno.meth)
anno.geno <- getAnno(meth.qtl.res,"geno")
head(anno.geno)

## ----plots,eval=TRUE, include=TRUE--------------------------------------------
result.table <- result.table[order(result.table$P.value,decreasing=FALSE),]
qtlPlotSNPCpGInteraction(imp.data,result.table$CpG[1],result.table$SNP[1])
qtlDistanceScatterplot(meth.qtl.res)

## ----interpretation, eval=TRUE, include=TRUE, warning=FALSE-------------------
res <- qtlBaseSubstitutionEnrichment(meth.qtl.res)

## ----list, eval=TRUE, warning=FALSE-------------------------------------------
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())
qtlUpsetPlot(meth.qtl.list,type = "cor.block")
spec.first <- getSpecificQTL(meth.qtl.list$First,meth.qtl.list[-1])

## ----cluster, eval=FALSE------------------------------------------------------
#  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)