inst/doc/gQTLstats.R

In gQTLstats: gQTLstats: computationally efficient analysis for eQTL and allied studies

## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()

## ----setup,echo=FALSE---------------------------------------------------------
suppressPackageStartupMessages({
library(SummarizedExperiment)
library(Homo.sapiens)
library(org.Hs.eg.db)
library(geuvStore2)
library(gQTLBase)
library(gQTLstats)
})

## ----contset------------------------------------------------------------------
library(geuvStore2)
library(gQTLBase)
library(gQTLstats)
library(parallel)
nco = detectCores()
library(doParallel)
registerDoSEQ()
if (.Platform$OS.type != "windows") {
  registerDoParallel(cores=max(c(1, floor(nco/2))))
}
prst = makeGeuvStore2() 

## ----getqs, cache=TRUE--------------------------------------------------------
qassoc = storeToQuantiles(prst, field="chisq",
    probs=c(seq(0,.999,.001), 1-(c(1e-4,1e-5,1e-6))))
tail(qassoc)

## ----gethist, cache=TRUE------------------------------------------------------
hh = storeToHist( prst, breaks= c(0,qassoc,1e9) )
tail(hh$counts)

## ----twoFDRs, cache=TRUE------------------------------------------------------
rawFDR = storeToFDR(prst, 
   xprobs=c(seq(.05,.95,.05),.975,.990,.995,.9975,.999, 
   .9995, .9999, .99999) )

## ----makefilt-----------------------------------------------------------------
dmfilt = function(x)  # define the filtering function
     x[ which(x$MAF >= 0.05 & x$mindist <= 500000) ] 

## ----runfilt, cache=TRUE------------------------------------------------------
filtFDR = storeToFDR(prst, 
   xprobs=c(seq(.05,.95,.05),.975,.990,.995,.9975,.999, 
   .9995, .9999, .99999), filter = dmfilt )

## ----lktails------------------------------------------------------------------
rawFDR
filtFDR

## ----showfd2, plot=TRUE-------------------------------------------------------
rawtab = getTab(rawFDR)
filttab = getTab(filtFDR)
 plot(rawtab[-(1:10),"assoc"], 
      -log10(rawtab[-(1:10),"fdr"]+1e-6), log="x", axes=FALSE,
  xlab="Observed association", ylab="-log10 plugin FDR")
 axis(1, at=c(seq(0,10,1),100,200))
 axis(2)
 points(filttab[-(1:10),1], -log10(filttab[-(1:10),2]+1e-6), pch=2)
 legend(1, 5, pch=c(1,2), legend=c("all loci", "MAF >= 0.05 & dist <= 500k"))

## ----shobfu,cache=FALSE,eval=FALSE--------------------------------------------
#  fdbp = storeToFDRByProbe( prst, xprobs=c(seq(.025,.975,.025),.99))
#  tail(getTab(fdbp),5)

## ----shobpf,cache=FALSE,eval=FALSE--------------------------------------------
#  fdAtM05bp = storeToFDRByProbe( prst, filter=function(x) x[which(x$MAF > .05)],
#     xprobs=c(seq(.025,.975,.025),.99))
#  tail(getTab(fdAtM05bp),5)

## ----lkgam,fig=TRUE-----------------------------------------------------------
library(mgcv)
rawFDR = setFDRfunc(rawFDR)
filtFDR = setFDRfunc(filtFDR)
par(mfrow=c(2,2))
txsPlot(rawFDR)
txsPlot(filtFDR)
directPlot(rawFDR)
directPlot(filtFDR)

## ----doenums,cache=TRUE-------------------------------------------------------
rawEnum = enumerateByFDR(prst, rawFDR, threshold=.05) 
rawEnum[order(rawEnum$chisq,decreasing=TRUE)[1:3]]
length(rawEnum)

## ----lkmd---------------------------------------------------------------------
names(metadata(rawEnum))

## ----dofenum,cache=TRUE-------------------------------------------------------
filtEnum = enumerateByFDR(prst, filtFDR, threshold=.05,
   filter=dmfilt) 
filtEnum[order(filtEnum$chisq,decreasing=TRUE)[1:3]]
length(filtEnum)

## ----dosens,fig=TRUE----------------------------------------------------------
data(sensByProbe) # see example(senstab) for construction approach
tab = senstab( sensByProbe )
plot(tab)

## ----counts,cache=TRUE--------------------------------------------------------
flens = storeApply( prst, function(x) {
    length(x[ which(x$MAF >= .08 & x$mindist <= 25000), ] )
})

## ----lklen--------------------------------------------------------------------
sum(unlist(flens))

## ----bindback-----------------------------------------------------------------
library(geuvPack)
data(geuFPKM)
basic = mcols(rowRanges(geuFPKM))[, c("gene_id", "gene_status", "gene_type",
    "gene_name")]
rownames(basic) = basic$gene_id
extr = basic[ filtEnum$probeid, ]
mcols(filtEnum) = cbind(mcols(filtEnum), extr)
stopifnot(all.equal(filtEnum$probeid, filtEnum$gene_id))
filtEnum[1:3]

## ----lkscores,fig=TRUE--------------------------------------------------------
data(hmm878)
library(geuvStore2)
prst = makeGeuvStore2() 
myg = "ENSG00000183814.10" # LIN9
data(filtFDR)
library(ggplot2)
manhWngr( store = prst, probeid = myg, sym="LIN9",
  fdrsupp=filtFDR, namedGR=hmm878 )

## ----dovaranno,cache=TRUE,eval=FALSE------------------------------------------
#  suppressPackageStartupMessages({
#  library(VariantAnnotation)
#  library(TxDb.Hsapiens.UCSC.hg19.knownGene)
#  })
#  txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
#  seqlevelsStyle(filtEnum) = "UCSC"
#  #seqinfo(filtEnum) = seqinfo(txdb)
#  seqlengths(filtEnum)[paste0("chr", c(1:22,"M"))] =
#   seqlengths(txdb)[paste0("chr", c(1:22,"M"))]
#  filtEnum = keepStandardChromosomes(filtEnum)
#  suppressWarnings({
#  allv = locateVariants(filtEnum, txdb, AllVariants()) # multiple recs per eQTL
#  })
#  table(allv$LOCATION)
#  hits = findOverlaps( filtEnum, allv )
#  filtEex = filtEnum[ queryHits(hits) ]
#  mcols(filtEex) = cbind(mcols(filtEex), mcols(allv[subjectHits(hits)])[,1:7])
#  filtEex[1:3]

## ----lkall--------------------------------------------------------------------
args(AllAssoc)

## ----demoit-------------------------------------------------------------------
require(GenomeInfoDb)
require(geuvPack)
require(Rsamtools)
data(geuFPKM)  # get a ranged summarized expt
lgeu = geuFPKM[ which(seqnames(geuFPKM)=="chr20"), ] # limit to chr20
seqlevelsStyle(lgeu) = "NCBI"
tf20 = TabixFile(system.file("vcf/c20exch.vcf.gz", package="gQTLstats"))
if (require(VariantAnnotation)) scanVcfHeader(tf20)
set.seed(1234)
mysr = GRanges("20", IRanges(33.099e6, 33.52e6))
lita = AllAssoc(geuFPKM[1:10,], tf20, mysr)
names(mcols(lita))

## ----dem2---------------------------------------------------------------------
litb = AllAssoc(geuFPKM[11:20,], tf20, mysr)
litc = AllAssoc(geuFPKM[21:30,], tf20, mysr)

## ----docoll-------------------------------------------------------------------
buf = gQTLstats:::collapseToBuf(lita, litb, frag="_obs")
buf
buf = gQTLstats:::collapseToBuf(buf, litc, frag="_obs")
buf

## ----doperm-------------------------------------------------------------------
pbuf = gQTLstats:::collapseToBuf(lita, litb, frag="_permScore_1")
pbuf = gQTLstats:::collapseToBuf(pbuf, litc, frag="_permScore_1")
pbuf

## ----dof,fig=TRUE-------------------------------------------------------------
plot(density(buf$scorebuf[,1]))
lines(density(pbuf$scorebuf[,1]), lty=2)