add hicdiff

coolbox/core/track.py

@@ -10,7 +10,7 @@
     "BigWig", "ABCompartment", "BedGraph",
     "Arcs", "BEDPE", "Pairs",
     "HiCMat", "Cool", "DotHiC",
-    "HicCompare", "Virtual4C",
+    "Virtual4C", "HiCDiff",
     "Ideogram", "GTF", "BAM"
 ]
 
@@ -823,7 +823,7 @@ def HiCMat(file_, *args, **kwargs):
 
 
 class HicCompare(Track, PlotHicCompare):
-    """
+    """Not in use for now!
     Track for express the comparison between two HiC Track.
 
     Parameters
@@ -1116,6 +1116,89 @@ def __init__(self, file_, **kwargs):
         super().__init__(properties_dict)
 
 
+class HiCDiff(Track, FetchHiCDiff, PlotHiCDiff):
+    """
+    Track for express the comparison between two HiC Track.
+
+    Parameters
+    ----------
+    hic1 : coolbox.api.track.Cool
+        First HiC Track or hic file path(.cool, .mcool, .hic).
+
+    hic2 : coolbox.api.track.Cool
+        Second HiC Track or hic file path(.cool, .mcool, .hic).
+
+    args_hic : dict, optional
+        Argument to create Hi-C instance, only in use
+        when first or second argument is a path.
+
+    style : {'triangular', 'window', 'matrix'}, optional
+        Matrix style, default 'triangular'.
+
+    depth_ratio : float, optional
+        Depth ratio of triangular matrix, use 'full' for full depth. default 'full'.
+
+    orientation : str, optional
+        Track orientation, use 'inverted' for inverted track plot.
+
+    normalize : str
+        Normalization method ('none', 'zscore', 'total', 'expect'), default 'expect'
+
+    diff_method : str
+        Difference method ('diff', 'log2fc'), default 'diff'
+
+    resolution : int, str
+        Resolution of sub two sample. default 'auto'
+
+    cmap : {str, matplotlib.colors.Colormap}, optional
+        A diverging colormap, positive color represent the first HiC file,
+        and negative represent the second HiC file.
+
+    color_bar : bool, optional
+        Show color bar or not.
+
+    max_value : {float, 'auto'}, optional
+        Max value of hic matrix, use 'auto' for specify max value automatically, default 'auto'.
+
+    min_value : {float, 'auto'}, optional
+        Min value of hic matrix, use 'auto' for specify min value automatically, default 'auto'.
+
+    title : str, optional
+        Label text, default ''.
+
+    name : str, optional
+        Track's name
+
+    """
+
+    DEFAULT_COLOR = "RdYlBu"
+
+    def __init__(self, hic1, hic2, args_hic=None, **kwargs):
+        args_hic = args_hic or {}
+        if isinstance(hic1, str):
+            hic1 = HiCMat(hic1, **args_hic)
+        if isinstance(hic2, str):
+            hic2 = HiCMat(hic2, **args_hic)
+        properties_dict = {
+            "hic1": hic1,
+            "hic2": hic2,
+            "resolution": "auto",
+            "normalize": "expect",
+            "diff_method": "diff",
+            "style": "triangular",
+            "depth_ratio": "full",
+            "cmap": HiCDiff.DEFAULT_COLOR,
+            "color_bar": True,
+            "max_value": "auto",
+            "min_value": "auto",
+            "title": '',
+        }
+        properties_dict.update(kwargs)
+        properties_dict['color'] = properties_dict['cmap']  # change key word
+
+        super().__init__(properties_dict)
+
+
 if __name__ == "__main__":
     import doctest
     doctest.testmod()

coolbox/fetchdata/__init__.py

@@ -16,14 +16,16 @@
 from .v4c import FetchVirtual4C
 from .gtf import FetchGTF
 from .bam import FetchBAM
+from .hicdiff import FetchHiCDiff
 
 
 __all__ = [
     "FetchBigWig", "FetchBedGraph", "FetchBed",
     "FetchBEDPE", "FetchPairs",
     "FetchCool", "FetchDotHiC",
     "FetchVirtual4C", "FetchGTF",
-    "FetchFrame", "FetchBAM"
+    "FetchFrame", "FetchBAM",
+    "FetchHiCDiff",
 ]
 
 

coolbox/fetchdata/hicdiff.py

@@ -0,0 +1,59 @@
+import numpy as np
+from scipy.linalg import toeplitz
+
+from .base import FetchTrackData
+
+
+class FetchHiCDiff(FetchTrackData):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def fetch_related_tracks(self, genome_range, resolution=None):
+        if resolution:
+            reso = resolution
+        else:
+            reso = self.properties['resolution']
+        hic1 = self.properties['hic1']
+        hic2 = self.properties['hic2']
+        mat1 = hic1.fetch_matrix(genome_range, reso)
+        mat2 = hic2.fetch_matrix(genome_range, reso)
+        return mat1, mat2
+
+    def __normalize_data(self, mat):
+        norm_mth = self.properties['normalize']
+        res = mat
+        if norm_mth == 'total':
+            total = np.sum(mat)
+            if total != 0:
+                res = mat / total
+        elif norm_mth == 'expect':
+            means = [np.diagonal(mat, i).mean() for i in range(mat.shape[0])]
+            expect = toeplitz(means)
+            res = mat / expect
+        elif norm_mth == 'zscore':
+            means = []
+            stds = []
+            for i in range(mat.shape[0]):
+                diagonal = np.diagonal(mat, i)
+                means.append(diagonal.mean())
+                stds.append(diagonal.std())
+            stds = np.array(stds)
+            stds[stds == 0] = stds[stds > 0].min()
+            mat_mean = toeplitz(means)
+            mat_std = toeplitz(stds)
+            res = (mat - mat_mean) / mat_std
+        return res
+
+    def __diff_data(self, mat1, mat2):
+        diff_mth = self.properties['diff_method']
+        if diff_mth == 'log2fc':
+            return np.log2((mat1 + 1)/(mat2 + 1))
+        else:
+            return mat1 - mat2
+
+    def fetch_data(self, genome_range, resolution=None):
+        mat1, mat2 = self.fetch_related_tracks(genome_range, resolution)
+        mat1, mat2 = self.__normalize_data(mat1), self.__normalize_data(mat2)
+        diff = self.__diff_data(mat1, mat2)
+        return diff

coolbox/plots/track/__init__.py

@@ -5,6 +5,7 @@
 from .arcs import PlotBEDPE, PlotPairs
 from coolbox.plots.track.hicmatrix import PlotCool, PlotDotHiC
 from .hiccompare import PlotHicCompare
+from .hicdiff import PlotHiCDiff
 from .virtual4C import PlotVirtual4C
 from .ideogram import PlotIdeogram
 from .gtf import PlotGTF

coolbox/plots/track/hiccompare.py

@@ -46,9 +46,9 @@ def plot(self, ax, chrom_region, start_region, end_region):
         self.ax = ax
 
         genome_range = GenomeRange(chrom_region, start_region, end_region)
-        arr1 = self.hic1.__fetch_matrix(genome_range)
+        arr1 = self.hic1.fetch_matrix(genome_range)
         self.hic1.matrix = arr1
-        arr2 = self.hic2.__fetch_matrix(genome_range)
+        arr2 = self.hic2.fetch_matrix(genome_range)
         self.hic2.matrix = arr2
         self.matrix = np.triu(arr1 * (-1), 1) + np.tril(arr2, -1)
 
@@ -110,7 +110,7 @@ def __plot_colorbar(self, img):
         cax = ax_divider.append_axes("bottom", size=0.09, pad=0.2)
         plt.colorbar(img, cax=cax, orientation='horizontal')
 
-    def get_tracks_height(self, frame_width):
+    def get_track_height(self, frame_width):
         """
         calculate track height dynamically.
         """

coolbox/plots/track/hicdiff.py

@@ -0,0 +1,18 @@
+from .hicmatrix import PlotHiCMatrix
+
+
+class PlotHiCDiff(PlotHiCMatrix):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.properties['transform'] = 'no'
+        self.properties['norm'] = 'no'
+
+    def fetch_matrix(self, genome_range, resolution='auto'):
+        diff = self.fetch_data(genome_range, None)
+        try:
+            self.small_value = diff[diff > 0].min()
+        except:
+            self.small_value = 1e-12
+        return diff
+

coolbox/plots/track/hicmatrix.py

@@ -7,7 +7,6 @@
 
 from coolbox.utilities import (
     GenomeRange,
-    change_chrom_names,
     get_logger
 )
 
@@ -39,9 +38,6 @@ def __init__(self, *args, **kwargs):
         self.matrix = None
         self._out_of_bound = False
 
-        from coolbox.utilities.hic.tools import file_type
-        self.file_type = file_type(self.properties['file'])
-
         self.fetched_binsize = None
 
     def __set_default_properties(self):
@@ -85,7 +81,9 @@ def balance(self):
         if self.properties['balance'] == 'no':
             return False
         else:
-            if self.file_type == '.hic':
+            file = self.properties['file']
+            from coolbox.utilities.hic.tools import file_type
+            if file_type(file) == '.hic':
                 if self.properties['balance'] == 'yes':
                     return 'KR'  # default use KR balance
                 else:
@@ -133,7 +131,7 @@ def matrix_val_range(self):
 
         return min_, max_
 
-    def __fetch_matrix(self, genome_range, resolution='auto'):
+    def fetch_matrix(self, genome_range, resolution='auto'):
         """
         Fetch the matrix.
 
@@ -150,7 +148,9 @@ def __fetch_matrix(self, genome_range, resolution='auto'):
         from coolbox.utilities.hic.wrap import StrawWrap, CoolerWrap
 
         path = self.properties['file']
-        if self.file_type == '.hic':
+
+        from coolbox.utilities.hic.tools import file_type
+        if file_type(self.properties['file']) == '.hic':
             wrap = StrawWrap(path, normalization=self.balance, binsize=resolution)
         else:
             wrap = CoolerWrap(path, balance=self.balance, binsize=resolution)
@@ -321,34 +321,30 @@ def __fetch_window_matrix(self, genome_range):
             self._out_of_bound = 'left'
 
         try:
-            arr = self.__fetch_matrix(fetch_range)
+            arr = self.fetch_matrix(fetch_range)
         except ValueError as e:
             if self._out_of_bound == 'left':
                 self._out_of_bound = 'both'
-                arr = self.__fetch_matrix(genome_range)
+                arr = self.fetch_matrix(genome_range)
             else:
                 self._out_of_bound = 'right'
                 fetch_range.end = genome_range.end
-                arr = self.__fetch_matrix(fetch_range)
+                arr = self.fetch_matrix(fetch_range)
         return arr, fetch_range
 
     def plot(self, ax, chrom_region, start_region, end_region):
         self.ax = ax
 
         self._out_of_bound = False
 
-        log.debug("plotting {}".format(self.properties['file']))
-
         genome_range = GenomeRange(chrom_region, start_region, end_region)
 
-        self.ax = ax
-
         # fetch matrix and perform transform process
         if self.style == STYLE_WINDOW:
             arr, fetch_region = self.__fetch_window_matrix(genome_range)
             self.fetch_region = fetch_region
         else:
-            arr = self.__fetch_matrix(genome_range)
+            arr = self.fetch_matrix(genome_range)
 
         self.matrix = arr
 

docker/readme.md

@@ -0,0 +1,13 @@
+Usage:
+
+see https://gangcaolab.github.io/CoolBox/installation.html#docker
+
+Build image:
+
+```bash
+$ docker build . -t nanguage/coolbox:${version}
+```
+
+DockerHub: 
+
+https://hub.docker.com/repository/docker/nanguage/coolbox