Fix chromsizes in scaling

pairtools/cli/scaling.py

@@ -53,7 +53,9 @@
     help="Number of bins to split the distance range in log10-space, specified per a factor of 10 difference.",
 )
 @common_io_options
-def scaling(input_path, output, view, chunksize, dist_range, n_dist_bins_decade, **kwargs):
+def scaling(
+    input_path, output, view, chunksize, dist_range, n_dist_bins_decade, **kwargs
+):
     """Calculate pairs scalings.
 
     INPUT_PATH : by default, a .pairs/.pairsam file to calculate statistics.
@@ -63,10 +65,14 @@ def scaling(input_path, output, view, chunksize, dist_range, n_dist_bins_decade,
 
     Output is .tsv file with scaling stats (both cis scalings and trans levels).
     """
-    scaling_py(input_path, output, view, chunksize, dist_range, n_dist_bins_decade, **kwargs)
+    scaling_py(
+        input_path, output, view, chunksize, dist_range, n_dist_bins_decade, **kwargs
+    )
 
 
-def scaling_py(input_path, output, view, chunksize, dist_range, n_dist_bins_decade, **kwargs):
+def scaling_py(
+    input_path, output, view, chunksize, dist_range, n_dist_bins_decade, **kwargs
+):
 
     if len(input_path) == 0:
         raise ValueError(f"No input paths: {input_path}")
@@ -87,11 +93,11 @@ def scaling_py(input_path, output, view, chunksize, dist_range, n_dist_bins_deca
     if view is not None:
         view = pd.read_table(view)
 
-    # Pass the header to the instream so that it can parse the header automatically
+    # chromsizes will be automatically extracted from the pairs file header
     cis_scalings, trans_levels = compute_scaling(
         instream,
         regions=view,
-        chromsizes=None,
+        chromsizes=None,  # Will be extracted from pairs header
         dist_range=dist_range,
         n_dist_bins_decade=n_dist_bins_decade,
         chunksize=chunksize,

pairtools/lib/scaling.py

@@ -52,9 +52,8 @@ def assign_regs(chroms, pos, regs):
 
     regs_dict = {
         chrom.encode(): regs_per_chrom[["start", "end"]]
-            .values
-            .flatten()
-            .astype(np.int64)
+        .values.flatten()
+        .astype(np.int64)
         for chrom, regs_per_chrom in gb_regs
     }
 
@@ -135,11 +134,17 @@ def make_empty_cross_region_table(
 
 
 def bins_pairs_by_distance(
-    pairs_df, dist_bins, regions=None, chromsizes=None, ignore_trans=False,
+    pairs_df,
+    dist_bins,
+    regions=None,
+    chromsizes=None,
+    ignore_trans=False,
     keep_unassigned=False,
 ):
 
     dist_bins = np.r_[dist_bins, np.iinfo(np.int64).max]
+    # Explicitly remove unmapped pair
+    pairs_df = pairs_df[(pairs_df["chrom1"] != "!") & (pairs_df["chrom2"] != "!")]
     if regions is None:
         if chromsizes is None:
             chroms = sorted(
@@ -188,7 +193,6 @@ def bins_pairs_by_distance(
             pairs_df.chrom2.values, pairs_df.pos2.values, regions
         ).T
 
-
     pairs_reduced_df = pd.DataFrame(
         {
             "chrom1": pairs_df.chrom1.values,
@@ -208,10 +212,11 @@ def bins_pairs_by_distance(
     )
 
     if not keep_unassigned:
-        pairs_reduced_df = (pairs_reduced_df
-            .query('(start1 >= 0) and (start2 >= 0)') 
+        pairs_reduced_df = (
+            pairs_reduced_df.query("(start1 >= 0) and (start2 >= 0)")
             # do not test for end1 and end2, as they can be -1 if regions and not specified
-            .reset_index(drop=True))
+            .reset_index(drop=True)
+        )
 
     pairs_reduced_df["min_dist"] = np.where(
         pairs_reduced_df["dist_bin_idx"] > 0,
@@ -220,7 +225,7 @@ def bins_pairs_by_distance(
     )
 
     pairs_reduced_df["max_dist"] = np.where(
-        pairs_reduced_df["dist_bin_idx"] < len(dist_bins)-1,
+        pairs_reduced_df["dist_bin_idx"] < len(dist_bins) - 1,
         dist_bins[pairs_reduced_df["dist_bin_idx"]],
         np.iinfo(np.int64).max,
     )
@@ -349,7 +354,7 @@ def compute_scaling(
     Parameters
     ----------
     pairs : pd.DataFrame or str or file-like object
-        A table with pairs of genomic coordinates representing contacts. 
+        A table with pairs of genomic coordinates representing contacts.
         It can be a pandas DataFrame, a path to a pairs file, or a file-like object.
     regions : bioframe viewframe or None, optional
         Genomic regions of interest. It can be anything that can serve as input to bioframe.from_any,
@@ -380,16 +385,21 @@ def compute_scaling(
     """
 
     dist_bins = geomspace(
-        dist_range[0], 
+        dist_range[0],
         dist_range[1],
-        int(np.round(np.log10(dist_range[1]/dist_range[0])*n_dist_bins_decade))
+        int(np.round(np.log10(dist_range[1] / dist_range[0]) * n_dist_bins_decade)),
     )
 
     if isinstance(pairs, pd.DataFrame):
         pairs_df = pairs
 
     elif isinstance(pairs, str) or hasattr(pairs, "buffer") or hasattr(pairs, "peek"):
-        pairs_df, _, _ = pairsio.read_pairs(pairs, nproc=nproc_in, chunksize=chunksize)
+        pairs_df, _, chromsizes_extracted = pairsio.read_pairs(
+            pairs, nproc=nproc_in, chunksize=chunksize
+        )
+        # Use chromsizes from the header if not provided explicitly
+        if chromsizes is None:
+            chromsizes = chromsizes_extracted
     else:
         raise ValueError(
             "pairs must be either a path to a pairs file or a pd.DataFrame"
@@ -405,7 +415,7 @@ def compute_scaling(
             regions=regions,
             chromsizes=chromsizes,
             ignore_trans=ignore_trans,
-            keep_unassigned=keep_unassigned
+            keep_unassigned=keep_unassigned,
         )
 
         sc = sc_chunk if sc is None else sc.add(sc_chunk, fill_value=0)
@@ -416,7 +426,6 @@ def compute_scaling(
             else trans_counts.add(trans_counts_chunk, fill_value=0)
         )
 
-
     #         if not (isinstance(regions, pd.DataFrame) and
     #                  (set(regions.columns) == set(['chrom', 'start','end']))):
     #             raise ValueError('regions must be provided as a dict or chrom-indexed Series of chromsizes or as a bedframe.')
@@ -428,10 +437,9 @@ def compute_scaling(
 
     if not ignore_trans:
         trans_counts.reset_index(inplace=True)
-        trans_counts["n_bp2"] = (
-            (trans_counts["end1"] - trans_counts["start1"]) * (
+        trans_counts["n_bp2"] = (trans_counts["end1"] - trans_counts["start1"]) * (
             trans_counts["end2"] - trans_counts["start2"]
-        ))
+        )
 
     return sc, trans_counts
 
@@ -451,12 +459,12 @@ def norm_scaling_factor(bins, cfreqs, norm_window):
     """
 
     lo, hi = np.searchsorted(bins, norm_window)
-    return cfreqs[lo:hi+1].mean()
+    return cfreqs[lo : hi + 1].mean()
 
 
 def norm_scaling(bins, cfreqs, norm_window, log_input=False):
     """
-    Normalize a contact-frequency-vs-distance curve, by setting the average contact frequency 
+    Normalize a contact-frequency-vs-distance curve, by setting the average contact frequency
     in a given window to 1.0.
 
     Args:
@@ -468,14 +476,14 @@ def norm_scaling(bins, cfreqs, norm_window, log_input=False):
     Returns:
         float or array-like: The normalized contact frequencies.
     """
-    
+
     norm = norm_scaling_factor(bins, cfreqs, norm_window)
     if log_input:
         return cfreqs - norm
     else:
         return cfreqs / norm
 
-    
+
 def unity_norm_scaling(bins, cfreqs, norm_range=(1e4, 1e9)):
     bin_lens = np.diff(bins)
     bin_mids = np.sqrt(bins[1:] * bins[:-1])

tests/test_scaling.py

@@ -22,6 +22,4 @@ def test_scaling():
 
     output = pd.read_csv(io.StringIO(result), sep="\t", header=0)
 
-    assert (
-        output["n_pairs"].sum() == 9
-    )  # double unmapped pairs are currently ignored by lib.scaling
+    assert output["n_pairs"].sum() == 7  # unmapped pairs are ignored by lib.scaling