Merge pull request #38 from deeptools/develop

Develop

hicmatrix/HiCMatrix.py

@@ -50,7 +50,7 @@ def __init__(self, pMatrixFile=None, pChrnameList=None, pDistance=None, pNoInter
             fileType = 'cool'
             if pMatrixFile.endswith('.h5'):
                 fileType = 'h5'
-            self.matrixFileHandler = MatrixFileHandler(pFileType=fileType, pApplyCorrectionCoolerLoad=False, pMatrixFile=pMatrixFile, pChrnameList=pChrnameList, pDistance=pDistance, pMatrixFormat=pMatrixFormat, pLoadMatrixOnly=pLoadMatrixOnly)
+            self.matrixFileHandler = MatrixFileHandler(pFileType=fileType, pMatrixFile=pMatrixFile, pChrnameList=pChrnameList, pDistance=pDistance, pMatrixFormat=pMatrixFormat, pLoadMatrixOnly=pLoadMatrixOnly)
             log.debug('init time: {}'.format(time.time() - start_time))
             matrixFileHandler_load = self.matrixFileHandler.load()
             # check if there was any exception thrown in the load function
@@ -960,7 +960,25 @@ def printchrtoremove(self, to_remove, label="Number of poor regions to remove",
         log.info('{}: {} {}'.format(label, len(to_remove), cnt))
         self.prev_to_remove = to_remove
 
+    def get_chromosome_sizes_real(self):
+        '''
+        Function returns the size of a chromosome as it is stored in the  matrix.
+        The size can differ if e.g. some area from the start or end of a chromosome is not present in the interaction matrix.
+        '''
+        if self.chrBinBoundaries and len(self.chrBinBoundaries) > 0:
+            chrom_sizes = OrderedDict()
+            # for chrom, (start_bin, end_bin) in iteritems(self.chrBinBoundaries):
+            for chrom, (start_bin, end_bin) in self.chrBinBoundaries.items():
+                chrom, start0, end0, _ = self.cut_intervals[start_bin]
+                chrom, start1, end1, _ = self.cut_intervals[end_bin - 1]
+                chrom_sizes[chrom] = end1 - start0 + 1
+
+            return chrom_sizes
+
     def get_chromosome_sizes(self):
+        '''
+        Function returns the size of a chromosome as it is stored in the  matrix, assuming the chromosome starts is always at its genomic position 0.
+        '''
         if self.chrBinBoundaries and len(self.chrBinBoundaries) > 0:
             chrom_sizes = OrderedDict()
             # for chrom, (start_bin, end_bin) in iteritems(self.chrBinBoundaries):

hicmatrix/_version.py

@@ -1,2 +1,2 @@
-__version__ = '14'
+__version__ = '15'
 # Version number differs from HiCExplorer!

hicmatrix/test/test_HiCMatrix.py

@@ -43,6 +43,21 @@ def test_load_h5_save_and_load_cool():
     nt.assert_equal(end_cool, end)
 
 
+def test_load_h5_load_cool_weight():
+    hic_h5 = hm.hiCMatrix(ROOT + 'Li_et_al_2015.h5')
+    hic_cool = hm.hiCMatrix(ROOT + 'Li_et_al_2015.cool')
+
+    # there is always a small gap due to rounding errors and inaccurate floating operations
+    # test if it is equal for up to 10 decimal positions
+    nt.assert_almost_equal(hic_cool.matrix.data, hic_h5.matrix.data, decimal=10)
+    chrom_cool, start_cool, end_cool, _ = list(zip(*hic_cool.cut_intervals))
+    chrom, start, end, _ = list(zip(*hic_cool.cut_intervals))
+
+    nt.assert_equal(chrom_cool, chrom)
+    nt.assert_equal(start_cool, start)
+    nt.assert_equal(end_cool, end)
+
+
 def test_load_h5_save_and_load_cool_2():
     hic = hm.hiCMatrix(ROOT + 'small_test_matrix.h5')
 
@@ -291,7 +306,7 @@ def test_save():
     outfile_cool = NamedTemporaryFile(suffix='.cool', delete=False)
     outfile_cool.close()
 
-    outfile_h5 = NamedTemporaryFile(suffix='.cool', delete=False)
+    outfile_h5 = NamedTemporaryFile(suffix='.h5', delete=False)
     outfile_h5.close()
 
     hic = hm.hiCMatrix()
@@ -315,6 +330,7 @@ def test_save():
     h5_test = hm.hiCMatrix(outfile_h5.name)
 
     # test cool
+    hic.matrixFileHandler = None
     hic.save(outfile_cool.name)
     cool_test = hm.hiCMatrix(outfile_cool.name)
 
@@ -931,6 +947,45 @@ def test_printchrtoremove(capsys):
     nt.assert_equal(hic.prev_to_remove, np.array(to_remove))
 
 
+def test_get_chromosome_sizes_real():
+    # get matrix
+    hic = hm.hiCMatrix()
+    cut_intervals = [('a', 0, 10, 1), ('a', 10, 20, 1),
+                     ('a', 20, 30, 1), ('b', 30, 40, 1), ('b', 40, 50, 1)]
+
+    hic.nan_bins = []
+
+    matrix = np.array([[1, 8, 5, 3, 0],
+                       [0, 4, 15, 5, 1],
+                       [0, 0, 0, 0, 2],
+                       [0, 0, 0, 0, 1],
+                       [0, 0, 0, 0, 0]])
+
+    hic.matrix = csr_matrix(matrix)
+    hic.setMatrix(hic.matrix, cut_intervals)
+
+    nt.assert_equal(hic.getMatrix(), matrix)
+
+    # define expected outcome
+    expected_sizes = OrderedDict([('a', 31), ('b', 21)])
+
+    chrom_sizes = hic.get_chromosome_sizes_real()
+
+    nt.assert_equal(chrom_sizes, expected_sizes)
+
+    # define new intervals and test again
+    new_cut_intervals = [('a', 0, 10, 1), ('b', 10, 20, 1),
+                         ('b', 20, 30, 1), ('c', 30, 40, 1), ('c', 40, 90, 1)]
+
+    expected_sizes = OrderedDict([('a', 11), ('b', 21), ('c', 61)])
+
+    hic.setMatrix(hic.matrix, new_cut_intervals)
+
+    chrom_sizes = hic.get_chromosome_sizes_real()
+
+    nt.assert_equal(chrom_sizes, expected_sizes)
+
+
 def test_get_chromosome_sizes():
     # get matrix
     hic = hm.hiCMatrix()