In [27]:
from pkg_resources import resource_filename
from limmbo.io.reader import ReadData
from limmbo.io.utils import file_type

import numpy as np

# data container
data = ReadData(verbose=False)

# Read covariates
file_covs = "./data/limmbo/limmbo_covariates.csv"
data.getCovariates(file_covariates=file_covs)

# Read genotypes in delim-format
file_geno =  './data/plink_files/atchely_imputed'
data.getGenotypes(file_genotypes=file_geno)
data.genotypes

# Read phenotypes
file_pheno = './data/limmbo/limmbo_phenotypes.csv'
data.getPhenotypes(file_pheno=file_pheno)
data.phenotypes

# Read relatedness
file_relatedness = './data/limmbo/limmbo_relatedness.csv'
data.getRelatedness(file_relatedness=file_relatedness)

In [20]:
from limmbo.io import input

indata = input.InputData(verbose=False)
indata.addGenotypes(genotypes=data.genotypes,
                    genotypes_info=data.genotypes_info)

indata = input.InputData(verbose=False)
indata.addGenotypes(genotypes=data.genotypes,
                    genotypes_info=data.genotypes_info,
                    geno_samples=data.geno_samples)
indata.addPhenotypes(phenotypes = data.phenotypes, 
                     pheno_samples = data.pheno_samples,
                     phenotype_ID = data.phenotype_ID)
indata.addRelatedness(relatedness = data.relatedness)
indata.addCovariates(covariates = data.covariates,
                     covs_samples = data.covs_samples)

In [21]:
from limmbo.core.vdsimple import vd_reml

Cg, Cn, ptime = vd_reml(indata, verbose=False)
indata.addVarianceComponents(Cg = Cg, Cn = Cn)
Cg


Out[21]:
array([[162.58949648,  23.87759381],
       [ 23.87759381,   4.86702787]])

In [22]:
np.var(data.phenotypes['Final_weight'])


Out[22]:
23.56776102459051

In [24]:
indata.regress()
indata.transform(transform="scale")

In [25]:
from limmbo.core.gwas import GWAS

gwas = GWAS(datainput=indata, verbose=False)
gwas.name = "test"

In [26]:
resultsAssociation = gwas.runAssociationAnalysis(setup="lmm", mode="multitrait")


---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-26-22596b011515> in <module>()
----> 1 resultsAssociation = gwas.runAssociationAnalysis(setup="lmm", mode="multitrait")

/home/diogro/bin/anaconda3/envs/limmbo/lib/python2.7/site-packages/limmbo/core/gwas.pyc in runAssociationAnalysis(self, mode, setup, adjustSingleTrait)
    201         if mode == "multitrait":
    202             associationResults = self.__multiTraitAssociation_anyeffect(
--> 203                 genotypes=self.genotypes)
    204 
    205         if mode == "singletrait":

/home/diogro/bin/anaconda3/envs/limmbo/lib/python2.7/site-packages/limmbo/core/gwas.pyc in __multiTraitAssociation_anyeffect(self, genotypes, empiricalP, computeFDR)
    279                                                  covs=self.covariates, K1r=K1r,
    280                                                  K1c=K1c, K2c=K2c,
--> 281                                                  searchDelta=self.searchDelta)
    282 
    283         if not empiricalP and not computeFDR:

/home/diogro/bin/anaconda3/envs/limmbo/lib/python2.7/site-packages/limix/qtl/qtl.pyc in qtl_test_lmm_kronecker(snps, phenos, covs, Acovs, Asnps, K1r, K1c, K2r, K2c, NumIntervalsDelta0, NumIntervalsDeltaAlt, searchDelta)
    365         # add SNP design
    366         lmm.setSNPcoldesign(Asnps[iA])
--> 367         lmm.process()
    368         pv[iA, :] = lmm.getPv()[0]
    369     return lmm, pv

/home/diogro/bin/anaconda3/envs/limmbo/lib/python2.7/site-packages/limix_legacy/deprecated/core.pyc in process(self)
  12046 
  12047         """
> 12048         return _core.CKroneckerLMM_process(self)
  12049 
  12050 

RuntimeError: LIMIX error: dimensions in background model inconsistent

In [118]:
gwas.computeFDR(0.05)


Out[118]:
(0.9926760802230591,
 array([9.78043793e-33, 9.78043793e-33, 9.78043793e-33, ...,
        1.00000000e+00, 1.00000000e+00, 1.00000000e+00]))