gQTLstats: computationally efficient analysis and interpretation of large eQTL, mQTL, etc. archives
Vincent J. Carey, stvjc at channing.harvard.edu
Jun 2015
Introduction
The software in this package aims to support refinements and functional interpretation of members of a collection of association statistics on a family of feature \(\times\) genome hypotheses. provide a basis for refinement or functional interpretation.
We take for granted the use of the gQTL* infrastructure for testing and management of test results. We use for examples elements of the geuvPack and geuvStore2 packages.
Basic infrastructure for statistics on a distributed store of eQTL results
We work with a ciseStore instance based on a small subset of transcriptome-wide cis-eQTL tests for GEUVADIS FPKM data. The overall testing procedure was conducted for all SNP:probe pairs for which SNP minor allele frequency (MAF) is at least 1% and for which the minimum distance between SNP and either boundary of the gene coding region for the probe is at most 1 million bp.
library(geuvStore2)
library(gQTLBase)
library(gQTLstats)
library(parallel)
nco = detectCores()
library(doParallel)## Loading required package: foreach## Loading required package: iteratorsregisterDoSEQ()
if (.Platform$OS.type != "windows") {
registerDoParallel(cores=max(c(1, floor(nco/2))))
}
prst = makeGeuvStore2() Estimation of quantiles of the distribution of observed associations
Quantile estimation is very memory-efficient, based on a temporary ff representation of the vector of all association test results.
qassoc = storeToQuantiles(prst, field="chisq",
probs=c(seq(0,.999,.001), 1-(c(1e-4,1e-5,1e-6))))
tail(qassoc)## 99.7% 99.8% 99.9% 99.99% 99.999% 99.9999%
## 20.82392 29.06245 54.93786 193.82782 220.55134 235.75722Estimation of histogram of the distribution of association scores under permutation
Because we compute fixed breaks, contributions to the overall histogram can be assembled in parallel, with small footprint. This is a tremendous reduction of data.
hh = storeToHist( prst, breaks= c(0,qassoc,1e9) )
tail(hh$counts)## [1] 4333 1933 416 0 0 1Computing FDR from a gqtlStore
FDR computation is post-hoc relative to filtering that need not be specified prior to testing. For illustration, we survey the results in geuvStore2 to obtain FDRs for each SNP:probe pair in two forms. First, we obtain FDR without any filtering. Second, we compute an a FDR for those SNP:probe pairs separated by at most 500kb, and for which the MAF for the SNP is at least 5 per cent.
rawFDR = storeToFDR(prst,
xprobs=c(seq(.05,.95,.05),.975,.990,.995,.9975,.999,
.9995, .9999, .99999) )## counting tests...## counting #NA...## obtaining assoc quantiles...## computing perm_assoc histogram....dmfilt = function(x) # define the filtering function
x[ which(x$MAF >= 0.05 & x$mindist <= 500000) ] filtFDR = storeToFDR(prst,
xprobs=c(seq(.05,.95,.05),.975,.990,.995,.9975,.999,
.9995, .9999, .99999), filter = dmfilt )## counting tests...## counting #NA...## obtaining assoc quantiles...## computing perm_assoc histogram....rawFDR## FDRsupp instance with 27 rows.
## assoc fdr ncalls avg.false
## 5% 0.003666364 0.9503595 5874027 5582437
## 10% 0.014663963 0.9489049 5564867 5280530
## 15% 0.033307623 0.9466335 5255708 4975229
## ...
## assoc fdr ncalls avg.false
## 99.95% 157.9120 5.390964e-05 3091.59300 0.1666667
## 99.99% 193.8278 2.695482e-04 618.31860 0.1666667
## 99.999% 220.5513 0.000000e+00 61.83186 0.0000000
## No interpolating function is available; use 'setFDRfunc'.filtFDR## FDRsupp instance with 27 rows.
## assoc fdr ncalls avg.false
## 5% 0.004052362 0.9461793 2102872 1989694
## 10% 0.016469885 0.9415121 1992194 1875675
## 15% 0.037696938 0.9363941 1881517 1761841
## ...
## assoc fdr ncalls avg.false
## 99.95% 193.8278 0 1106.77450 0
## 99.99% 194.9828 0 221.35490 0
## 99.999% 228.3958 0 22.13549 0
## No interpolating function is available; use 'setFDRfunc'.The filtering leads to a lower FDR for a given strength of association. This is an inspiration for sensitivity analysis. Even with 5 million observations there is an effect of histogram bin selection in summarizing the permutation distribution of association. This can be seen fairly clearly in the wiggliness of the trace over the unfiltered association score:FDR plot.
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"))
We’ll address this below by fitting smooth functions for the score:FDR relationship.
Estimates of FDR at the probe level
The storeToFDRByProbe FDR function examines the maximal association score by gene, for observed and permuted measures.
Good performance of this procedure is obtained by using group_by and summarize utilities of dplyr. Iteration employs foreach.
fdbp = storeToFDRByProbe( prst, xprobs=c(seq(.025,.975,.025),.99))
tail(getTab(fdbp),5)fdAtM05bp = storeToFDRByProbe( prst, filter=function(x) x[which(x$MAF > .05)],
xprobs=c(seq(.025,.975,.025),.99))
tail(getTab(fdAtM05bp),5)Modeling the association-FDR relationship to support efficient variant selection and annotation
Choosing a smooth model for the association:FDR relationship
We’ll focus here on all-pairs analysis, with and without filtering.
Especially in this small example there will be some wiggling or even non-monotonicity in the trace of empirical FDR against association. We want to be able to compute the approximate FDR quickly and with minimal assumptions and pathology. To accomplish this, we will bind an interpolating model to the FDR estimates that we have. Interpolation will be accomplished with scatterplot smoothing in the mgcv framework.
The code that is used to fit the interpolating model is
fdrmod = gam(-log10(fdr+fudge)~s(assoc,bs="tp"), data=...,
subset=assoc<(1.1*maxch))where fudge defaults to 1e-6 and maxch defaults to 30
library(mgcv)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:BBmisc':
##
## collapse## The following object is masked from 'package:IRanges':
##
## collapse## This is mgcv 1.8-22. For overview type 'help("mgcv-package")'.rawFDR = setFDRfunc(rawFDR)
filtFDR = setFDRfunc(filtFDR)
par(mfrow=c(2,2))
txsPlot(rawFDR)
txsPlot(filtFDR)
directPlot(rawFDR)
directPlot(filtFDR)
More work is needed on assessing tolerability of relative error in FDR interpolation.
Enumerating significant cis-eQTL in a store
Recall that dmfilt is a function that obtains the SNP-probe pairs for which SNP has MAF at least five percent and SNP-probe distance at most 500kbp.
We use the FDRsupp instances with ciseStore to list the SNP-probe pairs with FDR lying beneath a given upper bound.
Unfiltered pairs:
rawEnum = enumerateByFDR(prst, rawFDR, threshold=.05)
rawEnum[order(rawEnum$chisq,decreasing=TRUE)[1:3]]## GRanges object with 3 ranges and 15 metadata columns:
## seqnames ranges strand | paramRangeID
## <Rle> <IRanges> <Rle> | <factor>
## 55 14 [106552724, 106552724] * | ENSG00000211968.2
## 15 1 [ 54683925, 54683925] * | ENSG00000231581.1
## 15 1 [ 54685855, 54685855] * | ENSG00000231581.1
## REF ALT chisq permScore_1 permScore_2
## <DNAStringSet> <CharacterList> <numeric> <numeric> <numeric>
## 55 C T 244.3467 0.05015518 0.06192041
## 15 G A 242.4429 3.80053240 0.07174190
## 15 G A 242.4429 3.80053240 0.07174190
## permScore_3 permScore_4 permScore_5 permScore_6 snp
## <numeric> <numeric> <numeric> <numeric> <character>
## 55 1.1187687 0.01283664 0.01086588 0.0003114708 rs587662269
## 15 0.3136501 0.07873665 4.98223625 0.0443687607 rs6621
## 15 0.3136501 0.07873665 4.98223625 0.0443687607 rs33988698
## MAF probeid mindist estFDR
## <numeric> <character> <numeric> <numeric>
## 55 0.004494382 ENSG00000211968.2 525649 0
## 15 0.125842697 ENSG00000231581.1 6256 0
## 15 0.125842697 ENSG00000231581.1 4326 0
## -------
## seqinfo: 86 sequences from hg19 genomelength(rawEnum)## [1] 44750A small quantity of metadata is bound into the resulting GRanges instance.
names(metadata(rawEnum))## [1] "enumCall" "enumSess" "fdrCall"Pairs meeting MAF and distance conditions are obtained with a filter setting to the enumerating function.
filtEnum = enumerateByFDR(prst, filtFDR, threshold=.05,
filter=dmfilt)
filtEnum[order(filtEnum$chisq,decreasing=TRUE)[1:3]]## GRanges object with 3 ranges and 15 metadata columns:
## seqnames ranges strand | paramRangeID
## <Rle> <IRanges> <Rle> | <factor>
## 15 1 [54683925, 54683925] * | ENSG00000231581.1
## 15 1 [54685855, 54685855] * | ENSG00000231581.1
## 15 1 [54683014, 54683014] * | ENSG00000231581.1
## REF ALT chisq permScore_1 permScore_2
## <DNAStringSet> <CharacterList> <numeric> <numeric> <numeric>
## 15 G A 242.4429 3.800532 0.07174190
## 15 G A 242.4429 3.800532 0.07174190
## 15 C G 241.9734 4.225007 0.02231639
## permScore_3 permScore_4 permScore_5 permScore_6 snp MAF
## <numeric> <numeric> <numeric> <numeric> <character> <numeric>
## 15 0.3136501 0.07873665 4.982236 0.04436876 rs6621 0.1258427
## 15 0.3136501 0.07873665 4.982236 0.04436876 rs33988698 0.1258427
## 15 0.2092431 0.02744535 4.326537 0.05104122 rs1410896 0.1224719
## probeid mindist estFDR
## <character> <numeric> <numeric>
## 15 ENSG00000231581.1 6256 0
## 15 ENSG00000231581.1 4326 0
## 15 ENSG00000231581.1 7167 0
## -------
## seqinfo: 86 sequences from hg19 genomelength(filtEnum)## [1] 81837Sensitivity analysis for eQTL enumeration
The yield of an enumeration procedure depends on filtering based on SNP-gene distance and SNP MAF. This can be illustrated as follows, with minimal computational effort owing to the retention of genome-scale permutations and the use of the plug-in FDR algorithm.
data(sensByProbe) # see example(senstab) for construction approach
tab = senstab( sensByProbe )
plot(tab)
If we wish to maximize the yield of eQTL enumeration at FDR at most 0.05, we can apply a filter to the store.
flens = storeApply( prst, function(x) {
length(x[ which(x$MAF >= .08 & x$mindist <= 25000), ] )
})sum(unlist(flens))## [1] 175988This is a count of gene-snp pairs satisfying structural and genetic criteria.
Visualizing and annotating significant loci
Re-binding probe annotation from RangedSummarizedExperiment
In the case of geuFPKM there is some relevant metadata in the rowRanges element. We will bind that into the collection of significant findings.
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]## GRanges object with 3 ranges and 19 metadata columns:
## seqnames ranges strand | paramRangeID REF
## <Rle> <IRanges> <Rle> | <factor> <DNAStringSet>
## 1 1 [ 940005, 940005] * | ENSG00000215915.5 A
## 1 1 [ 941539, 941539] * | ENSG00000215915.5 C
## 1 1 [1357992, 1357992] * | ENSG00000215915.5 C
## ALT chisq permScore_1 permScore_2 permScore_3
## <CharacterList> <numeric> <numeric> <numeric> <numeric>
## 1 G 6.019906 0.10806376 0.0003276232 0.08657792
## 1 T 6.410836 0.09361884 0.0072445068 0.07561893
## 1 T 6.103567 1.37837091 0.1842075347 0.67166832
## permScore_4 permScore_5 permScore_6 snp MAF
## <numeric> <numeric> <numeric> <character> <numeric>
## 1 0.004372199 0.07662322 0.35807525 rs2799056 0.41460674
## 1 0.393588577 0.07313992 0.34126371 rs9778087 0.41910112
## 1 3.379391274 1.35579611 0.01747091 rs3737716 0.06516854
## probeid mindist estFDR gene_id gene_status
## <character> <numeric> <numeric> <character> <character>
## 1 ENSG00000215915.5 445064 0.04880831 ENSG00000215915.5 KNOWN
## 1 ENSG00000215915.5 443530 0.03146240 ENSG00000215915.5 KNOWN
## 1 ENSG00000215915.5 27077 0.04455828 ENSG00000215915.5 KNOWN
## gene_type gene_name
## <character> <character>
## 1 protein_coding ATAD3C
## 1 protein_coding ATAD3C
## 1 protein_coding ATAD3C
## -------
## seqinfo: 86 sequences from hg19 genomeStatic visualization of FDR patterns
We have a utility to create an annotated Manhattan plot for a search cis to a gene. The basic ingredients are
- a
ciseStoreinstance for basic location and association information - a gene identifier that works for that store
- an
FDRsuppinstance that includes the function that maps from association scores to FDR, and the filter employed during FDR estimation - an annotation resource; here we use ChromHMM labeling based on NA12878, in the
hmm878GRanges instance in gQTLstats/data.
It is important to recognize that, given an FDRsupp instance we can compute the FDR for any association score, but validity of the FDR attribution requires that we refrain from computing it for any locus excluded by filtering. the manhWngr executes the FDRsupp-resident filter by default.
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 )## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called '<S4 object of class
## structure("OrganismDb", package = "OrganismDbi")>'## 'select()' returned 1:1 mapping between keys and columns