finemo.main

   1"""Main CLI module for the Fi-NeMo motif instance calling pipeline.
   2
   3This module provides the command-line interface for all Fi-NeMo operations:
   4- Data preprocessing from various genomic formats
   5- Motif hit calling using the Fi-NeMo algorithm
   6- Report generation and result visualization
   7- Post-processing operations (hit collapsing, intersection)
   8
   9The CLI supports multiple input formats including bigWig, HDF5 (ChromBPNet/BPNet),
  10and TF-MoDISco format.
  11"""
  12
  13from . import data_io
  14
  15import os
  16import argparse
  17import warnings
  18import inspect
  19from typing import Optional, List
  20
  21import polars as pl
  22
  23
  24def extract_regions_bw(
  25    peaks_path: str,
  26    chrom_order_path: Optional[str],
  27    fa_path: str,
  28    bw_paths: List[str],
  29    out_path: str,
  30    region_width: int = 1000,
  31) -> None:
  32    """Extract genomic regions and contribution scores from bigWig and FASTA files.
  33
  34    Parameters
  35    ----------
  36    peaks_path : str
  37        Path to ENCODE NarrowPeak format file.
  38    chrom_order_path : str, optional
  39        Path to chromosome ordering file.
  40    fa_path : str
  41        Path to genome FASTA file.
  42    bw_paths : List[str]
  43        List of bigWig file paths containing contribution scores.
  44    out_path : str
  45        Output path for NPZ file.
  46    region_width : int, default 1000
  47        Width of regions to extract around peak summits.
  48
  49    Notes
  50    -----
  51    BigWig files only provide projected contribution scores.
  52    """
  53    half_width = region_width // 2
  54
  55    # Load peak regions and extract sequences/contributions
  56    peaks_df = data_io.load_peaks(peaks_path, chrom_order_path, half_width)
  57    sequences, contribs = data_io.load_regions_from_bw(
  58        peaks_df, fa_path, bw_paths, half_width
  59    )
  60
  61    # Save processed data to NPZ format
  62    data_io.write_regions_npz(sequences, contribs, out_path, peaks_df=peaks_df)
  63
  64
  65def extract_regions_chrombpnet_h5(
  66    peaks_path: Optional[str],
  67    chrom_order_path: Optional[str],
  68    h5_paths: List[str],
  69    out_path: str,
  70    region_width: int = 1000,
  71) -> None:
  72    """Extract genomic regions and contribution scores from ChromBPNet HDF5 files.
  73
  74    Parameters
  75    ----------
  76    peaks_path : str, optional
  77        Path to ENCODE NarrowPeak format file. If None, lacks absolute coordinates.
  78    chrom_order_path : str, optional
  79        Path to chromosome ordering file.
  80    h5_paths : List[str]
  81        List of ChromBPNet HDF5 file paths.
  82    out_path : str
  83        Output path for NPZ file.
  84    region_width : int, default 1000
  85        Width of regions to extract around peak summits.
  86    """
  87    half_width = region_width // 2
  88
  89    if peaks_path is not None:
  90        peaks_df = data_io.load_peaks(peaks_path, chrom_order_path, half_width)
  91    else:
  92        peaks_df = None
  93
  94    sequences, contribs = data_io.load_regions_from_chrombpnet_h5(h5_paths, half_width)
  95
  96    data_io.write_regions_npz(sequences, contribs, out_path, peaks_df=peaks_df)
  97
  98
  99def extract_regions_bpnet_h5(
 100    peaks_path: Optional[str],
 101    chrom_order_path: Optional[str],
 102    h5_paths: List[str],
 103    out_path: str,
 104    region_width: int = 1000,
 105) -> None:
 106    """Extract genomic regions and contribution scores from BPNet HDF5 files.
 107
 108    Parameters
 109    ----------
 110    peaks_path : str, optional
 111        Path to ENCODE NarrowPeak format file. If None, output lacks absolute coordinates.
 112    chrom_order_path : str, optional
 113        Path to chromosome ordering file.
 114    h5_paths : List[str]
 115        List of BPNet HDF5 file paths.
 116    out_path : str
 117        Output path for NPZ file.
 118    region_width : int, default 1000
 119        Width of regions to extract around peak summits.
 120    """
 121    half_width = region_width // 2
 122
 123    if peaks_path is not None:
 124        peaks_df = data_io.load_peaks(peaks_path, chrom_order_path, half_width)
 125    else:
 126        peaks_df = None
 127
 128    sequences, contribs = data_io.load_regions_from_bpnet_h5(h5_paths, half_width)
 129
 130    data_io.write_regions_npz(sequences, contribs, out_path, peaks_df=peaks_df)
 131
 132
 133def extract_regions_modisco_fmt(
 134    peaks_path: Optional[str],
 135    chrom_order_path: Optional[str],
 136    shaps_paths: List[str],
 137    ohe_path: str,
 138    out_path: str,
 139    region_width: int = 1000,
 140) -> None:
 141    """Extract genomic regions and contribution scores from TF-MoDISco format files.
 142
 143    Parameters
 144    ----------
 145    peaks_path : str, optional
 146        Path to ENCODE NarrowPeak format file. If None, output lacks absolute coordinates.
 147    chrom_order_path : str, optional
 148        Path to chromosome ordering file.
 149    shaps_paths : List[str]
 150        List of paths to .npy/.npz files containing SHAP/attribution scores.
 151    ohe_path : str
 152        Path to .npy/.npz file containing one-hot encoded sequences.
 153    out_path : str
 154        Output path for NPZ file.
 155    region_width : int, default 1000
 156        Width of regions to extract around peak summits.
 157    """
 158    half_width = region_width // 2
 159
 160    if peaks_path is not None:
 161        peaks_df = data_io.load_peaks(peaks_path, chrom_order_path, half_width)
 162    else:
 163        peaks_df = None
 164
 165    sequences, contribs = data_io.load_regions_from_modisco_fmt(
 166        shaps_paths, ohe_path, half_width
 167    )
 168
 169    data_io.write_regions_npz(sequences, contribs, out_path, peaks_df=peaks_df)
 170
 171
 172def call_hits(
 173    regions_path: str,
 174    peaks_path: Optional[str],
 175    modisco_h5_path: str,
 176    chrom_order_path: Optional[str],
 177    motifs_include_path: Optional[str],
 178    motif_names_path: Optional[str],
 179    motif_lambdas_path: Optional[str],
 180    out_dir: str,
 181    cwm_trim_coords_path: Optional[str] = None,
 182    cwm_trim_thresholds_path: Optional[str] = None,
 183    cwm_trim_threshold_default: float = 0.3,
 184    lambda_default: float = 0.7,
 185    step_size_max: float = 3.0,
 186    step_size_min: float = 0.08,
 187    sqrt_transform: bool = False,
 188    convergence_tol: float = 0.0005,
 189    max_steps: int = 10000,
 190    batch_size: int = 2000,
 191    step_adjust: float = 0.7,
 192    device: Optional[str] = None,
 193    mode: str = "pp",
 194    no_post_filter: bool = False,
 195    compile_optimizer: bool = False,
 196) -> None:
 197    """Call motif hits using the Fi-NeMo algorithm on preprocessed genomic regions.
 198
 199    This function implements the core Fi-NeMo hit calling pipeline, which identifies
 200    motif instances by solving a sparse reconstruction problem using proximal gradient
 201    descent. The algorithm represents contribution scores as weighted combinations of
 202    motif CWMs at specific positions.
 203
 204    Parameters
 205    ----------
 206    regions_path : str
 207        Path to NPZ file containing preprocessed regions (sequences, contributions,
 208        and optional peak coordinates).
 209    peaks_path : str, optional
 210        DEPRECATED. Path to ENCODE NarrowPeak format file. Peak data should be
 211        included during preprocessing instead.
 212    modisco_h5_path : str
 213        Path to TF-MoDISco H5 file containing motif CWMs.
 214    chrom_order_path : str, optional
 215        DEPRECATED. Path to chromosome ordering file.
 216    motifs_include_path : str, optional
 217        Path to file listing motif names to include in analysis.
 218    motif_names_path : str, optional
 219        Path to file mapping motif IDs to custom names.
 220    motif_lambdas_path : str, optional
 221        Path to file specifying per-motif lambda values.
 222    out_dir : str
 223        Output directory for results.
 224    cwm_trim_coords_path : str, optional
 225        Path to file specifying custom motif trimming coordinates.
 226    cwm_trim_thresholds_path : str, optional
 227        Path to file specifying custom motif trimming thresholds.
 228    cwm_trim_threshold_default : float, default 0.3
 229        Default threshold for motif trimming.
 230    lambda_default : float, default 0.7
 231        Default L1 regularization weight.
 232    step_size_max : float, default 3.0
 233        Maximum optimization step size.
 234    step_size_min : float, default 0.08
 235        Minimum optimization step size.
 236    sqrt_transform : bool, default False
 237        Whether to apply signed square root transform to contributions.
 238    convergence_tol : float, default 0.0005
 239        Convergence tolerance for duality gap.
 240    max_steps : int, default 10000
 241        Maximum number of optimization steps.
 242    batch_size : int, default 2000
 243        Batch size for GPU processing.
 244    step_adjust : float, default 0.7
 245        Step size adjustment factor on divergence.
 246    device : str, optional
 247        DEPRECATED. Use CUDA_VISIBLE_DEVICES environment variable instead.
 248    mode : str, default "pp"
 249        Contribution type mode ('pp', 'ph', 'hp', 'hh') where 'p'=projected, 'h'=hypothetical.
 250    no_post_filter : bool, default False
 251        If True, skip post-hit-calling similarity filtering.
 252    compile_optimizer : bool, default False
 253        Whether to JIT-compile the optimizer for speed.
 254
 255    Notes
 256    -----
 257    The Fi-NeMo algorithm solves the optimization problem:
 258    minimize_c: ||contribs - reconstruction(c)||²₂ + λ||c||₁
 259    subject to: c ≥ 0
 260
 261    where c represents motif hit coefficients and reconstruction uses convolution
 262    with motif CWMs.
 263    """
 264
 265    params = locals()
 266    import torch
 267    from . import hitcaller
 268
 269    if device is not None:
 270        warnings.warn(
 271            "The `--device` flag is deprecated and will be removed in a future version. Please use the `CUDA_VISIBLE_DEVICES` environment variable to specify the GPU device."
 272        )
 273
 274    sequences, contribs, peaks_df, has_peaks = data_io.load_regions_npz(regions_path)
 275
 276    region_width = sequences.shape[2]
 277    if region_width % 2 != 0:
 278        raise ValueError(f"Region width of {region_width} is not divisible by 2.")
 279
 280    half_width = region_width // 2
 281    num_regions = contribs.shape[0]
 282
 283    if peaks_path is not None:
 284        warnings.warn(
 285            "Providing a peaks file to `call-hits` is deprecated, and this option will be removed in a future version. Peaks should instead be provided in the preprocessing step to be included in `regions.npz`."
 286        )
 287        peaks_df = data_io.load_peaks(peaks_path, chrom_order_path, half_width)
 288        has_peaks = True
 289
 290    if not has_peaks:
 291        warnings.warn(
 292            "No peak region data provided. Output hits will lack absolute genomic coordinates."
 293        )
 294
 295    if mode == "pp":
 296        motif_type = "cwm"
 297        use_hypothetical_contribs = False
 298    elif mode == "ph":
 299        motif_type = "cwm"
 300        use_hypothetical_contribs = True
 301    elif mode == "hp":
 302        motif_type = "hcwm"
 303        use_hypothetical_contribs = False
 304    elif mode == "hh":
 305        motif_type = "hcwm"
 306        use_hypothetical_contribs = True
 307    else:
 308        raise ValueError(
 309            f"Invalid mode: {mode}. Must be one of 'pp', 'ph', 'hp', 'hh'."
 310        )
 311
 312    if motifs_include_path is not None:
 313        motifs_include = data_io.load_txt(motifs_include_path)
 314    else:
 315        motifs_include = None
 316
 317    if motif_names_path is not None:
 318        motif_name_map = data_io.load_mapping(motif_names_path, str)
 319    else:
 320        motif_name_map = None
 321
 322    if motif_lambdas_path is not None:
 323        motif_lambdas = data_io.load_mapping(motif_lambdas_path, float)
 324    else:
 325        motif_lambdas = None
 326
 327    if cwm_trim_coords_path is not None:
 328        trim_coords = data_io.load_mapping_tuple(cwm_trim_coords_path, int)
 329    else:
 330        trim_coords = None
 331
 332    if cwm_trim_thresholds_path is not None:
 333        trim_thresholds = data_io.load_mapping(cwm_trim_thresholds_path, float)
 334    else:
 335        trim_thresholds = None
 336
 337    motifs_df, cwms, trim_masks, _ = data_io.load_modisco_motifs(
 338        modisco_h5_path,
 339        trim_coords,
 340        trim_thresholds,
 341        cwm_trim_threshold_default,
 342        motif_type,
 343        motifs_include,
 344        motif_name_map,
 345        motif_lambdas,
 346        lambda_default,
 347        True,
 348    )
 349    num_motifs = cwms.shape[0]
 350    motif_width = cwms.shape[2]
 351    lambdas = motifs_df.get_column("lambda").to_numpy(writable=True)
 352
 353    device_obj = torch.device(device) if device is not None else None
 354    hits_df, qc_df = hitcaller.fit_contribs(
 355        cwms,
 356        contribs,
 357        sequences,
 358        trim_masks,
 359        use_hypothetical_contribs,
 360        lambdas,
 361        step_size_max,
 362        step_size_min,
 363        sqrt_transform,
 364        convergence_tol,
 365        max_steps,
 366        batch_size,
 367        step_adjust,
 368        not no_post_filter,
 369        device_obj,
 370        compile_optimizer,
 371    )
 372
 373    os.makedirs(out_dir, exist_ok=True)
 374    out_path_qc = os.path.join(out_dir, "peaks_qc.tsv")
 375    data_io.write_hits(hits_df, peaks_df, motifs_df, qc_df, out_dir, motif_width)
 376    data_io.write_qc(qc_df, peaks_df, out_path_qc)
 377
 378    out_path_motif_df = os.path.join(out_dir, "motif_data.tsv")
 379    data_io.write_motifs_df(motifs_df, out_path_motif_df)
 380
 381    out_path_motif_cwms = os.path.join(out_dir, "motif_cwms.npy")
 382    data_io.write_motif_cwms(cwms, out_path_motif_cwms)
 383
 384    params |= {
 385        "region_width": region_width,
 386        "num_regions": num_regions,
 387        "untrimmed_motif_width": motif_width,
 388        "num_motifs": num_motifs,
 389    }
 390
 391    out_path_params = os.path.join(out_dir, "parameters.json")
 392    data_io.write_params(params, out_path_params)
 393
 394
 395def report(
 396    regions_path: str,
 397    hits_dir: str,
 398    modisco_h5_path: Optional[str],
 399    peaks_path: Optional[str],
 400    motifs_include_path: Optional[str],
 401    motif_names_path: Optional[str],
 402    out_dir: str,
 403    modisco_region_width: int = 400,
 404    cwm_trim_threshold: float = 0.3,
 405    compute_recall: bool = True,
 406    use_seqlets: bool = True,
 407) -> None:
 408    """Generate comprehensive HTML report with statistics and visualizations.
 409
 410    This function creates detailed analysis reports comparing Fi-NeMo hit calling
 411    results with TF-MoDISco seqlets, including performance metrics, distribution
 412    plots, and motif visualization. The report provides insights into hit calling
 413    quality and motif discovery accuracy.
 414
 415    Parameters
 416    ----------
 417    regions_path : str
 418        Path to NPZ file containing the same regions used for hit calling.
 419    hits_dir : str
 420        Path to directory containing Fi-NeMo hit calling outputs.
 421    modisco_h5_path : str, optional
 422        Path to TF-MoDISco H5 file. If None, seqlet comparisons are skipped.
 423    peaks_path : str, optional
 424        DEPRECATED. Peak coordinates should be included in regions file.
 425    motifs_include_path : str, optional
 426        DEPRECATED. This information is inferred from hit calling outputs.
 427    motif_names_path : str, optional
 428        DEPRECATED. This information is inferred from hit calling outputs.
 429    out_dir : str
 430        Output directory for report files.
 431    modisco_region_width : int, default 400
 432        Width of regions used by TF-MoDISco (needed for coordinate conversion).
 433    cwm_trim_threshold : float, default 0.3
 434        DEPRECATED. This information is inferred from hit calling outputs.
 435    compute_recall : bool, default True
 436        Whether to compute recall metrics against TF-MoDISco seqlets.
 437    use_seqlets : bool, default True
 438        Whether to include seqlet-based comparisons in the report.
 439
 440    Notes
 441    -----
 442    The generated report includes:
 443    - Hit vs seqlet count comparisons
 444    - Motif CWM visualizations
 445    - Hit statistic distributions
 446    - Co-occurrence heatmaps
 447    - Confusion matrices for overlapping motifs
 448    """
 449    from . import evaluation, visualization
 450
 451    sequences, contribs, peaks_df, _ = data_io.load_regions_npz(regions_path)
 452    if len(contribs.shape) == 3:
 453        regions = contribs * sequences
 454    elif len(contribs.shape) == 2:
 455        regions = contribs[:, None, :] * sequences
 456    else:
 457        raise ValueError(f"Unexpected contribs shape: {contribs.shape}")
 458
 459    half_width = regions.shape[2] // 2
 460    modisco_half_width = modisco_region_width // 2
 461
 462    if peaks_path is not None:
 463        warnings.warn(
 464            "Providing a peaks file to `report` is deprecated, and this option will be removed in a future version. Peaks should instead be provided in the preprocessing step to be included in `regions.npz`."
 465        )
 466        peaks_df = data_io.load_peaks(peaks_path, None, half_width)
 467
 468    if hits_dir.endswith(".tsv"):
 469        warnings.warn(
 470            "Passing a hits.tsv file to `finemo report` is deprecated. Please provide the directory containing the hits.tsv file instead."
 471        )
 472
 473        hits_path = hits_dir
 474
 475        hits_df = data_io.load_hits(hits_path, lazy=True)
 476
 477        if motifs_include_path is not None:
 478            motifs_include = data_io.load_txt(motifs_include_path)
 479        else:
 480            motifs_include = None
 481
 482        if motif_names_path is not None:
 483            motif_name_map = data_io.load_mapping(motif_names_path, str)
 484        else:
 485            motif_name_map = None
 486
 487        if modisco_h5_path is not None:
 488            motifs_df, cwms_modisco, _, motif_names = data_io.load_modisco_motifs(
 489                modisco_h5_path,
 490                None,
 491                None,
 492                cwm_trim_threshold,
 493                "cwm",
 494                motifs_include,
 495                motif_name_map,
 496                None,
 497                1.0,
 498                True,
 499            )
 500        else:
 501            # When no modisco_h5_path is provided in legacy TSV mode, we can't compute motifs
 502            # This will cause an error later, but that's expected behavior
 503            raise ValueError(
 504                "modisco_h5_path is required when providing a hits.tsv file directly"
 505            )
 506
 507    else:
 508        hits_df_path = os.path.join(hits_dir, "hits.tsv")
 509        hits_df = data_io.load_hits(hits_df_path, lazy=True)
 510
 511        motifs_df_path = os.path.join(hits_dir, "motif_data.tsv")
 512        motifs_df, motif_names = data_io.load_motifs_df(motifs_df_path)
 513
 514        cwms_modisco_path = os.path.join(hits_dir, "motif_cwms.npy")
 515        cwms_modisco = data_io.load_motif_cwms(cwms_modisco_path)
 516
 517        params_path = os.path.join(hits_dir, "parameters.json")
 518        params = data_io.load_params(params_path)
 519        cwm_trim_threshold = params["cwm_trim_threshold_default"]
 520
 521    if not use_seqlets:
 522        warnings.warn(
 523            "Usage of the `--no-seqlets` flag is deprecated and will be removed in a future version. Please omit the `--modisco-h5` argument instead."
 524        )
 525        seqlets_df = None
 526    elif modisco_h5_path is None:
 527        compute_recall = False
 528        seqlets_df = None
 529    else:
 530        seqlets_df = data_io.load_modisco_seqlets(
 531            modisco_h5_path,
 532            peaks_df,
 533            motifs_df,
 534            half_width,
 535            modisco_half_width,
 536            lazy=True,
 537        )
 538
 539    motif_width = cwms_modisco.shape[2]
 540
 541    # Convert to LazyFrame if needed and ensure motif_names is a list
 542    if isinstance(hits_df, pl.LazyFrame):
 543        hits_df_lazy: pl.LazyFrame = hits_df
 544    else:
 545        hits_df_lazy: pl.LazyFrame = hits_df.lazy()
 546
 547    motif_names_list: List[str] = list(motif_names)
 548
 549    occ_df, coooc = evaluation.get_motif_occurences(hits_df_lazy, motif_names_list)
 550
 551    report_data, report_df, motifs, trim_bounds = evaluation.tfmodisco_comparison(
 552        regions,
 553        sequences,
 554        hits_df,
 555        peaks_df,
 556        seqlets_df,
 557        motifs_df,
 558        cwms_modisco,
 559        motif_names_list,
 560        modisco_half_width,
 561        motif_width,
 562        compute_recall,
 563    )
 564
 565    if seqlets_df is not None:
 566        confusion_df, confusion_mat = evaluation.seqlet_confusion(
 567            hits_df, seqlets_df, peaks_df, motif_names_list, motif_width
 568        )
 569    else:
 570        confusion_df, confusion_mat = None, None
 571
 572    os.makedirs(out_dir, exist_ok=True)
 573
 574    occ_path = os.path.join(out_dir, "motif_occurrences.tsv")
 575    data_io.write_occ_df(occ_df, occ_path)
 576
 577    data_io.write_report_data(report_df, motifs, out_dir)
 578
 579    visualization.plot_hit_stat_distributions(hits_df_lazy, motif_names_list, out_dir)
 580    visualization.plot_hit_peak_distributions(occ_df, motif_names_list, out_dir)
 581    visualization.plot_peak_motif_indicator_heatmap(coooc, motif_names_list, out_dir)
 582
 583    plot_dir = os.path.join(out_dir, "motifs")
 584    visualization.plot_motifs(motifs, trim_bounds, plot_dir)
 585
 586    if seqlets_df is not None:
 587        seqlets_collected = (
 588            seqlets_df.collect() if isinstance(seqlets_df, pl.LazyFrame) else seqlets_df
 589        )
 590        seqlets_path = os.path.join(out_dir, "seqlets.tsv")
 591        data_io.write_modisco_seqlets(seqlets_collected, seqlets_path)
 592
 593        if confusion_df is not None and confusion_mat is not None:
 594            seqlet_confusion_path = os.path.join(out_dir, "seqlet_confusion.tsv")
 595            data_io.write_seqlet_confusion_df(confusion_df, seqlet_confusion_path)
 596
 597            visualization.plot_hit_vs_seqlet_counts(report_data, out_dir)
 598            visualization.plot_seqlet_confusion_heatmap(
 599                confusion_mat, motif_names_list, out_dir
 600            )
 601
 602    report_path = os.path.join(out_dir, "report.html")
 603    visualization.write_report(
 604        report_df, motif_names_list, report_path, compute_recall, seqlets_df is not None
 605    )
 606
 607
 608def collapse_hits(hits_path: str, out_path: str, overlap_frac: float = 0.2) -> None:
 609    """Collapse overlapping hits by selecting the best hit per overlapping group.
 610
 611    This function processes a set of motif hits and identifies overlapping hits,
 612    keeping only the hit with the highest similarity score within each overlapping
 613    group. This reduces redundancy in hit calls while preserving the most confident
 614    predictions.
 615
 616    Parameters
 617    ----------
 618    hits_path : str
 619        Path to input TSV file containing hit data (hits.tsv or hits_unique.tsv).
 620    out_path : str
 621        Path to output TSV file with additional 'is_primary' column.
 622    overlap_frac : float, default 0.2
 623        Minimum fractional overlap for considering hits as overlapping.
 624        For hits of lengths x and y, minimum overlap = overlap_frac * (x + y) / 2.
 625
 626    Notes
 627    -----
 628    The algorithm uses a sweep line approach with a heap data structure to
 629    efficiently identify overlapping intervals and select the best hit based
 630    on similarity scores.
 631    """
 632    from . import postprocessing
 633
 634    hits_df = data_io.load_hits(hits_path, lazy=False)
 635    hits_collapsed_df = postprocessing.collapse_hits(hits_df, overlap_frac)
 636
 637    data_io.write_hits_processed(
 638        hits_collapsed_df, out_path, schema=data_io.HITS_COLLAPSED_DTYPES
 639    )
 640
 641
 642def intersect_hits(hits_paths: List[str], out_path: str, relaxed: bool = False) -> None:
 643    """Find intersection of hits across multiple Fi-NeMo runs.
 644
 645    This function identifies motif hits that are consistently called across
 646    multiple independent runs, providing a way to assess reproducibility and
 647    identify high-confidence hits.
 648
 649    Parameters
 650    ----------
 651    hits_paths : List[str]
 652        List of paths to input TSV files from different runs.
 653    out_path : str
 654        Path to output TSV file containing intersection results.
 655        Duplicate columns are suffixed with run index.
 656    relaxed : bool, default False
 657        If True, uses relaxed intersection criteria based only on motif names
 658        and untrimmed coordinates. If False, assumes consistent region definitions
 659        and motif trimming across runs.
 660
 661    Notes
 662    -----
 663    The strict intersection mode requires consistent input regions and motif
 664    processing parameters across all runs. The relaxed mode is more permissive
 665    but may not be suitable when genomic coordinates are unavailable.
 666    """
 667    from . import postprocessing
 668
 669    hits_dfs = [data_io.load_hits(hits_path, lazy=False) for hits_path in hits_paths]
 670    hits_df = postprocessing.intersect_hits(hits_dfs, relaxed)
 671
 672    data_io.write_hits_processed(hits_df, out_path, schema=None)
 673
 674
 675def cli() -> None:
 676    """Command-line interface for the Fi-NeMo motif instance calling pipeline.
 677
 678    This function provides the main entry point for all Fi-NeMo operations including:
 679    - Data preprocessing from various formats (bigWig, HDF5, TF-MoDISco)
 680    - Motif hit calling using the Fi-NeMo algorithm
 681    - Report generation and visualization
 682    - Post-processing operations (hit collapsing, intersection)
 683    """
 684    parser = argparse.ArgumentParser()
 685    subparsers = parser.add_subparsers(required=True, dest="cmd")
 686
 687    extract_regions_bw_parser = subparsers.add_parser(
 688        "extract-regions-bw",
 689        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 690        help="Extract sequences and contributions from FASTA and bigwig files.",
 691    )
 692
 693    extract_regions_bw_parser.add_argument(
 694        "-p",
 695        "--peaks",
 696        type=str,
 697        required=True,
 698        help="A peak regions file in ENCODE NarrowPeak format.",
 699    )
 700    extract_regions_bw_parser.add_argument(
 701        "-C",
 702        "--chrom-order",
 703        type=str,
 704        default=None,
 705        help="A tab-delimited file with chromosome names in the first column to define sort order of chromosomes. Missing chromosomes are ordered as they appear in -p/--peaks.",
 706    )
 707    extract_regions_bw_parser.add_argument(
 708        "-f",
 709        "--fasta",
 710        type=str,
 711        required=True,
 712        help="A genome FASTA file. If an .fai index file doesn't exist in the same directory, it will be created.",
 713    )
 714    extract_regions_bw_parser.add_argument(
 715        "-b",
 716        "--bigwigs",
 717        type=str,
 718        required=True,
 719        nargs="+",
 720        help="One or more bigwig files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 721    )
 722
 723    extract_regions_bw_parser.add_argument(
 724        "-o",
 725        "--out-path",
 726        type=str,
 727        required=True,
 728        help="The path to the output .npz file.",
 729    )
 730
 731    extract_regions_bw_parser.add_argument(
 732        "-w",
 733        "--region-width",
 734        type=int,
 735        default=inspect.signature(extract_regions_bw)
 736        .parameters["region_width"]
 737        .default,
 738        help="The width of the input region centered around each peak summit.",
 739    )
 740
 741    extract_chrombpnet_regions_h5_parser = subparsers.add_parser(
 742        "extract-regions-chrombpnet-h5",
 743        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 744        help="Extract sequences and contributions from ChromBPNet contributions H5 files.",
 745    )
 746
 747    extract_chrombpnet_regions_h5_parser.add_argument(
 748        "-p",
 749        "--peaks",
 750        type=str,
 751        default=None,
 752        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 753    )
 754    extract_chrombpnet_regions_h5_parser.add_argument(
 755        "-C",
 756        "--chrom-order",
 757        type=str,
 758        default=None,
 759        help="A tab-delimited file with chromosome names in the first column to define sort order of chromosomes. Missing chromosomes are ordered as they appear in -p/--peaks.",
 760    )
 761
 762    extract_chrombpnet_regions_h5_parser.add_argument(
 763        "-c",
 764        "--h5s",
 765        type=str,
 766        required=True,
 767        nargs="+",
 768        help="One or more H5 files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 769    )
 770
 771    extract_chrombpnet_regions_h5_parser.add_argument(
 772        "-o",
 773        "--out-path",
 774        type=str,
 775        required=True,
 776        help="The path to the output .npz file.",
 777    )
 778
 779    extract_chrombpnet_regions_h5_parser.add_argument(
 780        "-w",
 781        "--region-width",
 782        type=int,
 783        default=inspect.signature(extract_regions_chrombpnet_h5)
 784        .parameters["region_width"]
 785        .default,
 786        help="The width of the input region centered around each peak summit.",
 787    )
 788
 789    extract_regions_h5_parser = subparsers.add_parser(
 790        "extract-regions-h5",
 791        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 792        help="Extract sequences and contributions from ChromBPNet contributions H5 files. DEPRECATED: Use `extract-regions-chrombpnet-h5` instead.",
 793    )
 794
 795    extract_regions_h5_parser.add_argument(
 796        "-p",
 797        "--peaks",
 798        type=str,
 799        default=None,
 800        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 801    )
 802    extract_regions_h5_parser.add_argument(
 803        "-C",
 804        "--chrom-order",
 805        type=str,
 806        default=None,
 807        help="A tab-delimited file with chromosome names in the first column to define sort order of chromosomes. Missing chromosomes are ordered as they appear in -p/--peaks.",
 808    )
 809
 810    extract_regions_h5_parser.add_argument(
 811        "-c",
 812        "--h5s",
 813        type=str,
 814        required=True,
 815        nargs="+",
 816        help="One or more H5 files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 817    )
 818
 819    extract_regions_h5_parser.add_argument(
 820        "-o",
 821        "--out-path",
 822        type=str,
 823        required=True,
 824        help="The path to the output .npz file.",
 825    )
 826
 827    extract_regions_h5_parser.add_argument(
 828        "-w",
 829        "--region-width",
 830        type=int,
 831        default=inspect.signature(extract_regions_chrombpnet_h5)
 832        .parameters["region_width"]
 833        .default,
 834        help="The width of the input region centered around each peak summit.",
 835    )
 836
 837    extract_bpnet_regions_h5_parser = subparsers.add_parser(
 838        "extract-regions-bpnet-h5",
 839        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 840        help="Extract sequences and contributions from BPNet contributions H5 files.",
 841    )
 842
 843    extract_bpnet_regions_h5_parser.add_argument(
 844        "-p",
 845        "--peaks",
 846        type=str,
 847        default=None,
 848        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 849    )
 850    extract_bpnet_regions_h5_parser.add_argument(
 851        "-C",
 852        "--chrom-order",
 853        type=str,
 854        default=None,
 855        help="A tab-delimited file with chromosome names in the first column to define sort order of chromosomes. Missing chromosomes are ordered as they appear in -p/--peaks.",
 856    )
 857
 858    extract_bpnet_regions_h5_parser.add_argument(
 859        "-c",
 860        "--h5s",
 861        type=str,
 862        required=True,
 863        nargs="+",
 864        help="One or more H5 files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 865    )
 866
 867    extract_bpnet_regions_h5_parser.add_argument(
 868        "-o",
 869        "--out-path",
 870        type=str,
 871        required=True,
 872        help="The path to the output .npz file.",
 873    )
 874
 875    extract_bpnet_regions_h5_parser.add_argument(
 876        "-w",
 877        "--region-width",
 878        type=int,
 879        default=inspect.signature(extract_regions_bpnet_h5)
 880        .parameters["region_width"]
 881        .default,
 882        help="The width of the input region centered around each peak summit.",
 883    )
 884
 885    extract_regions_modisco_fmt_parser = subparsers.add_parser(
 886        "extract-regions-modisco-fmt",
 887        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 888        help="Extract sequences and contributions from tfmodisco-lite input files.",
 889    )
 890
 891    extract_regions_modisco_fmt_parser.add_argument(
 892        "-p",
 893        "--peaks",
 894        type=str,
 895        default=None,
 896        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 897    )
 898    extract_regions_modisco_fmt_parser.add_argument(
 899        "-C",
 900        "--chrom-order",
 901        type=str,
 902        default=None,
 903        help="A tab-delimited file with chromosome names in the first column to define sort order of chromosomes. Missing chromosomes are ordered as they appear in -p/--peaks.",
 904    )
 905
 906    extract_regions_modisco_fmt_parser.add_argument(
 907        "-s",
 908        "--sequences",
 909        type=str,
 910        required=True,
 911        help="A .npy or .npz file containing one-hot encoded sequences.",
 912    )
 913
 914    extract_regions_modisco_fmt_parser.add_argument(
 915        "-a",
 916        "--attributions",
 917        type=str,
 918        required=True,
 919        nargs="+",
 920        help="One or more .npy or .npz files of hypothetical contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 921    )
 922
 923    extract_regions_modisco_fmt_parser.add_argument(
 924        "-o",
 925        "--out-path",
 926        type=str,
 927        required=True,
 928        help="The path to the output .npz file.",
 929    )
 930
 931    extract_regions_modisco_fmt_parser.add_argument(
 932        "-w",
 933        "--region-width",
 934        type=int,
 935        default=inspect.signature(extract_regions_modisco_fmt)
 936        .parameters["region_width"]
 937        .default,
 938        help="The width of the input region centered around each peak summit.",
 939    )
 940
 941    call_hits_parser = subparsers.add_parser(
 942        "call-hits",
 943        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 944        help="Call hits on provided sequences, contributions, and motif CWM's.",
 945    )
 946
 947    call_hits_parser.add_argument(
 948        "-M",
 949        "--mode",
 950        type=str,
 951        default=inspect.signature(call_hits).parameters["mode"].default,
 952        choices={"pp", "ph", "hp", "hh"},
 953        help="The type of attributions to use for CWM's and input contribution scores, respectively. 'h' for hypothetical and 'p' for projected.",
 954    )
 955
 956    call_hits_parser.add_argument(
 957        "-r",
 958        "--regions",
 959        type=str,
 960        required=True,
 961        help="A .npz file of input sequences, contributions, and coordinates. Can be generated using `finemo extract-regions-*` subcommands.",
 962    )
 963    call_hits_parser.add_argument(
 964        "-m",
 965        "--modisco-h5",
 966        type=str,
 967        required=True,
 968        help="A tfmodisco-lite output H5 file of motif patterns.",
 969    )
 970
 971    call_hits_parser.add_argument(
 972        "-p",
 973        "--peaks",
 974        type=str,
 975        default=None,
 976        help="DEPRECATED: Please provide this file to a preprocessing `finemo extract-regions-*` subcommand instead.",
 977    )
 978    call_hits_parser.add_argument(
 979        "-C",
 980        "--chrom-order",
 981        type=str,
 982        default=None,
 983        help="DEPRECATED: Please provide this file to a preprocessing `finemo extract-regions-*` subcommand instead.",
 984    )
 985
 986    call_hits_parser.add_argument(
 987        "-I",
 988        "--motifs-include",
 989        type=str,
 990        default=None,
 991        help="A tab-delimited file with tfmodisco motif names (e.g pos_patterns.pattern_0) in the first column to include in hit calling. If omitted, all motifs in the modisco H5 file are used.",
 992    )
 993    call_hits_parser.add_argument(
 994        "-N",
 995        "--motif-names",
 996        type=str,
 997        default=None,
 998        help="A tab-delimited file with tfmodisco motif names (e.g pos_patterns.pattern_0) in the first column and custom names in the second column. Omitted motifs default to tfmodisco names.",
 999    )
1000
1001    call_hits_parser.add_argument(
1002        "-o",
1003        "--out-dir",
1004        type=str,
1005        required=True,
1006        help="The path to the output directory.",
1007    )
1008
1009    call_hits_parser.add_argument(
1010        "-t",
1011        "--cwm-trim-threshold",
1012        type=float,
1013        default=inspect.signature(call_hits)
1014        .parameters["cwm_trim_threshold_default"]
1015        .default,
1016        help="The default threshold to determine motif start and end positions within the full CWMs.",
1017    )
1018    call_hits_parser.add_argument(
1019        "-T",
1020        "--cwm-trim-thresholds",
1021        type=str,
1022        default=None,
1023        help="A tab-delimited file with tfmodisco motif names (e.g pos_patterns.pattern_0) in the first column and custom trim thresholds in the second column. Omitted motifs default to the `--cwm-trim-threshold` value.",
1024    )
1025    call_hits_parser.add_argument(
1026        "-R",
1027        "--cwm-trim-coords",
1028        type=str,
1029        default=None,
1030        help="A tab-delimited file with tfmodisco motif names (e.g pos_patterns.pattern_0) in the first column and custom trim start and end coordinates in the second and third columns, respectively. Omitted motifs default to `--cwm-trim-thresholds` values.",
1031    )
1032
1033    call_hits_parser.add_argument(
1034        "-l",
1035        "--global-lambda",
1036        type=float,
1037        default=inspect.signature(call_hits).parameters["lambda_default"].default,
1038        help="The default L1 regularization weight determining the sparsity of hits.",
1039    )
1040    call_hits_parser.add_argument(
1041        "-L",
1042        "--motif-lambdas",
1043        type=str,
1044        default=None,
1045        help="A tab-delimited file with tfmodisco motif names (e.g pos_patterns.pattern_0) in the first column and motif-specific lambdas in the second column. Omitted motifs default to the `--global-lambda` value.",
1046    )
1047    call_hits_parser.add_argument(
1048        "-a",
1049        "--alpha",
1050        type=float,
1051        default=None,
1052        help="DEPRECATED: Please use the `--lambda` argument instead.",
1053    )
1054    call_hits_parser.add_argument(
1055        "-A",
1056        "--motif-alphas",
1057        type=str,
1058        default=None,
1059        help="DEPRECATED: Please use the `--motif-lambdas` argument instead.",
1060    )
1061
1062    call_hits_parser.add_argument(
1063        "-f",
1064        "--no-post-filter",
1065        action="store_true",
1066        help="Do not perform post-hit-calling filtering. By default, hits are filtered based on a minimum cosine similarity of `lambda` with the input contributions.",
1067    )
1068    call_hits_parser.add_argument(
1069        "-q",
1070        "--sqrt-transform",
1071        action="store_true",
1072        help="Apply a signed square root transform to the input contributions and CWMs before hit calling.",
1073    )
1074    call_hits_parser.add_argument(
1075        "-s",
1076        "--step-size-max",
1077        type=float,
1078        default=inspect.signature(call_hits).parameters["step_size_max"].default,
1079        help="The maximum optimizer step size.",
1080    )
1081    call_hits_parser.add_argument(
1082        "-i",
1083        "--step-size-min",
1084        type=float,
1085        default=inspect.signature(call_hits).parameters["step_size_min"].default,
1086        help="The minimum optimizer step size.",
1087    )
1088    call_hits_parser.add_argument(
1089        "-j",
1090        "--step-adjust",
1091        type=float,
1092        default=inspect.signature(call_hits).parameters["step_adjust"].default,
1093        help="The optimizer step size adjustment factor. If the optimizer diverges, the step size is multiplicatively adjusted by this factor",
1094    )
1095    call_hits_parser.add_argument(
1096        "-c",
1097        "--convergence-tol",
1098        type=float,
1099        default=inspect.signature(call_hits).parameters["convergence_tol"].default,
1100        help="The tolerance for determining convergence. The optimizer exits when the duality gap is less than the tolerance.",
1101    )
1102    call_hits_parser.add_argument(
1103        "-S",
1104        "--max-steps",
1105        type=int,
1106        default=inspect.signature(call_hits).parameters["max_steps"].default,
1107        help="The maximum number of optimization steps.",
1108    )
1109    call_hits_parser.add_argument(
1110        "-b",
1111        "--batch-size",
1112        type=int,
1113        default=inspect.signature(call_hits).parameters["batch_size"].default,
1114        help="The batch size used for optimization.",
1115    )
1116    call_hits_parser.add_argument(
1117        "-d",
1118        "--device",
1119        type=str,
1120        default=None,
1121        help="DEPRECATED: Please use the `CUDA_VISIBLE_DEVICES` environment variable to specify the GPU device.",
1122    )
1123    call_hits_parser.add_argument(
1124        "-J",
1125        "--compile",
1126        action="store_true",
1127        help="JIT-compile the optimizer for faster performance. This may not be supported on older GPUs.",
1128    )
1129
1130    report_parser = subparsers.add_parser(
1131        "report",
1132        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
1133        help="Generate statistics and visualizations from hits and tfmodisco-lite motif data.",
1134    )
1135
1136    report_parser.add_argument(
1137        "-r",
1138        "--regions",
1139        type=str,
1140        required=True,
1141        help="A .npz file containing input sequences, contributions, and coordinates. Must be the same as that used for `finemo call-hits`.",
1142    )
1143    report_parser.add_argument(
1144        "-H",
1145        "--hits",
1146        type=str,
1147        required=True,
1148        help="The output directory generated by the `finemo call-hits` command on the regions specified in `--regions`.",
1149    )
1150    report_parser.add_argument(
1151        "-p",
1152        "--peaks",
1153        type=str,
1154        default=None,
1155        help="DEPRECATED: Please provide this file to a preprocessing `finemo extract-regions-*` subcommand instead.",
1156    )
1157    report_parser.add_argument(
1158        "-m",
1159        "--modisco-h5",
1160        type=str,
1161        default=None,
1162        help="The tfmodisco-lite output H5 file of motif patterns. Must be the same as that used for hit calling unless `--no-recall` is set. If omitted, seqlet-derived metrics will not be computed.",
1163    )
1164    report_parser.add_argument(
1165        "-I",
1166        "--motifs-include",
1167        type=str,
1168        default=None,
1169        help="DEPRECATED: This information is now inferred from the outputs of `finemo call-hits`.",
1170    )
1171    report_parser.add_argument(
1172        "-N",
1173        "--motif-names",
1174        type=str,
1175        default=None,
1176        help="DEPRECATED: This information is now inferred from the outputs of `finemo call-hits`.",
1177    )
1178
1179    report_parser.add_argument(
1180        "-o",
1181        "--out-dir",
1182        type=str,
1183        required=True,
1184        help="The path to the report output directory.",
1185    )
1186
1187    report_parser.add_argument(
1188        "-W",
1189        "--modisco-region-width",
1190        type=int,
1191        default=inspect.signature(report).parameters["modisco_region_width"].default,
1192        help="The width of the region around each peak summit used by tfmodisco-lite.",
1193    )
1194    report_parser.add_argument(
1195        "-t",
1196        "--cwm-trim-threshold",
1197        type=float,
1198        default=inspect.signature(report).parameters["cwm_trim_threshold"].default,
1199        help="DEPRECATED: This information is now inferred from the outputs of `finemo call-hits`.",
1200    )
1201    report_parser.add_argument(
1202        "-n",
1203        "--no-recall",
1204        action="store_true",
1205        help="Do not compute motif recall metrics.",
1206    )
1207    report_parser.add_argument(
1208        "-s",
1209        "--no-seqlets",
1210        action="store_true",
1211        help="DEPRECATED: Please omit the `--modisco-h5` argument instead.",
1212    )
1213
1214    collapse_hits_parser = subparsers.add_parser(
1215        "collapse-hits",
1216        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
1217        help="Identify best hit by motif similarity among sets of overlapping hits.",
1218    )
1219
1220    collapse_hits_parser.add_argument(
1221        "-i",
1222        "--hits",
1223        type=str,
1224        required=True,
1225        help="The `hits.tsv` or `hits_unique.tsv` file from `call-hits`.",
1226    )
1227    collapse_hits_parser.add_argument(
1228        "-o",
1229        "--out-path",
1230        type=str,
1231        required=True,
1232        help='The path to the output .tsv file with an additional "is_primary" column.',
1233    )
1234    collapse_hits_parser.add_argument(
1235        "-O",
1236        "--overlap-frac",
1237        type=float,
1238        default=inspect.signature(collapse_hits).parameters["overlap_frac"].default,
1239        help="The threshold for determining overlapping hits. For two hits with lengths x and y, the minimum overlap is defined as `overlap_frac * (x + y) / 2`. The default value of 0.2 means that two hits must overlap by at least 20% of their average lengths to be considered overlapping.",
1240    )
1241
1242    intersect_hits_parser = subparsers.add_parser(
1243        "intersect-hits",
1244        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
1245        help="Intersect hits across multiple runs.",
1246    )
1247
1248    intersect_hits_parser.add_argument(
1249        "-i",
1250        "--hits",
1251        type=str,
1252        required=True,
1253        nargs="+",
1254        help="One or more hits.tsv or hits_unique.tsv files, with paths delimited by whitespace.",
1255    )
1256    intersect_hits_parser.add_argument(
1257        "-o",
1258        "--out-path",
1259        type=str,
1260        required=True,
1261        help="The path to the output .tsv file. Duplicate columns are suffixed with the positional index of the input file.",
1262    )
1263    intersect_hits_parser.add_argument(
1264        "-r",
1265        "--relaxed",
1266        action="store_true",
1267        help="Use relaxed intersection criteria, using only motif names and untrimmed coordinates. By default, the intersection assumes consistent region definitions and motif trimming. This option is not recommended if genomic coordinates are unavailable.",
1268    )
1269
1270    args = parser.parse_args()
1271
1272    if args.cmd == "extract-regions-bw":
1273        extract_regions_bw(
1274            args.peaks,
1275            args.chrom_order,
1276            args.fasta,
1277            args.bigwigs,
1278            args.out_path,
1279            args.region_width,
1280        )
1281
1282    elif args.cmd == "extract-regions-chrombpnet-h5":
1283        extract_regions_chrombpnet_h5(
1284            args.peaks, args.chrom_order, args.h5s, args.out_path, args.region_width
1285        )
1286
1287    elif args.cmd == "extract-regions-h5":
1288        print(
1289            "WARNING: The `extract-regions-h5` command is deprecated. Use `extract-regions-chrombpnet-h5` instead."
1290        )
1291        extract_regions_chrombpnet_h5(
1292            args.peaks, args.chrom_order, args.h5s, args.out_path, args.region_width
1293        )
1294
1295    elif args.cmd == "extract-regions-bpnet-h5":
1296        extract_regions_bpnet_h5(
1297            args.peaks, args.chrom_order, args.h5s, args.out_path, args.region_width
1298        )
1299
1300    elif args.cmd == "extract-regions-modisco-fmt":
1301        extract_regions_modisco_fmt(
1302            args.peaks,
1303            args.chrom_order,
1304            args.attributions,
1305            args.sequences,
1306            args.out_path,
1307            args.region_width,
1308        )
1309
1310    elif args.cmd == "call-hits":
1311        if args.alpha is not None:
1312            warnings.warn(
1313                "The `--alpha` flag is deprecated and will be removed in a future version. Please use the `--global-lambda` flag instead."
1314            )
1315            args.global_lambda = args.alpha
1316        if args.motif_alphas is not None:
1317            warnings.warn(
1318                "The `--motif-alphas` flag is deprecated and will be removed in a future version. Please use the `--motif-lambdas` flag instead."
1319            )
1320            args.motif_lambdas = args.motif_alphas
1321
1322        call_hits(
1323            args.regions,
1324            args.peaks,
1325            args.modisco_h5,
1326            args.chrom_order,
1327            args.motifs_include,
1328            args.motif_names,
1329            args.motif_lambdas,
1330            args.out_dir,
1331            args.cwm_trim_coords,
1332            args.cwm_trim_thresholds,
1333            args.cwm_trim_threshold,
1334            args.global_lambda,
1335            args.step_size_max,
1336            args.step_size_min,
1337            args.sqrt_transform,
1338            args.convergence_tol,
1339            args.max_steps,
1340            args.batch_size,
1341            args.step_adjust,
1342            args.device,
1343            args.mode,
1344            args.no_post_filter,
1345            args.compile,
1346        )
1347
1348    elif args.cmd == "report":
1349        report(
1350            args.regions,
1351            args.hits,
1352            args.modisco_h5,
1353            args.peaks,
1354            args.motifs_include,
1355            args.motif_names,
1356            args.out_dir,
1357            args.modisco_region_width,
1358            args.cwm_trim_threshold,
1359            not args.no_recall,
1360            not args.no_seqlets,
1361        )
1362
1363    elif args.cmd == "collapse-hits":
1364        collapse_hits(args.hits, args.out_path, args.overlap_frac)
1365
1366    elif args.cmd == "intersect-hits":
1367        intersect_hits(args.hits, args.out_path, args.relaxed)