Single-variant testing

Linear models

limix.qtl.qtl_test_lm(snps, pheno, covs=None, test=’lrt’, verbose=None)[source]

Wrapper function for univariate single-variant association testing using a linear model.

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • covs (ndarray, optional) – (N, D) ndarray of D covariates for N individuals. By default, covs is a (N, 1) array of ones.
  • test ({'lrt', 'f'}, optional) – test statistic. ‘lrt’ for likelihood ratio test (default) or ‘f’ for F-test.
  • verbose (bool, optional) – if True, details such as runtime as displayed.
Returns:

LIMIX LMM object

Return type:

limix.qtl.LMM

Examples

>>> from numpy.random import RandomState
>>> from limix.qtl import qtl_test_lm
>>> from numpy import set_printoptions
>>> set_printoptions(4)
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>>
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, 1)
>>>
>>> lm = qtl_test_lm(snps, pheno)
>>> print(lm.getPv()[:,:4])
[[ 0.8796  0.5065  0.5666  0.6016]]

Linear mixed models

limix.qtl.qtl_test_lmm(snps, pheno, K=None, covs=None, test=’lrt’, NumIntervalsDelta0=100, NumIntervalsDeltaAlt=100, searchDelta=False, verbose=None)[source]

Wrapper function for univariate single-variant association testing using a linear mixed model.

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • K (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients. If not provided, then standard linear regression is considered.
  • covs (ndarray, optional) – (N, D) ndarray of D covariates for N individuals. By default, covs is a (N, 1) array of ones.
  • test ({'lrt', 'f'}, optional) – test statistic. ‘lrt’ for likelihood ratio test (default) or ‘f’ for F-test.
  • NumIntervalsDelta0 (int, optional) – number of steps for delta optimization on the null model. By default NumIntervalsDelta0 is 100.
  • NumIntervalsDeltaAlt (int, optional) – number of steps for delta optimization on the alternative model. Requires searchDelta=True to have an effect.
  • searchDelta (bool, optional) – if True, delta optimization on the alternative model is carried out. By default searchDelta is False.
  • verbose (bool, optional) – if True, details such as runtime as displayed.
Returns:

LIMIX LMM object

Return type:

limix.qtl.LMM

Examples

>>> from numpy.random import RandomState
>>> from numpy import dot
>>> from limix.qtl import qtl_test_lmm
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>>
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, 1)
>>> W = random.randn(N, 10)
>>> kinship = dot(W, W.T) / float(10)
>>>
>>> lmm = qtl_test_lmm(snps, pheno, kinship)
>>> print(lmm.getPv()[:,:4])
[[ 0.8571  0.4668  0.5872  0.5589]]
limix.qtl.qtl_test_interaction_lmm(snps, pheno, Inter, Inter0=None, covs=None, K=None, test=’lrt’)[source]

Wrapper function for single-variant interaction test using LMMs.

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • Inter (ndarray, optional) – (N, I) ndarray of I interaction variables to be tested for N individuals.
  • Inter0 (ndarray, optional) – (N, I0) ndarray of I0 interaction variables to be included for N individuals to be included both in the alternative and in the null model. By default Inter0 is a (N, 1) array of ones.
  • covs (ndarray, optional) – (N, D) ndarray of D covariates for N individuals. By default, covs is a (N, 1) array of ones.
  • K (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients (optional) If not provided, then standard linear regression is considered.
  • test ({'lrt', 'f'}, optional) – test statistic. ‘lrt’ for likelihood ratio test (default) or ‘f’ for F-test.
Returns:

limix CInteractLMM object

Return type:

limix_legacy.deprecated.CInteractLMM()

Examples

>>> from numpy.random import RandomState
>>> from numpy import dot, eye, ones, concatenate
>>> from limix.qtl import qtl_test_interaction_lmm
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>>
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, 1)
>>> W = random.randn(N, 10)
>>> kinship = dot(W, W.T) / float(10)
>>> kinship+= 1e-4 * eye(N)
>>>
>>> #environments to include in null and alt models
>>> Inter0 = ones((N,1))
>>>
>>> #environments to include in alt model
>>> Inter = random.randn(N,2)
>>>
>>> #environments are also included as fixed effect covariates
>>> mean = ones([N, 1])
>>> covs = concatenate([mean, Inter], 1)
>>>
>>> #interaction test
>>> lmi = qtl_test_interaction_lmm(snps, pheno, Inter=Inter,
...                                Inter0=Inter0)
>>> pvi = lmi.getPv()
>>>
>>> print(pvi.shape)
(1, 1000)
>>> print(pvi[:,:4])
[[ 0.8179  0.2185  0.3338  0.077 ]]
limix.qtl.qtl_test_lmm_kronecker(snps, phenos, covs=None, Acovs=None, Asnps=None, K1r=None, K1c=None, K2r=None, K2c=None, NumIntervalsDelta0=100, NumIntervalsDeltaAlt=100, searchDelta=False)[source]

Wrapper function for single-variant multi-trait association testing class using linear mixed models.

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • covs (ndarray or list of ndarrays, optional) – sample design matrices of the different fixed effect terms. In the case of a single fixed effect term, covs can be specified as a (N, D) ndarray of D covariates for N individuals. In this case Acovs should also be a (K, P) ndarray, where K is the rank of the trait design and P is the number of traits. For multiple terms, covs is a list of ndarray, where each element covs[i] has shape (N, Di). In this case, Acovs should also be list with the same number of ndarray elements. Each element Acovs[i] has shape (Ki, P), where P is the number of analyzed phenotypes. By default covs is (N, 1) ndarray of ones and Acovs[i] is a (P, P) identity ndarray.
  • Acovs (ndarray or list of ndarrays, optional) – trait design matrices of the different fixed effect terms. See covs for more info.
  • Asnps (ndarray, optional) – (K, P) trait design of rank K of the snp fixed effect on the P traits. By default Asnps is (1, P) ndrray of ones.
  • K1r (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients. If not provided, K1r = numpy.dot(snps, snps.T) is assumed.
  • K1c (ndarray, optional) – (P, P) ndarray of LMM-covariance/kinship coefficients. If either K1c or K2c is not provided, K1c is estimated by (restricted) maximum likelihood from the model without the snp effect.
  • K2r (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients. By default K2r is a (N, N) identity ndarray.
  • K2c (ndarray, optional) – (P, P) ndarray of residual covariances. If either K1c or K2c is not provided, K1c is estimated by (restricted) maximum likelihood from the model without the snp effect.
  • NumIntervalsDelta0 (int, optional) – number of steps for delta optimization on the null model. By default NumIntervalsDelta0 is 100.
  • NumIntervalsDeltaAlt (int, optional) – number of steps for delta optimization on the alternative model. Requires searchDelta=True to have an effect.
  • searchDelta (bool, optional) – if True, delta optimization on the alternative model is carried out. By default searchDelta is False.
Returns:

tuple containing:
  • lmm (:class:`limix_legacy.deprecated.CKroneckerLMM()`): CKroneckerLMM limix_legacy object
  • pv (ndarray): P-values for all SNPs from liklelihood ratio test

Return type:

(tuple)

Examples

Example of multi-trait single-variant association testing using a linear mixed model and specifying a common effects.

>>> from numpy.random import RandomState
>>> from numpy import dot, eye, ones
>>> from limix.qtl import qtl_test_lmm_kronecker
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>> P = 2
>>>
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, P)
>>> W = random.randn(N, 10)
>>> kinship = dot(W, W.T) / float(10)
>>> kinship+= 1e-4 * eye(N)
>>>
>>> #common effect test
>>> Asnps = ones((1, P))
>>> lmm, pv = qtl_test_lmm_kronecker(snps, pheno, K1r=kinship,
...                                 Asnps=Asnps)
>>> beta = lmm.getBetaSNP()
>>>
>>> print(pv.shape)
(1, 1000)
>>> print(pv[:,:4])
[[ 0.8722  0.6509  0.3593  0.6816]]
>>>
>>> print(beta.shape)
(1, 1000)
>>> print(beta[:,:4])
[[ 0.0281  0.0758  0.1458  0.0701]]

Example showing how to set an any effect test. For more information on effect designs see limix_tutorials.

>>> #any effect test
>>> Asnps = eye(P)
>>> lmm, pv = qtl_test_lmm_kronecker(snps, pheno, K1r=kinship,
...                                 Asnps=Asnps)
>>> beta = lmm.getBetaSNP()
>>>
>>> print(pv.shape)
(1, 1000)
>>> print(pv[:,:4])
[[ 0.9838  0.8664  0.3555  0.5244]]
>>>
>>> print(beta.shape)
(2, 1000)
>>> print(beta[:,:4])
[[ 0.0413  0.0312 -0.0177 -0.098 ]
 [ 0.0139  0.1234  0.3204  0.2495]]
limix.qtl.qtl_test_interaction_lmm_kronecker(snps, phenos, covs=None, Acovs=None, Asnps1=None, Asnps0=None, K1r=None, K1c=None, K2r=None, K2c=None, NumIntervalsDelta0=100, NumIntervalsDeltaAlt=100, searchDelta=False, return_lmm=False)[source]

Wrapper function for single-variant multi-trait interaction test

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • covs (ndarray or list of ndarrays, optional) – sample design matrices of the different fixed effect terms. See limix.qtl.qtl_test_lmm_kronecker for more info.
  • Acovs (ndarray or list of ndarrays, optional) – trait design matrices of the different fixed effect terms. See limix.qtl.qtl_test_lmm_kronecker for more info.
  • Asnps1 (ndarray or list of ndarray, optional) – (K, P) trait design of rank K of the snp fixed effect on the P traits for the alternative model. Asnps1 is tested against the null trait design Asnps0. Multiple designs can be tested vs the null trait design If a list of ndarray trait designs is specified for Asnps1, the multiple designs are tested vs the null trait design Asnps0 one-by-one. By default Asnps1 is a (P, P) identity ndarray.
  • Asnps0 (ndarray, optional) – (K0, P) trait design of rank K of the snp fixed effect on the P traits for the null model. By default Asnps0 is (1, P) ndrray of ones.
  • K1r (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients. If not provided, K1r = numpy.dot(snps, snps.T) is assumed.
  • K1c (ndarray, optional) – (P, P) ndarray of LMM-covariance/kinship coefficients. If either K1c or K2c is not provided, K1c is estimated by (restricted) maximum likelihood from the model without the snp effect.
  • K2r (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients. By default K2r is a (N, N) identity ndarray.
  • K2c (ndarray, optional) – (P, P) ndarray of residual covariances. If either K1c or K2c is not provided, K1c is estimated by (restricted) maximum likelihood from the model without the snp effect.
  • NumIntervalsDelta0 (int, optional) – number of steps for delta optimization on the null model. By default NumIntervalsDelta0 is 100.
  • NumIntervalsDeltaAlt (int, optional) – number of steps for delta optimization on the alternative model. Requires searchDelta=True to have an effect.
  • searchDelta (bool, optional) – if True, delta optimization on the alternative model is carried out. By default searchDelta is False.
Returns:

tuple containing:
  • pv (ndarray): P-values of the interaction test (Asnps1 vs Asnps0) If Asnps1 is a list of trait designs, pv contains the P values for the different trait designs (as rows).
  • pv0 (ndarray): P-values of the association test Asnps0 vs 0.
  • pvAlt (ndarray): P-values of the association test Asnps1 vs 0. If Asnps1 is a list of trait designs, pvAlt contains the P values for the different trait designs (as rows).

Return type:

(tuple)

Examples

Example of how to perform a specific, common and any effect test. For more information see limix_tutorials.

>>> from numpy.random import RandomState
>>> from numpy import dot, eye, ones
>>> from limix.qtl import qtl_test_interaction_lmm_kronecker
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>> P = 2
>>>
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, P)
>>> W = random.randn(N, 10)
>>> kinship = dot(W, W.T) / float(10)
>>> kinship+= 1e-4 * eye(N)
>>>
>>> #interaction effect test
>>> Asnps0 = ones((1,P))
>>> Asnps1 = eye(P)
>>> pv, pv0, pvAlt = qtl_test_interaction_lmm_kronecker(snps, pheno,
...                                                     K1r=kinship,
...                                                     Asnps0=Asnps0,
...                                                     Asnps1=Asnps1)
>>> #interaction P value
>>> print(pv.shape)
(1, 1000)
>>> print(pv[:,:4])
[[ 0.9347  0.7744  0.2678  0.2893]]
>>>
>>> #common effect P value
>>> print(pv0.shape)
(1, 1000)
>>> print(pv0[:,:4])
[[ 0.8722  0.6509  0.3593  0.6816]]
>>>
>>> #any effect P value
>>> print(pvAlt.shape)
(1, 1000)
>>> print(pvAlt[:,:4])
[[ 0.9838  0.8664  0.3555  0.5244]]
limix.qtl.forward_lmm(snps, pheno, K=None, covs=None, qvalues=False, threshold=5e-08, maxiter=2, test=’lrt’, verbose=None, **kw_args)[source]

Wrapper function for univariate single-variant test with forward selection

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • K (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients (optional) If not provided, then standard linear regression is considered.
  • covs (ndarray, optional) – (N, D) ndarray of D covariates for N individuals. By default, covs is a (N, 1) array of ones.
  • qvalues (bool, optional) – if True, threshold is set to Storey qvalues to control FDR. (default is False)
  • threshold (float, optional) – P-value thrashold for inclusion in forward selection (default 5e-8). If qvalues=True, the threshold is set to to Storey qvalues.
  • maxiter (int, optional) – maximum number of interaction scans. First scan is without inclusion, so maxiter-1 inclusions can be performed. (default 2)
  • test ({'lrt', 'f'}, optional) – test statistic. ‘lrt’ for likelihood ratio test (default) or ‘f’ for F-test.
  • verbose (bool, optional) – print verbose output. (default is False)
Returns:

tuple containing:
  • lmm (:class:`limix.qtl.LMM`): LIMIX LMM object
  • res (dict): {iadded, pvadded, pvall}. iadded is an ndarray of indices of SNPs included in order of inclusion, pvadded is an ndarray of Pvalues obtained by the included SNPs in iteration and pvall is a (Nadded, S) ndarray of Pvalues for all iterations

Return type:

(tuple)

Examples

>>> from numpy.random import RandomState
>>> from numpy import dot, eye, ones
>>> from limix.qtl import forward_lmm
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>>
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, 1)
>>> W = random.randn(N, 10)
>>> kinship = dot(W, W.T) / float(10)
>>> kinship+= 1e-4 * eye(N)
>>>
>>> #forward lmm
>>> lmm, res = forward_lmm(snps, pheno, K=kinship, threshold=1.)
>>>
>>> print(res['pvall'].shape)
(2, 1000)
>>> print(res['pvall'][:,:4])
[[ 0.8571  0.4668  0.5872  0.5589]
 [ 0.77    0.4226  0.6165  0.8727]]
class limix.qtl.LMM(snps, pheno, K=None, covs=None, test=’lrt’, NumIntervalsDelta0=100, NumIntervalsDeltaAlt=100, searchDelta=False, verbose=None)[source]

Basic class for univariate single-variant association testing with LMMs.

Parameters:
  • snps (ndarray) – (N, S) ndarray of S SNPs for N individuals.
  • pheno (ndarray) – (N, P) ndarray of P phenotype sfor N individuals. If phenotypes have missing values, then the subset of individuals used for each phenotype column will be subsetted.
  • K (ndarray, optional) – (N, N) ndarray of LMM-covariance/kinship coefficients (optional) If not provided, then standard linear regression is considered.
  • covs (ndarray, optional) – (N, D) ndarray of D covariates for N individuals. By default, covs is a (N, 1) array of ones.
  • test ({'lrt', 'f'}, optional) – test statistic. ‘lrt’ for likelihood ratio test (default) or ‘f’ for F-test.
  • NumIntervalsDelta0 (int, optional) – number of steps for delta optimization on the null model. By default NumIntervalsDelta0 is 100.
  • NumIntervalsDeltaAlt (int, optional) – number of steps for delta optimization on the alternative model. Requires searchDelta=True to have an effect.
  • searchDelta (bool, optional) – if True, delta optimization on the alternative model is carried out. By default searchDelta is False.
  • verbose (bool, optional) – if True, details such as runtime as displayed.

Examples

Example of single-variant association testing using a linear mixed model.

>>> from numpy.random import RandomState
>>> from limix.qtl import LMM
>>> from numpy import dot
>>> from numpy import set_printoptions
>>> set_printoptions(4)
>>> random = RandomState(1)
>>>
>>> N = 100
>>> S = 1000
>>>
>>> # generate data
>>> snps = (random.rand(N, S) < 0.2).astype(float)
>>> pheno = random.randn(N, 1)
>>> W = random.randn(N, 10)
>>> kinship = dot(W, W.T) / float(10)
>>>
>>> # run single-variant associaiton testing with LMM
>>> lmm = LMM(snps, pheno, kinship)
>>> pv = lmm.getPv()
>>> beta = lmm.getBetaSNP()
>>> beta_ste = lmm.getBetaSNPste()
>>>
>>> print(pv.shape)
(1, 1000)
>>> print(pv[:,:4])
[[ 0.8571  0.4668  0.5872  0.5589]]
>>> print(beta[:,:4])
[[ 0.0436 -0.1695 -0.12    0.1388]]
>>> print(beta_ste[:,:4])
[[ 0.2422  0.2329  0.221   0.2375]]

As shown in the exmaple below, the method can be applied directly on multiple phenotypes. Each phenotype is indipendently tested against all variants one-by-one.

>>> random = RandomState(1)
>>>
>>> P = 3
>>> phenos = random.randn(N, P)
>>>
>>> lmm = LMM(snps, phenos, kinship)
>>> pv = lmm.getPv()
>>>
>>> print(pv.shape)
(3, 1000)
>>> print(pv[:,:4])
[[ 0.4712  0.694   0.3369  0.5221]
 [ 0.7336  0.8799  0.8109  0.1071]
 [ 0.0662  0.9203  0.2873  0.8268]]
getBetaSNP()[source]
Returns:(P, S) ndarray of SNP effect sizes.
Return type:ndarray
getBetaSNPste()[source]
Returns:(P, S) ndarray of standard errors over SNP effects.
Return type:ndarray
getPv()[source]
Returns:(P, S) ndarray of P-values.
Return type:ndarray

Generalised linear mixed models

limix.qtl.qtl_test_glmm(snps, pheno, lik, K, covs=None, test=’lrt’, NumIntervalsDeltaAlt=100, searchDelta=False, verbose=True)[source]

Wrapper function for univariate single-variant association testing using a generalised linear mixed model.

Parameters:
  • snps (array_like) – N individuals by S SNPs.
  • pheno (tuple, array_like) – Either a tuple of two arrays of N individuals each (Binomial phenotypes) or an array of N individuals (Poisson or Bernoulli phenotypes). It does not support missing values yet.
  • lik ({'bernoulli', 'binomial', 'poisson'}) – Sample likelihood describing the residual distribution.
  • K (array_like) – N by N covariance matrix (e.g., kinship coefficients).
  • covs (array_like, optional) – N individuals by D covariates. By default, covs is a (N, 1) array of ones.
  • test ({'lrt'}, optional) – Likelihood ratio test (default).
  • NumIntervalsDeltaAlt (int, optional) – number of steps for delta optimization on the alternative model. Requires searchDelta=True to have an effect.
  • searchDelta (bool, optional) – if True, delta optimization on the alternative model is carried out. By default searchDelta is False.
  • verbose (bool, optional) – if True, details such as runtime are displayed.
Returns:

LIMIX LMM object

Return type:

limix.qtl.LMM

Examples

>>> from numpy import dot, exp, sqrt
>>> from numpy.random import RandomState
>>> from limix.qtl import qtl_test_glmm
>>>
>>> random = RandomState(0)
>>>
>>> G = random.randn(250, 500) / sqrt(500)
>>> beta = 0.01 * random.randn(500)
>>>
>>> z = dot(G, beta) + 0.1 * random.randn(250)
>>> z += dot(G[:, 0], 1) # causal SNP
>>>
>>> y = random.poisson(exp(z))
>>>
>>> candidates = G[:, :5]
>>> K = dot(G[:, 5:], G[:, 5:].T)
>>> lm = qtl_test_glmm(candidates, y, 'poisson', K, verbose=False)
>>>
>>> print(lm.getPv())
[[ 0.0694  0.3336  0.5899  0.7388  0.7796]]