Merge pull request #73 from zhanxw/fan-bare

Fan bare
zhanxw · Dec 26, 2018 · ccc16bb · ccc16bb
2 parents ddaf2c8 + 2359352
commit ccc16bb
Show file tree

Hide file tree

Showing 4 changed files with 332 additions and 12 deletions.
diff --git a/misc/combineKinship.py b/misc/combineKinship.py
@@ -26,28 +26,46 @@ def getVariant(fn):
     return int(ret)
 
 def loadKinship(fn):
+    import numpy as np
     ids = []
-    kin = []
+    kin = None
     ncol = -1
     for i, ln in enumerate(myopen(fn)):
         fd = ln.strip().split()
         if i == 0:
             ncol = len(fd)
+            nsample = ncol - 2
+            kin = np.empty([nsample, nsample])
             continue
         ids.append(fd[:2])
-        kin.append([float(i) for i in fd[2:]])
+        kin[i - 1, ] = [float(ii) for ii in fd[2:]]
+        if i % 1000 == 0:
+            print 'i = ', i
+
     nsample = ncol - 2
     if len(ids) != nsample or \
        any([nsample != len(i) for i in kin]):
         print >> sys.stderr, "Dimension not match in ", fn
+
+    if kin.shape[0] != kin.shape[1] or kin.shape[0] != nsample:
+        print >> sys.stderr, "Dimension not match in ", fn
+
+    # save to disk
+    fnDisk = 'tmp.' + os.path.basename(fn) + '.mmap.npy'
+    kinDisk = np.memmap(fnDisk, dtype='float32', mode='w+', shape=kin.shape)
+    kinDisk[:] = kin[:]
+    del kinDisk
+    kin = np.memmap(fnDisk, dtype='float32', mode='r', shape=kin.shape)
+
     print >> sys.stderr, "Kinship file %s with %d samples loaded" % (fn, nsample)
     return ids, kin
 
 def combineKinship(kin1, cnt1, kin2, cnt2):
-    kin = kin1
-    for i in xrange(len(kin1)):
-        for j in xrange(len(kin1[0])):
-            kin[i][j] = (kin1[i][j] * cnt1 + kin2[i][j] * cnt2) / (cnt1 + cnt2)
+    # kin = kin1
+    # for i in xrange(len(kin1)):
+    #     for j in xrange(len(kin1[0])):
+    #         kin[i][j] = (kin1[i][j] * cnt1 + kin2[i][j] * cnt2) / (cnt1 + cnt2)
+    kin = (kin1 * cnt1 + kin2 * cnt2 ) / (cnt1 + cnt2)
     cnt = cnt1 + cnt2
     print >> sys.stderr, "Combined %d variants from %d and %d variants" % (cnt, cnt1, cnt2)
     return kin, cnt
@@ -99,11 +117,13 @@ def isSameId(id1, id2):
     print >> sys.stderr, "Total %d variants from %d files are recognized." % (sum(numAuto), len(fn))
 
     # read first kinship
+    print >> sys.stderr, "Process kinship file [ %s ]" % (fn[0])    
     ids, kin = loadKinship(fn[0])
     cnt = numAuto[0]
 
     # gradually combine rest kinships
     for i in xrange(1, len(fn)):
+        print >> sys.stderr, "Process kinship file [ %s ]" % (fn[i])
         ids2, kin2 = loadKinship(fn[i])
         cnt2 = numAuto[i]
         if not isSameId(ids, ids2):

diff --git a/src/Model.cpp b/src/Model.cpp
@@ -903,6 +903,7 @@ int MetaCovTest::fitWithGivenGenotype(const Matrix& genotype,
     // copyCovariateAndIntercept(genotype.rows, covariate, &cov);
     fitOK = (0 == model->FitNullModel(genotype, dc));
     if (!fitOK) return -1;
+    // always prepare covZZ, covZZInv when null model is fitted
     model->calculateZZ(&this->covZZ);
     CholeskyInverseMatrix(this->covZZ, &this->covZZInv);
     model->needToFitNullModel = false;
@@ -925,11 +926,6 @@ int MetaCovTest::fitWithGivenGenotype(const Matrix& genotype,
     if (nCovariate) {
       model->calculateXZ(x, xz);
     }
-    // if (model->needToFitNullModel || dc->isPhenotypeUpdated() ||
-    //     dc->isCovariateUpdated()) {
-    //   model->calculateZZ(&this->covZZ);
-    //   CholeskyInverseMatrix(this->covZZ, &this->covZZInv);
-    // }
   }
   fitOK = true;
   return 0;

diff --git a/vcfUtils/Makefile b/vcfUtils/Makefile
@@ -16,7 +16,8 @@ UTIL_EXEC = vcf2plink vcfSummary vcfConcordance \
             vcf2kinship \
             vcfPeek \
             vcf2ld_neighbor \
-            kinshipDecompose 
+            kinshipDecompose \
+            combineKinship
             # vcf2merlin 
 
 DIR_EXEC = ../executable