Commits

mimakaev committed 189d882

Bugfix in highResBinnedData

Comments (0)

Files changed (2)

src/hiclib/highResBinnedData.py

 
     def _checkConsistency(self):
         "Checks if shapes of all datasets are correct"
-
+        self._hasData()
         for i in self.data:
             shape = self.data[i].shape
             chr1, chr2 = i
                 raise ValueError("Wrong dimensions of a chromosomes {2}: expected {0},"
                                  "got {1}".format(expected, shape, i))
 
+    def _hasData(self):
+        mykeys = self.data.keys()
+        for i in xrange(self.genome.chrmCount):
+            if ((i, i)) not in mykeys:
+                print "Not all cis keys found in data!"
+                print "Please load data!"
+                raise ValueError("Please load data first")
+
+
     def _marginalError(self, marginals=None):
         "Checks after each pass of IC, if marginals are close enough to 1"
         if marginals is None:
 
         self._checkConsistency()
 
+
     def getMarginals(self, normalizeForIC=False):
         """
         Returns a sum over each row/column, and saves them to self.marginals
 
         normalizeForIc=True will normalize them to mean 1.
         """
-
+        self._hasData()
         marginals = [np.zeros(i, float) for i in self.genome.chrmLensBin]
         for chr1, chr2 in self.data:
             m2, m1 = self.data[(chr1, chr2)].getSums()
         If the vector is not provided, will divide it by a
         marginals calculated previously.
         """
+        self._hasData()
         if marginals is None:
             marginals = self.marginals
 
 
         Tries to set biases to self.biases after IC.
         """
-
+        self._hasData()
         curPass = 0
         marginals = np.zeros(self.genome.numBins, float)
         while self._marginalError() > tolerance:
         Removes diagonal from the data.
         m=0 is main diagonal, m=1 is main +-1, etc.
         """
+        self._hasData()
         for i in self.cisKeys:
             data = self.data[i].getData()
             removeDiagonals(data, m)
         """
         Removes "percent" percent of regions with low coverage.
         """
+        self._hasData()
         marg = self.getMarginals(normalizeForIC=False)
         allMarg = np.concatenate(marg)
         cut = np.percentile(allMarg[allMarg > 0], percent)
         mydict = h5dict(filename)
         for i in self.allKeys:
             data = self.data[i].getData()
-            mydict["heatmap%d %d"] = data
+            mydict["heatmap%d %d" % (i, i)] = data
         mydict["resolution"] = self.resolution
 
 

tests/fragmentHiC/testFragmentHiC.py

 
 if os.path.exists("test-1M.hm"):
     os.remove("test-1M.hm")
+if os.path.exists("test-1M-byChr.hm"):
+    os.remove("test-1M-byChr.hm")
+
 workingGenome = "hg19"
 genomeFolder = "../../fasta/hg19"
 if not os.path.exists(genomeFolder):
     setExceptionHook()
     print "----> saving by chromosome heatmap"
     TR.saveByChromosomeHeatmap(
-        filename[1] + "-1M.hm", resolution=1000000, includeTrans=True,
+        filename[1] + "-1M-byChr.hm", resolution=1000000, includeTrans=True,
         countDiagonalReads="once")
 
-    Tb = h5dict(filename[1] + "-1M.hm")["1 1"]
-    Tbb = h5dict(filename[1] + "-1M.hm")["1 2"]
+    Tb = h5dict(filename[1] + "-1M-byChr.hm")["1 1"]
+    Tbb = h5dict(filename[1] + "-1M-byChr.hm")["1 2"]
     assert ((Tb - chrom1) == 0).all()
     assert ((Tbb - chrom12) == 0).all()
     assert ((Tb + np.diag(np.diag(Tb))) == b).all()
 #os.remove("test_breaks.frag")
 os.remove("test-hg19.hdf5_parsed.frag")
 os.remove("test_merged.frag")
-os.remove("test-1M.hm")
+#os.remove("test-1M.hm")
 #os.remove("test_refined.frag")
 print "Test finished successfully!"