Merge pull request #312 from deeptools/expectedMatrix

Expected matrix

README.rst

@@ -31,6 +31,8 @@ Citation:
 Fidel Ramirez, Vivek Bhardwaj, Jose Villaveces, Laura Arrigoni, Bjoern A Gruening, Kin Chung Lam, Bianca Habermann, Asifa Akhtar, Thomas Manke.
 **"High-resolution TADs reveal DNA sequences underlying genome organization in flies". Nature Communications**, Volume 9, Article number: 189 (2018), doi: https://doi.org/10.1038/s41467-017-02525-w
 
+Joachim Wolff, Vivek Bhardwaj, Stephan Nothjunge, Gautier Richard, Gina Renschler, Ralf Gilsbach, Thomas Manke, Rolf Backofen, Fidel Ramírez, Björn A Grüning. 
+**"Galaxy HiCExplorer: a web server for reproducible Hi-C data analysis, quality control and visualization", Nucleic Acids Research**, Volume 46, Issue W1, 2 July 2018, Pages W11–W16, doi: https://doi.org/10.1093/nar/gky504
 
 .. image:: ./docs/images/hicex2.png
 

docs/content/example_usage.rst

@@ -219,7 +219,7 @@ the covariance matrix are computed. All these steps are computed with the comman
 
 .. code-block:: bash
 
-   $ hicPCA -m hic_corrected.h5 --outFileName pca1.bedgraph pca2.bedgraph 
+   $ hicPCA -m hic_corrected.h5 --outFileName pca1.bw pca2.bw --format bigwig 
 
 If the intermediate matrices of this process should be used for plotting run:
 
@@ -233,7 +233,7 @@ The A / B compartments can be plotted with :ref:`hicPlotMatrix`.
 
 .. code-block:: bash
 
-   $ hicPlotMatrix -m pearson_all.h5 --outFileName pca1.png --perChr --pca pca1.bedgraph
+   $ hicPlotMatrix -m pearson_all.h5 --outFileName pca1.png --perChr --bigwig pca1.bw
 
 //.. figure:: ../images/eigenvector1_lieberman.png
 //    :scale: 90 %

docs/index.rst

@@ -40,13 +40,12 @@ tool                            description
 :ref:`hicPlotMatrix`             Plots a Hi-C matrix as a heatmap
 :ref:`hicPlotTADs`               Plots TADs as a track that can be combined with other tracks (genes, signal, interactions)
 :ref:`hicPlotViewpoint`          A plot with the interactions around a reference point or region.
-:ref:`hicAggreagteContacts`      A tool that allows plotting of aggregated Hi-C sub-matrices of a specified list of positions.
+:ref:`hicAggregateContacts`      A tool that allows plotting of aggregated Hi-C sub-matrices of a specified list of positions.
 :ref:`hicSumMatrices`            Adds Hi-C matrices of the same size
 :ref:`hicPlotDistVsCounts`       Plots distance vs. Hi-C counts of corrected data
 :ref:`hicExport`                 Export matrix to text formats
 :ref:`hicInfo`                   Shows information about a Hi-C matrix file (no. of bins, bin length, sum, max, min, etc)
 :ref:`hicCompareMatrices`        Computes difference or ratio between two matrices
-:ref:`hicLog2Ratio`              Computes the log2 ratio between two matrices.
 =============================== ==========================================================================================================================================================
 
 

hicexplorer/hicExpectedMatrix.py

@@ -0,0 +1,62 @@
+from __future__ import division
+
+import sys
+import argparse
+import numpy as np
+import os
+from unidecode import unidecode
+from hicmatrix import HiCMatrix as hm
+import warnings
+
+def parse_arguments():
+    parser = argparse.ArgumentParser(description = "Making the expected matrix")
+    #Required
+    parser.add_argument("-m",
+                        dest = "count_matrix",
+                        required = True)
+    parser.add_argument("-o",
+                        dest = "output_matrix",
+                        required = True)
+    return parser
+def expected_interactions(matrix):
+    """
+        Computes the expected number of interactions per distance
+    """
+
+    expected_interactions = np.zeros(matrix.shape[0])
+    row, col = matrix.nonzero()
+    distance = np.absolute(row - col)
+    occurences = np.arange(matrix.shape[0], 0, -1)
+
+    for i, distance_ in enumerate(distance):
+        expected_interactions[distance_] += matrix.data[i]
+    expected_interactions /= occurences
+
+    mask = np.isnan(expected_interactions)
+    expected_interactions[mask] = 0
+    mask = np.isinf(expected_interactions)
+    expected_interactions[mask] = 0
+    return expected_interactions
+
+
+def _expected_matrix(matrix):
+    expected_interactions_in_distance_ = expected_interactions(matrix)
+    row, col = matrix.nonzero()
+    distance = np.absolute(row - col).astype(np.int32)
+
+    if len(matrix.data) > 0:
+        data_type = type(matrix.data[0])
+        expected = expected_interactions_in_distance_[distance]
+        matrix.data = matrix.data.astype(np.float32)
+        matrix.data = expected
+        matrix.data = np.nan_to_num(matrix.data).astype(data_type)
+    return matrix
+
+
+def main():
+    parser = parse_arguments()
+    args = parser.parse_args()
+    ma = hm.hiCMatrix(args.count_matrix)
+    expected_matrix = _expected_matrix(ma.matrix)
+    ma.setMatrixValues(expected_matrix)
+    ma.save(args.output_matrix, pApplyCorrection=False)

hicexplorer/utilities.py

@@ -3,6 +3,8 @@
 import numpy as np
 import argparse
 from matplotlib import use as mplt_use
+from hicexplorer.hicExpectedMatrix import expected_interactions
+
 mplt_use('Agg')
 from unidecode import unidecode
 import cooler
@@ -260,28 +262,6 @@ def expected_interactions_non_zero(pSubmatrix):
     return expected_interactions
 
 
-def expected_interactions(pSubmatrix):
-    """
-        Computes the expected number of interactions per distance
-    """
-
-    expected_interactions = np.zeros(pSubmatrix.shape[0])
-    row, col = pSubmatrix.nonzero()
-    distance = np.absolute(row - col)
-    occurrences = np.arange(pSubmatrix.shape[0] + 1, 1, -1)
-    # occurences = np.zeros(pSubmatrix.shape[0])
-    for i, distance_ in enumerate(distance):
-        expected_interactions[distance_] += pSubmatrix.data[i]
-        # occurences[distance_] += 1
-    expected_interactions /= occurrences
-
-    mask = np.isnan(expected_interactions)
-    expected_interactions[mask] = 0
-    mask = np.isinf(expected_interactions)
-    expected_interactions[mask] = 0
-
-    return expected_interactions
-
 
 def obs_exp_matrix_lieberman(pSubmatrix, pLength_chromosome, pChromosome_count):
     """
@@ -366,18 +346,16 @@ def obs_exp_matrix(pSubmatrix):
         exp_i,j = sum(interactions at distance abs(i-j)) / number of non-zero interactions at abs(i-j)
 
     """
-
     expected_interactions_in_distance_ = expected_interactions(pSubmatrix)
     row, col = pSubmatrix.nonzero()
-    distance = np.ceil(np.absolute(row - col) / 2).astype(np.int32)
+    distance = np.absolute(row - col).astype(np.int32)
 
     if len(pSubmatrix.data) > 0:
-        data_type = type(pSubmatrix.data[0])
-
         expected = expected_interactions_in_distance_[distance]
         pSubmatrix.data = pSubmatrix.data.astype(np.float32)
         pSubmatrix.data /= expected
-        pSubmatrix.data = convertInfsToZeros_ArrayFloat(pSubmatrix.data).astype(data_type)
+        pSubmatrix.data = convertInfsToZeros_ArrayFloat(pSubmatrix.data)
+
     return pSubmatrix
 
 

setup.py

@@ -128,7 +128,7 @@ def checkProgramIsInstalled(self, program, args, where_to_download,
              'bin/hicPlotTADs', 'bin/hicSumMatrices', 'bin/hicExport', 'bin/hicInfo', 'bin/hicexplorer',
              'bin/hicQC', 'bin/hicCompareMatrices', 'bin/hicPCA', 'bin/hicTransform', 'bin/hicPlotViewpoint',
              'bin/hicConvertFormat', 'bin/hicAdjustMatrix', 'bin/hicNormalize',
-             'bin/hicAverageRegions', 'bin/hicPlotAverageRegions'
+             'bin/hicAverageRegions', 'bin/hicPlotAverageRegions', 'bin/hicExpectedMatrix'
              ],
     include_package_data=True,
     package_dir={'hicexplorer': 'hicexplorer'},