identifying and plotting an meQTL

This chunk illustrates use of MultiAssayExperiment to combine local methylation and remote VCF to find and visualize an meQTL

library(MultiAssayExperiment) # need https://github.com/vjcitn/MultiAssayExperiment/issues/171 resolved
library(ldblock) # get stack from cloud
st = stack1kg()
library(yriMulti) # get 450k data on YRI
data(banovichSE)
library(erma)  # have a look at gene address for  fig 3 of
genemodel("HLA-DQB1") # http://dx.plos.org/10.1371/journal.pgen.1004663
genome(banovichSE) = "b37" # the publication uses b36 (quietly) but we won't
seqlevelsStyle(banovichSE) = "NCBI"
ee = ExperimentList(list(meth=banovichSE, snp=st))
mee = MultiAssayExperiment(ee, pData=colData(banovichSE))
foc = subsetByRow(mee, rowRanges(banovichSE)["cg04793911"]) # one locus
library(gQTLstats)
assoc = cisAssoc( experiments(foc)[["meth"]],
   TabixFile( files(experiments(foc)[["snp"]] ) ) )
eqBox2( "cg04793911", experiments(foc)[["meth"]],
  TabixFile( files(experiments(foc)[["snp"]] ) ),
  assoc[ which.max(assoc$chisq) ] )