Fix Hi-C bin coverage calculation across sparse batches

scripts/hicres.py

@@ -12,24 +12,37 @@ def calculate_valid_bins_percentage(file_path, threshold=1000, num_threads=1, ch
 
 def calculate_cool_valid_bins_percentage(file_path, threshold=1000, chunk_size=100000):
     clr = cooler.Cooler(file_path)
-    bins = clr.bins()[:]
-    total_valid_bins = 0
-    total_bins = 0
-
-    for start in range(0, len(bins), chunk_size):
-        end = min(start + chunk_size, len(bins))
-        matrix = clr.matrix(balance=False)[start:end, start:end]
-
-        bin_sums = np.sum(matrix, axis=0) + np.sum(matrix, axis=1) - np.diag(matrix)
-        valid_bins = np.sum(bin_sums > threshold)
-        total_valid_bins += valid_bins
-        total_bins += bin_sums.size
-
-    valid_bin_percentage = (total_valid_bins / total_bins) * 100 if total_bins > 0 else 0
-
     return {
         'File': file_path,
-        'Valid Bin Percentage': valid_bin_percentage,
+        **calculate_bin_statistics(clr, threshold, chunk_size)
+    }
+
+def iter_bin_contact_sums(clr, chunk_size=100000):
+    """Yield raw contact coverage, counting each matrix row exactly once.
+
+    Cooler reconstructs symmetric-upper storage automatically. Read every
+    column, including contacts outside the row batch, and count the diagonal
+    once. Bound internal batches to keep sparse selections manageable.
+    """
+    if chunk_size <= 0:
+        raise ValueError('chunk_size must be positive')
+    row_batch_size = min(chunk_size, 4096)
+    total_bins = clr.info['nbins']
+    selector = clr.matrix(balance=False, sparse=True)
+    for start in range(0, total_bins, row_batch_size):
+        end = min(start + row_batch_size, total_bins)
+        matrix = selector[start:end, :]
+        # Accumulate fractional input counts in float64, not float32.
+        yield np.asarray(matrix.sum(axis=1, dtype=np.float64)).ravel()
+
+def calculate_bin_statistics(clr, threshold=1000, chunk_size=100000):
+    total_bins = clr.info['nbins']
+    total_valid_bins = sum(
+        int(np.count_nonzero(bin_sums >= threshold))
+        for bin_sums in iter_bin_contact_sums(clr, chunk_size)
+    )
+    return {
+        'Valid Bin Percentage': (total_valid_bins / total_bins) * 100 if total_bins else 0,
         'Total Valid Bins': total_valid_bins,
         'Total Bins': total_bins
     }
@@ -70,29 +83,11 @@ def calculate_mcool_valid_bins_percentage(file_path, threshold=1000, num_threads
 
 def process_resolution(file_path, resolution_path, threshold, chunk_size):
     clr_res = cooler.Cooler(f"{file_path}::{resolution_path}")
-    bins = clr_res.bins()[:]
 
     print(f"Processing resolution: {resolution_path}")
-
-    total_valid_bins = 0
-    total_bins = 0
-
-    for start in range(0, len(bins), chunk_size):
-        end = min(start + chunk_size, len(bins))
-        matrix = clr_res.matrix(balance=False)[start:end, start:end]
-
-        bin_sums = np.sum(matrix, axis=0) + np.sum(matrix, axis=1) - np.diag(matrix)
-        valid_bins = np.sum(bin_sums > threshold)
-        total_valid_bins += valid_bins
-        total_bins += bin_sums.size
-
-    valid_bin_percentage = (total_valid_bins / total_bins) * 100 if total_bins > 0 else 0
-
     return {
         'Resolution': resolution_path.split('/')[-1],
-        'Valid Bin Percentage': valid_bin_percentage,
-        'Total Valid Bins': total_valid_bins,
-        'Total Bins': total_bins
+        **calculate_bin_statistics(clr_res, threshold, chunk_size)
     }
 
 def save_results_to_csv(results, output_path):
@@ -110,7 +105,7 @@ def main():
     parser.add_argument("file_path", type=str, help="Path to the .cool or .mcool file")
     parser.add_argument("--output", type=str, required=True, help="Path to save the results as CSV file")
     parser.add_argument("--threads", type=int, default=1, help="Number of threads to use for processing")
-    parser.add_argument("--chunk_size", type=int, default=100000, help="Chunk size for processing large matrices")
+    parser.add_argument("--chunk_size", type=int, default=100000, help="Maximum rows per sparse batch (internally capped at 4096; does not affect counts)")
     parser.add_argument("--resolutions", type=str, help="Comma-separated list of resolutions to process in the .mcool file")
 
     args = parser.parse_args()
@@ -123,4 +118,3 @@ def main():
 
 if __name__ == "__main__":
     main()
-