Association Scan: Every Column of x Against Every Column of y, With a GPU

Tests every column of x against every column of y by linear regression, lm(y[, j] ~ x[, i] + covariates), and returns the significant associations. It is the scan of expression quantitative trait loci (eQTL) studies, genotypes against gene expression, and of other many-by-many association screens; it computes what Matrix eQTL's linear model does (Shabalin 2012), on the GPU.

Usage

gpuAssoc(x, y, covariates = NULL, pvalue = 1e-5, top = NULL,
    limit = 1e7, tile = 2048L,
    precision = getOption("gputools.precision", "single"))

Arguments

x: a numeric matrix, a column per tested variable (a SNP's genotypes, say, coded 0, 1, 2), a row per sample.
y: a numeric matrix, a column per outcome (a gene's expression), the same rows.
covariates: NULL or a numeric matrix of covariates to adjust for (sex, age, principal components of ancestry, ...), the same rows. An intercept is always included.
pvalue: report every association with a p-value at most this.
top: or instead, report each column of y's top most significant columns of x; pvalue is then not used.
limit: the most associations to return: a pvalue that lets more through is an error, not a result of unexpected size.
tile: the block size the GPU works in; see gpuCorPairs. It has no effect on the result.
precision: "single" or "double": the floating point precision the GPU computes in; see gpuCor.

Details

For each pair, the slope of x[, i] in the least-squares fit of y[, j] on it, an intercept and the covariates is tested for zero with the t test summary(lm(...)) reports, on \(n - 2 - k\) degrees of freedom for \(n\) samples and \(k\) covariates. As in Matrix eQTL, the intercept and covariates are projected out of every column of x and y once, after which each test is a correlation: the scan is one large matrix product, computed a block at a time and filtered on the GPU.

A column of x or y that the covariates explain completely, or a constant one, has no test and is never reported. Missing values are not allowed: impute them first (Matrix eQTL uses the SNP's mean genotype) or drop the sample.

The whole of x and y is copied to the GPU. For more SNPs than fit, call gpuAssoc on a chromosome or a block of SNPs at a time: the results are the same apart from the false discovery rates, which then count only that call's tests.

Value

A data frame, a row per reported association:

i: the column of x.
j: the column of y.
x, y: the columns' names, if they have them.
beta: the slope of x[, i].
t: its t statistic.
p.value: its two-sided p-value.
fdr: the Benjamini-Hochberg false discovery rate over all ncol(x) * ncol(y) tests. It is computed from the reported p-values alone, as Matrix eQTL computes it, and so can overstate the rate but never understate it.

With pvalue, rows are ordered from the most significant; with top, by j and then from the most significant.

References

Shabalin, A. A. (2012) Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28, 1353–1358.

Examples

if (FALSE) { # gpuAvailable()
set.seed(1)
samples <- 200
snps <- matrix(rbinom(samples * 1000, 2, 0.3), samples, 1000,
               dimnames = list(NULL, paste0("rs", 1:1000)))
age <- rnorm(samples, 50, 10)
genes <- matrix(rnorm(samples * 300), samples, 300,
                dimnames = list(NULL, paste0("gene", 1:300)))
genes[, 7] <- genes[, 7] + 0.6 * snps[, 42] + 0.05 * age   # one real eQTL

gpuAssoc(snps, genes, covariates = age, pvalue = 1e-6)

# the same test by lm()
summary(lm(genes[, 7] ~ snps[, 42] + age))$coefficients[2, ]
}