Public Interface¶
LMM Core Classes¶
-
class
limix_lmm.lmm_core.
LMMCore
(y, F, Ki_dot=None)[source]¶ Core LMM for interaction and association testing.
The model is
\[\mathbf{y}\sim\mathcal{N}( \underbrace{\mathbf{F}\mathbf{b}}_{\text{covariates}}+ \underbrace{(\mathbf{X}\hat{\odot}\mathbf{I})\; \boldsymbol{\beta}}_{\text{genetics}}, \sigma^2\underbrace{\mathbf{K}_{ \boldsymbol{\theta}_0}}_{\text{covariance}})\]\(\hat{\odot}\) is defined as follows. Let
\[\mathbf{A}=\left[\mathbf{a}_1, \dots, \mathbf{a}_{m} \right]\in\mathbb{R}^{n\times{m}}]\]and
\[\mathbf{B}=\left[\mathbf{b}_1, \dots, \mathbf{b}_{l} \right]\in\mathbb{R}^{n\times{l}}\]then
\[\mathbf{A}\hat{\odot}\mathbf{B} = \left[\mathbf{a}_1\odot\mathbf{b}_1 \dots, \mathbf{a}_1\odot\mathbf{b}_{l}, \mathbf{a}_2\odot\mathbf{b}_1, \dots, \mathbf{a}_2\odot\mathbf{b}_{l}, \dots, \mathbf{a}_{m}\odot\mathbf{b}_{l}, \dots,\mathbf{a}_{m}\odot\mathbf{b}_{l} \right]\in\mathbb{R}^{n\times{(ml)}}\]where \(\odot\) is the element-wise product.
Importantly, \(\mathbf{K}_{\boldsymbol{\theta}_0}\) is not re-estimated under the laternative model, i.e. only \(\sigma^2\) is learnt. The covariance parameters \(\boldsymbol{\theta}_0\) need to be learnt using another module, e.g. limix-core. Note that as a consequence of this, such an implementation of single-variant analysis does not require the specification of the whole covariance but just of a method that specifies the product of its inverse by a vector.
The test \(\boldsymbol{\beta}\neq{0}\) is done in bocks of
step
variants, wherestep
can be specifed by the user.Parameters: - y ((N, 1) ndarray) – phenotype vector
- F ((N, L) ndarray) – fixed effect design for covariates.
- Ki_dot (function) – method that takes an array and returns the dot product of the inverse of the covariance and the input array.
Examples
Example with Inter and step=1
>>> from numpy.random import RandomState >>> import scipy as sp >>> from limix_lmm.lmm_core import LMMCore >>> from limix_core.gp import GP2KronSumLR >>> from limix_core.covar import FreeFormCov >>> random = RandomState(1) >>> from numpy import set_printoptions >>> set_printoptions(4) >>> >>> N = 100 >>> k = 1 >>> m = 2 >>> S = 1000 >>> y = random.randn(N,1) >>> E = random.randn(N,k) >>> G = 1.*(random.rand(N,S)<0.2) >>> F = sp.concatenate([sp.ones((N,1)), random.randn(N,1)], 1) >>> Inter = random.randn(N, m) >>> >>> gp = GP2KronSumLR(Y=y, Cn=FreeFormCov(1), G=E, F=F, A=sp.ones((1,1))) >>> gp.covar.Cr.setCovariance(0.5*sp.ones((1,1))) >>> gp.covar.Cn.setCovariance(0.5*sp.ones((1,1))) >>> info_null = gp.optimize(verbose=False) >>> >>> lmm = LMMCore(y, F, gp.covar.solve) >>> lmm.process(G, Inter) >>> pv = lmm.getPv() >>> beta = lmm.getBetaSNP() >>> beta_ste = lmm.getBetaSNPste() >>> lrt = lmm.getLRT() >>> >>> print(pv.shape) (1000,) >>> print(beta.shape) (2, 1000) >>> print(pv[:4]) [0.1428 0.3932 0.4749 0.3121] >>> print(beta[:,:4]) [[ 0.0535 0.2571 0.122 -0.2328] [ 0.3418 0.3179 -0.2099 0.0986]]
Example with step=4
>>> lmm.process(G, step=4) >>> pv = lmm.getPv() >>> beta = lmm.getBetaSNP() >>> lrt = lmm.getLRT() >>> >>> print(pv.shape) (250,) >>> print(beta.shape) (4, 250) >>> print(pv[:4]) [0.636 0.3735 0.7569 0.3282] >>> print(beta[:,:4]) [[-0.0176 -0.4057 0.0925 -0.4175] [ 0.2821 0.2232 -0.2124 0.2433] [-0.0575 0.1528 -0.0384 0.019 ] [-0.2156 -0.2327 -0.1773 -0.1497]]
Example with step=4 and Inter
>>> lmm = LMMCore(y, F, gp.covar.solve) >>> lmm.process(G, Inter=Inter, step=4) >>> pv = lmm.getPv() >>> beta = lmm.getBetaSNP() >>> lrt = lmm.getLRT() >>> >>> print(pv.shape) (250,) >>> print(beta.shape) (8, 250) >>> print(pv[:4]) [0.1591 0.2488 0.4109 0.5877] >>> print(beta[:,:4]) [[ 0.0691 0.1977 0.435 -0.3214] [ 0.1701 -0.3195 0.0179 0.3344] [ 0.1887 -0.0765 -0.3847 0.1843] [-0.2974 0.2787 0.2427 -0.0717] [ 0.3784 0.0854 0.1566 0.0652] [ 0.4012 0.5255 -0.1572 0.1674] [-0.413 0.0278 0.1946 -0.1199] [ 0.0268 -0.0317 -0.1059 0.1414]]
Plot functions¶
-
limix_lmm.plot.
plot_manhattan
(ax, df, pv_thr=None, colors=None, offset=None, callback=None)[source]¶ Utility function to make manhattan plot
Parameters: - ax (pyplot plot) – subplot
- df (pandas.DataFrame) – pandas DataFrame with chrom, pos and pv
- colors (list) – colors to use in the manhattan plot
- offset (float) – offset between in chromosome expressed as fraction of the length of the longest chromosome (default is 0.2)
- callback (function) – callback function that takes as input df
Examples
>>> from matplotlib import pyplot as plt >>> from limix_lmm.plot import plot_manhattan >>> import scipy as sp >>> import pandas as pd >>> n_chroms = 5 >>> n_snps_per_chrom = 10000 >>> chrom = sp.kron(sp.arange(1, n_chroms + 1), sp.ones(n_snps_per_chrom)) >>> pos = sp.kron(sp.ones(n_chroms), sp.arange(n_snps_per_chrom)) >>> pv = sp.rand(n_chroms * n_snps_per_chrom) >>> df = pd.DataFrame({'chrom': chrom, 'pos': pos, 'pv': pv}) >>> >>> ax = plt.subplot(111) >>> plot_manhattan(ax, df)
-
limix_lmm.plot.
qqplot
(ax, pv, color='k', plot_method=None, line=False, pv_thr=0.01, plot_xyline=False, xy_labels=False, xyline_color='r')[source]¶ Utility function to make manhattan plot
Parameters: - ax (pyplot plot) – subplot
- df (pandas.DataFrame) – pandas DataFrame with chrom, pos and pv
- color (color) – colors to use in the manhattan plot (default is black)
- plot_method (function) – function that takes x and y and plots a custom scatter plot The default value is None.
- pv_thr (float) – threshold on pvs
- plot_xyline (bool) – if True, the x=y line is plotteed. The default value is False.
Examples
>>> from limix_lmm.plot import qqplot >>> import scipy as sp >>> from matplotlib import pyplot as plt >>> >>> pv1 = sp.rand(10000) >>> pv2 = sp.rand(10000) >>> pv3 = sp.rand(10000) >>> >>> ax = plt.subplot(111) >>> qqplot(ax, pv1, color='C0') >>> qqplot(ax, pv2, color='C1') >>> qqplot(ax, pv3, color='C2', plot_xyline=True, xy_labels=True)