Merge pull request #594 from deeptools/develop

Develop

Author: joachimwolff

hicexplorer/_version.py

@@ -2,4 +2,4 @@
 # This file is originally generated from Git information by running 'setup.py
 # version'. Distribution tarballs contain a pre-generated copy of this file.
 
-__version__ = '3.5'
+__version__ = '3.5.1'

hicexplorer/hicAdjustMatrix.py

@@ -137,6 +137,7 @@ def adjustMatrix(pArgs):
         elif pArgs.action == 'remove':
             hic_matrix.maskBins(matrix_indices_regions)
             hic_matrix.orig_bin_ids = []
+            hic_matrix.orig_cut_intervals = []
             hic_matrix.nan_bins = []
 
     elif pArgs.maskBadRegions:

hicexplorer/hicBuildMatrix.py

@@ -119,21 +119,19 @@ def parse_arguments(args=None):
                                 metavar='FOLDER',
                                 required=True)
     parserRequired.add_argument('--restrictionCutFile', '-rs',
-                                help=('BED file(s) with all restriction cut places '
-                                      '(output of "findRestSite" command). '
-                                      'Should contain only  mappable '
-                                      'restriction sites. If given, the bins are '
-                                      'set to match the restriction fragments (i.e. '
-                                      'the region between one restriction site and '
-                                      'the next). Alternatively, a fixed binSize can be defined instead. '
-                                      'However, either binSize or restrictionCutFile must be defined. '
-                                      'To use more than one restriction enzyme, generate for each one a restrictionCutFile and list them space seperated.'),
+                                help='BED file(s) with all restriction cut sites '
+                                '(output of "hicFindRestSite" command). '
+                                'Should only contain the restriction sites of the same genome which has been used '
+                                'to generate the input sam files. Using regions of a different genome version can '
+                                'generate false results! To use more than one restriction enzyme, generate '
+                                'a restrictionCutFile for each enzyne and list them space seperated.',
                                 type=argparse.FileType('r'),
                                 metavar='BED file',
                                 nargs='+',
                                 required=True)
     parserRequired.add_argument('--restrictionSequence', '-seq',
-                                help='Sequence of the restriction site, if multiple are used, please list them space seperated. If a dangling sequence '
+                                help='Sequence of the restriction site, if multiple are used, '
+                                'please list them space seperated. If a dangling sequence '
                                 'is listed at the same time, please preserve the same order.',
                                 type=str,
                                 nargs='+',
@@ -172,9 +170,8 @@ def parse_arguments(args=None):
     parserOpt.add_argument('--binSize', '-bs',
                            help='Size in bp for the bins. The bin size depends '
                            'on the depth of sequencing. Use a larger bin size for '
-                           'libraries sequenced with lower depth. Alternatively, the location of '
-                           'the restriction sites can be given (see --restrictionCutFile). However, either binSize or restrictionCutFile must be defined.'
-                           'Optional for mcool file format: Define multiple resolutions which are all a multiple of the first value. '
+                           'libraries sequenced with lower depth. If not given, matrices of restriction site resolution will be built. '
+                           'Optionally for mcool file format: Define multiple resolutions which are all a multiple of the first value. '
                            ' Example: --binSize 10000 20000 50000 will create a mcool file formate containing the three defined resolutions.',
                            type=int,
                            nargs='+')

hicexplorer/hicCorrelate.py

@@ -325,7 +325,7 @@ def main(args=None):
         hic_mat_list.append(_mat)
 
     # remove nan bins
-    rows_keep = cols_keep = np.delete(list(range(all_mat.shape[1])), all_nan)
+    rows_keep = cols_keep = np.delete(list(range(all_mat.shape[1])), all_nan.astype('int'))
     all_mat = all_mat[rows_keep, :][:, cols_keep]
 
     # make large matrix to correlate by

hicexplorer/hicMergeDomains.py

@@ -6,6 +6,7 @@
 import os
 import logging
 log = logging.getLogger(__name__)
+from hicexplorer._version import __version__
 
 
 def parse_arguments(args=None):
@@ -72,6 +73,9 @@ def parse_arguments(args=None):
                            help='File format of the relationship tree. Supported formats are listed on: https://www.graphviz.org/doc/info/output.html',
                            default='pdf',
                            required=False)
+    parserOpt.add_argument("--help", "-h", action="help", help="show this help message and exit")
+
+    parserOpt.add_argument('--version', action='version', version='%(prog)s {}'.format(__version__))
 
     return parser
 

hicexplorer/hicPlotMatrix.py

@@ -390,15 +390,14 @@ def plotHeatmap(ma, chrBinBoundaries, fig, position, args, cmap, xlabel=None,
                            pValueMin=args.vMinBigwig, pValueMax=args.vMaxBigwig, pResolution=pResolution)
 
 
-def translate_region(region_string):
+def translate_region(region_string, ma):
     """
     Takes an string and returns a list
     of chrom, start, end.
     If the region string only contains
     the chrom, then start and end
-    are set to a 0 and 1e15
+    are set to matrix start and end on this chromosome
     """
-
     # region_string = toBytes(region_string)
     region_string = region_string.replace(",", "")
     region_string = region_string.replace(";", "")
@@ -407,14 +406,24 @@ def translate_region(region_string):
 
     fields = region_string.split(":")
     chrom = fields[0]
+
+    chrom = check_chrom_str_bytes(ma.interval_trees, chrom)
+    if chrom not in list(ma.interval_trees):
+        chrom = change_chrom_names(chrom)
+        chrom = check_chrom_str_bytes(ma.interval_trees, chrom)
+        if chrom not in list(ma.interval_trees):
+            exit(
+                "Chromosome name {} in --region not in matrix".format(change_chrom_names(chrom)))
+
+    first_bin, last_bin = ma.getChrBinRange(chrom)
     try:
         region_start = int(fields[1])
     except IndexError:
-        region_start = 0
+        region_start = ma.getBinPos(first_bin)[1]
     try:
         region_end = int(fields[2])
     except IndexError:
-        region_end = 1e15  # vert large number
+        region_end = ma.getBinPos(last_bin - 1)[2]
 
     return chrom, region_start, region_end
 
@@ -559,43 +568,27 @@ def plotPerChr(hic_matrix, cmap, args, pBigwig, pResolution):
             args, hic_matrix)
         plotHeatmap(matrix, chr_bin_boundary, fig, None,
                     args, cmap, xlabel=chrname, ylabel=chrname,
-                    start_pos=start_pos1, start_pos2=start_pos2, pNorm=norm, pAxis=axis, pBigwig=bigwig_info, pChromsomeStartEndDict=chromosome_start_end(hic_matrix), pResolution=pResolution)
+                    start_pos=start_pos1, start_pos2=start_pos2, pNorm=norm, pAxis=axis, pBigwig=bigwig_info,
+                    pChromsomeStartEndDict=chromosome_start_end(hic_matrix), pResolution=pResolution)
     return fig
 
 
 def getRegion(args, ma):
     chrom = region_start = region_end = idx1 = start_pos1 = chrom2 = region_start2 = region_end2 = idx2 = start_pos2 = None
-    chrom, region_start, region_end = translate_region(args.region)
-
-    chrom = check_chrom_str_bytes(ma.interval_trees, chrom)
-    # if type(next(iter(ma.interval_trees))) in [np.bytes_, bytes]:
-    #     chrom = toBytes(chrom)
-
-    if chrom not in list(ma.interval_trees):
-
-        chrom = change_chrom_names(chrom)
-
-        chrom = check_chrom_str_bytes(ma.interval_trees, chrom)
-
-        # if type(next(iter(ma.interval_trees))) in [np.bytes_, bytes]:
-        #     chrom = toBytes(chrom)
-
-        if chrom not in list(ma.interval_trees):
-            exit(
-                "Chromosome name {} in --region not in matrix".format(change_chrom_names(chrom)))
+    chrom, region_start, region_end = translate_region(args.region, ma)
 
     args.region = [chrom, region_start, region_end]
     is_cooler = check_cooler(args.matrix)
     if is_cooler:
         idx1, start_pos1 = zip(*[(idx, x[1]) for idx, x in enumerate(ma.cut_intervals) if x[0] == chrom and  # noqa: W504
-                                 ((x[1] >= region_start and x[2] < region_end) or                            # noqa: W504
+                                 ((x[1] >= region_start and x[2] <= region_end) or                            # noqa: W504
                                   (x[1] < region_end and x[2] < region_end and x[2] > region_start) or       # noqa: W504
                                   (x[1] > region_start and x[1] < region_end))])
     else:
         idx1, start_pos1 = zip(*[(idx, x[1]) for idx, x in enumerate(ma.cut_intervals) if x[0] == chrom and  # noqa: W504
-                                 x[1] >= region_start and x[2] < region_end])
+                                 x[1] >= region_start and x[2] <= region_end])
     if hasattr(args, 'region2') and args.region2:
-        chrom2, region_start2, region_end2 = translate_region(args.region2)
+        chrom2, region_start2, region_end2 = translate_region(args.region2, ma)
         chrom2 = check_chrom_str_bytes(ma.interval_trees, chrom2)
 
         # if type(next(iter(ma.interval_trees))) in [np.bytes_, bytes]:
@@ -780,7 +773,6 @@ def main(args=None):
             matrix = np.asarray(ma.getMatrix().astype(float))
 
     resolution = ma.getBinSize()
-
     matrix_length = len(matrix[0])
     log.debug("Number of data points matrix: {}".format(matrix_length))
 

hicexplorer/test/general/test_hicAggregateContacts.py

@@ -25,7 +25,7 @@
 
 
 ROOT = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "test_data/")
-tolerance = 13  # default matplotlib pixed difference tolerance
+tolerance = 50  # default matplotlib pixed difference tolerance
 
 
 @pytest.mark.xfail(raises=ImageComparisonFailure, reason='Matplotlib plots for reasons a different image size.')

hicexplorer/test/general/test_hicCorrelate.py

@@ -22,8 +22,8 @@ def test_correlate():
         "--outFileNameHeatmap {} --outFileNameScatter {}".format(ROOT + "hicCorrectMatrix/small_test_matrix_ICEcorrected_chrUextra_chr3LHet.h5",
                                                                  ROOT + "hicCorrectMatrix/small_test_matrix_ICEcorrected_chrUextra_chr3LHet.h5",
                                                                  outfile_heatmap.name, outfile_scatter.name).split()
-    # hicCorrelate.main(args)
-    compute(hicCorrelate.main, args, 5)
+    hicCorrelate.main(args)
+    # compute(hicCorrelate.main, args, 5)
     res = compare_images(ROOT + "hicCorrelate" + '/heatmap.png', outfile_heatmap.name, tol=40)
     assert res is None, res
 
@@ -44,8 +44,8 @@ def test_correlate_chromosomes():
         "--chromosomes chrUextra chr3LHet".format(ROOT + "hicCorrectMatrix/small_test_matrix_ICEcorrected_chrUextra_chr3LHet.h5",
                                                   ROOT + "hicCorrectMatrix/small_test_matrix_ICEcorrected_chrUextra_chr3LHet.h5",
                                                   outfile_heatmap.name, outfile_scatter.name).split()
-    # hicCorrelate.main(args)
-    compute(hicCorrelate.main, args, 5)
+    hicCorrelate.main(args)
+    # compute(hicCorrelate.main, args, 5)
 
     res = compare_images(ROOT + "hicCorrelate" + '/heatmap_chrom.png', outfile_heatmap.name, tol=40)
     assert res is None, res

hicexplorer/test/general/test_hicPlotMatrix.py

@@ -440,7 +440,7 @@ def test_hicPlotMatrix_cool_perChr_log1p():
 
     outfile = NamedTemporaryFile(suffix='.png', prefix='hicexplorer_test', delete=False)
 
-    args = "--matrix {0}/small_test_matrix_50kb_res.h5 --perChr  --disable_tight_layout " \
+    args = "--matrix {0}/small_test_matrix_50kb_res.cool --perChr  --disable_tight_layout " \
            "--outFileName  {1} --log1 --vMax 10 ".format(ROOT, outfile.name).split()
     # hicexplorer.hicPlotMatrix.main(args)
     compute(hicexplorer.hicPlotMatrix.main, args, 5)

requirements.txt

@@ -11,12 +11,12 @@ pandas >= 0.25.*
 pybigwig >= 0.3.*
 jinja2 >= 2.10
 unidecode >= 1.1.*
-hicmatrix >= 13
+hicmatrix >= 14
 hic2cool >= 0.8.3
 psutil >= 5.6.*
 pygenometracks >= 3.2
 fit_nbinom >= 1.1
-cooler >= 0.8.5
+cooler >= 0.8.9
 krbalancing >= 0.0.5
 pybedtools >= 0.8.*
 future >= 0.17