## ----style, echo = FALSE, results = 'asis'------------------------------- suppressMessages({ suppressWarnings({ suppressPackageStartupMessages({ BiocStyle::markdown() }) }) }) ## ----bib, echo = FALSE, results = 'hide'--------------------------------- suppressMessages({ suppressWarnings({ suppressPackageStartupMessages({ library(knitcitations) library(bibtex) library(yriMulti) library(dsQTL) library(geuvPack) library(erma) library(VariantAnnotation) library(gQTLstats) library(Homo.sapiens) }) }) }) allbib = read.bibtex("allbib.bib") ## ----fakeload1,echo=FALSE,results="hide"--------------------------------- if (!exists("geuFPKM")) data(geuFPKM) ## ----lkexpr,eval=TRUE---------------------------------------------------- library(geuvPack) data(geuFPKM) ## ----do1----------------------------------------------------------------- geuFPKM ## ----fakeload2,echo=FALSE,results="hide"--------------------------------- if (!exists("banovichSE")) data(banovichSE) ## ----lkmul,eval=TRUE----------------------------------------------------- library(yriMulti) data(banovichSE) ## ----do2----------------------------------------------------------------- banovichSE ## ----lkdhsf,echo=FALSE,results="hide"------------------------------------ if (!exists("DHStop5_hg19")) data(DHStop5_hg19) ## ----lkdhs,eval=TRUE----------------------------------------------------- library(dsQTL) if (!exists("DHStop5_hg19")) data(DHStop5_hg19) ## ----lkddd--------------------------------------------------------------- DHStop5_hg19 ## ----lkgeno,cache=TRUE--------------------------------------------------- litvcf = readVcf(gtpath(20), param=ScanVcfParam(which=GRanges("20", IRanges(3.7e7,3.701e7))), genome="hg19") ## ----outca--------------------------------------------------------------- litvcf length(colnames(litvcf)) length(intersect(colnames(litvcf), colnames(banovichSE))) length(intersect(colnames(litvcf), colnames(geuFPKM))) length(intersect(colnames(litvcf), colnames(DHStop5_hg19))) ## ----lkmex,cache=TRUE---------------------------------------------------- m1 = mexGR(banovichSE, geuFPKM, symbol="ORMDL3") m1 mcols(m1)[1:4,1:4] table(mcols(m1)$type) ## ----lkbi,cache=TRUE----------------------------------------------------- b1 = bindelms(geuFPKM, banovichSE, symbol="BRCA2", ytx=log, gradius=20000) b1 mcols(b1)[1:3,] summary(mcols(b1)$t) mintind = which.min(mcols(b1)$t) mincpg = names(b1)[mintind] mincpg ## ----lkevm,fig=TRUE------------------------------------------------------ plotEvM(b1) ## ----lkbcma, eval=TRUE--------------------------------------------------- library(MultiAssayExperiment) myobs = list(geuvRNAseq=geuFPKM, yri450k=banovichSE, yriDHS=DHStop5_hg19) cold = colData(geuFPKM) suppressWarnings({ library(MultiAssayExperiment) mm = MultiAssayExperiment(myobs, as.data.frame(cold)) }) mm ## ----lkbrc, eval=TRUE---------------------------------------------------- library(erma) brr = range(genemodel("BRCA2")) brr ## ----dorsub, eval=TRUE--------------------------------------------------- .subsetByRanges = function(ma, r) { subsetByRow(ma,r) } newmm = .subsetByRanges(mm, brr+20000) newmm ## ----makeco-------------------------------------------------------------- library(doParallel) registerDoSEQ() allLM_pw = function(fmla, mae, xtx=force, ytx=force) { # # formula specifies dependent and independent assays # form all regressions of ytx(dep) on xtx(indep) for all # pairs of dependent and independent variables defined by # assays for samples held in common # lf = as.list(fmla) nms = lapply(lf, as.character) yel = experiments(mae)[[nms[[2]]]] xel = experiments(mae)[[nms[[3]]]] sy = colnames(yel) sx = colnames(xel) sb = intersect(sy,sx) yel = yel[,sb] xel = xel[,sb] vdf = as.matrix(expand.grid( rownames(yel), rownames(xel), stringsAsFactors=FALSE )) allf = apply(vdf, 1, function(x) as.formula(paste(x, collapse="~"))) alllm = foreach (i = 1:length(allf)) %dopar% { df = data.frame(ytx(assay(yel)[vdf[i,1],]), xtx(assay(xel)[vdf[i,2],])) names(df) = vdf[i,] lm(allf[[i]], data=df) } names(alllm) = apply(vdf,1,function(x) paste(x, collapse="~")) allts = lapply(alllm, function(x) summary(x)$coef[2, "t value"]) list(mods=alllm, tslopes=allts) } pwplot = function(fmla1, fmla2, mae, ytx=force, xtx=force, ...) { # # use fmla1 with assays as components to identify # two assays to regard as sources of y and x # fmla2 indicates which features to plot # lf = as.list(fmla1) nms = lapply(lf, as.character) yel = experiments(mae)[[nms[[2]]]] xel = experiments(mae)[[nms[[3]]]] sy = colnames(yel) sx = colnames(xel) sb = intersect(sy,sx) yel = yel[,sb] xel = xel[,sb] lf2 = lapply(as.list(fmla2), as.character) ndf = data.frame( ytx(assay(yel)[ lf2[[2]], ]), xtx(assay(xel)[ lf2[[3]], ]) ) names(ndf) = c(lf2[[2]], lf2[[3]]) plot(fmla2, ndf, ...) } ## ----lkpwl, eval=TRUE---------------------------------------------------- pp = allLM_pw(geuvRNAseq~yri450k, newmm, ytx=log) names(pp) summary(pp[[1]][[1]]) which.min(unlist(pp[[2]])) # not BRCA2 but FRY ## ----lkf,fig=TRUE, eval=TRUE--------------------------------------------- pwplot(geuvRNAseq~yri450k, ENSG00000139618.9~cg20073910, newmm, ytx=log) ## ----mkvs, eval=TRUE----------------------------------------------------- library(gQTLstats) library(VariantAnnotation) library(GenomicFiles) pa = paths1kg(paste0("chr", c(21:22))) #,"Y"))) sn = vcfSamples(scanVcfHeader(TabixFile(pa[1]))) library(Homo.sapiens) # necessary? stopifnot(requireNamespace("GenomeInfoDb")) ob = RangedVcfStack(VcfStack(pa, seqinfo(Homo.sapiens))) cd = DataFrame(id1kg=sn) rownames(cd) = sn colData(ob) = cd ## ----mkrvs, eval=TRUE---------------------------------------------------- myr = GRanges("chr22", IRanges(20e6,20.01e6)) rowRanges(ob) = myr colData(ob) = DataFrame(colData(ob), zz=runif(nrow(colData(ob)))) hasInternetConnectivity = function() !is.null(nsl("www.r-project.org")) #if (hasInternetConnectivity()) lka = assay(ob) myobs = list(geuvRNAseq=geuFPKM, yri450k=banovichSE, yriDHS=DHStop5_hg19, yriGeno=ob) suppressWarnings({ mm = MultiAssayExperiment(myobs) }) mm ## ----results='asis',echo=FALSE------------------------------------------- bibliography() #style="markdown")