Merge branch 'develop'

limix/io/_csv.py

@@ -5,10 +5,12 @@ def read_csv(filename, sep=None, header=True):
     ----------
     filename : str
         Path to a CSV file.
+    sep : str
+        Separator.
 
     Returns
     -------
-    data : dask dataframe
+    data : dask dataframes
 
     Examples
     --------

limix/io/plink.py

@@ -88,3 +88,18 @@ def see_kinship(filepath):
         K = read_grm_raw(filepath)
 
     limix.plot.plot_kinship(K)
+
+
+def _print_title(title, msg):
+    k = msg.find("\n") - len(title) - 2
+    left = ("-" * (k // 2)) + " "
+    right = " " + ("-" * (k // 2 + k % 2))
+    print(left + title + right)
+    print(msg)
+
+
+def see_bed(filepath):
+    (bim, fam, bed) = read_plink(filepath)
+
+    _print_title("Samples", repr(bim))
+    _print_title("Genotype", repr(fam))

limix/io/util.py

@@ -1,3 +1,5 @@
+from os.path import exists
+
 def file_type(filepath):
     imexts = ['.png', '.bmp', '.jpg', 'jpeg']
     if filepath.endswith('.hdf5') or filepath.endswith('.h5'):
@@ -8,6 +10,11 @@ def file_type(filepath):
         return 'grm.raw'
     if filepath.endswith('.npy'):
         return 'npy'
+    if _is_bed(filepath):
+        return 'bed'
     if any([filepath.endswith(ext) for ext in imexts]):
         return 'image'
     return 'unknown'
+
+def _is_bed(filepath):
+    return all([exists(filepath + ext) for ext in ['.bed', '.bim', '.fam']])

limix/plot/qqplot.py

@@ -36,7 +36,7 @@ def _qqplot_bar(M=1000000, alphaLevel=0.05, distr='log10'):
     return betaUp, betaDown, theoreticalPvals
 
 
-def qqplot(pv, distr='log10', alphaLevel=0.05, ax=None):
+def qqplot(pv, label='unknown', distr='log10', alphaLevel=0.05, ax=None):
     r"""Produces a Quantile-Quantile plot of the observed P value
         distribution against the theoretical one under the null.
 
@@ -104,7 +104,7 @@ def qqplot(pv, distr='log10', alphaLevel=0.05, ax=None):
 
     plt.plot(qnull, qemp, '.')
     # plt.plot([0,qemp.m0x()], [0,qemp.max()],'r')
-    plt.plot([0, qnull.max()], [0, qnull.max()], 'r')
+    plt.plot([0, qnull.max()], [0, qnull.max()], 'r', label=label)
     plt.ylabel(xl)
     plt.xlabel(yl)
     if alphaLevel is not None:

limix/scripts/limix.py

@@ -17,6 +17,8 @@ def do_see(args):
         limix.io.plink.see_kinship(args.file)
     elif ft == 'npy':
         limix.io.npy.see_kinship(args.file)
+    elif ft == 'bed':
+        limix.io.plink.see_bed(args.file)
     elif ft == 'image':
         limix.plot.see_image(args.file)
     else:

limix/stats/kinship.py

@@ -1,11 +1,43 @@
-from numpy import copyto, sqrt
+from numpy import ascontiguousarray, copyto, sqrt, zeros
+from tqdm import tqdm
 
 
-def linear_kinship(G):
-    G = G - G.mean(0)
-    G /= G.std(0)
-    G /= sqrt(G.shape[1])
-    return G.dot(G.T)
+def linear_kinship(G, out=None, progress=True):
+    r"""Estimate Kinship matrix via linear kernel.
+
+    Examples
+    --------
+    .. doctest::
+
+        >>> from numpy.random import RandomState
+        >>> from limix.stats import linear_kinship
+        >>>
+        >>> random = RandomState(1)
+        >>> X = random.randn(4, 100)
+        >>> K = linear_kinship(X, progress=False)
+        >>> print(K)
+        [[ 27.3944  -5.785  -10.2402 -11.3693]
+         [ -5.785   26.9655  -7.068  -14.1125]
+         [-10.2402  -7.068   28.7332 -11.425 ]
+         [-11.3693 -14.1125 -11.425   36.9068]]
+    """
+    (n, p) = G.shape
+    if out is None:
+        out = zeros((n, n))
+
+    nsteps = min(30, p)
+
+    for i in tqdm(range(nsteps), disable=not progress):
+        start = i * (p // nsteps)
+        stop = min(start + p // nsteps, p)
+
+        G = G - G.mean(0)
+        G /= G.std(0)
+        G /= sqrt(p)
+
+        out += ascontiguousarray(G.dot(G.T), float)
+
+    return out
 
 
 def gower_norm(K, out=None):
@@ -16,22 +48,22 @@ def gower_norm(K, out=None):
     Examples
     --------
 
-        .. doctest::
-
-            >>> from numpy.random import RandomState
-            >>> from limix.stats import gower_norm
-            >>> import scipy as sp
-            >>>
-            >>> random = RandomState(1)
-            >>> X = random.randn(4, 4)
-            >>> K = sp.dot(X,X.T)
-            >>> Z = random.multivariate_normal(sp.zeros(4), K, 50)
-            >>> print("%.3f" % sp.mean(Z.var(1,ddof=1)))
-            2.335
-            >>> Kn = gower_norm(K)
-            >>> Zn = random.multivariate_normal(sp.zeros(4), Kn, 50)
-            >>> print("%.3f" % sp.mean(Zn.var(1, ddof=1)))
-            0.972
+    .. doctest::
+
+        >>> from numpy.random import RandomState
+        >>> from limix.stats import gower_norm
+        >>> import scipy as sp
+        >>>
+        >>> random = RandomState(1)
+        >>> X = random.randn(4, 4)
+        >>> K = sp.dot(X,X.T)
+        >>> Z = random.multivariate_normal(sp.zeros(4), K, 50)
+        >>> print("%.3f" % sp.mean(Z.var(1,ddof=1)))
+        2.335
+        >>> Kn = gower_norm(K)
+        >>> Zn = random.multivariate_normal(sp.zeros(4), Kn, 50)
+        >>> print("%.3f" % sp.mean(Zn.var(1, ddof=1)))
+        0.972
     """
 
     c = (K.shape[0] - 1) / (K.trace() - K.mean(0).sum())

setup.py

@@ -48,7 +48,7 @@ def setup_package():
 
     metadata = dict(
         name='limix',
-        version='1.0.7',
+        version='1.0.8',
         maintainer="Limix Developers",
         maintainer_email="horta@ebi.ac.uk",
         author=("Christoph Lippert, Danilo Horta, " +