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, cwms, trim_bounds = evaluation.tfmodisco_comparison(
 552        regions,
 553        hits_df,
 554        peaks_df,
 555        seqlets_df,
 556        motifs_df,
 557        cwms_modisco,
 558        motif_names_list,
 559        modisco_half_width,
 560        motif_width,
 561        compute_recall,
 562    )
 563
 564    if seqlets_df is not None:
 565        confusion_df, confusion_mat = evaluation.seqlet_confusion(
 566            hits_df, seqlets_df, peaks_df, motif_names_list, motif_width
 567        )
 568    else:
 569        confusion_df, confusion_mat = None, None
 570
 571    os.makedirs(out_dir, exist_ok=True)
 572
 573    occ_path = os.path.join(out_dir, "motif_occurrences.tsv")
 574    data_io.write_occ_df(occ_df, occ_path)
 575
 576    data_io.write_report_data(report_df, cwms, out_dir)
 577
 578    visualization.plot_hit_stat_distributions(hits_df_lazy, motif_names_list, out_dir)
 579    visualization.plot_hit_peak_distributions(occ_df, motif_names_list, out_dir)
 580    visualization.plot_peak_motif_indicator_heatmap(coooc, motif_names_list, out_dir)
 581
 582    plot_dir = os.path.join(out_dir, "CWMs")
 583    visualization.plot_cwms(cwms, trim_bounds, plot_dir)
 584
 585    if seqlets_df is not None:
 586        seqlets_collected = (
 587            seqlets_df.collect() if isinstance(seqlets_df, pl.LazyFrame) else seqlets_df
 588        )
 589        seqlets_path = os.path.join(out_dir, "seqlets.tsv")
 590        data_io.write_modisco_seqlets(seqlets_collected, seqlets_path)
 591
 592        if confusion_df is not None and confusion_mat is not None:
 593            seqlet_confusion_path = os.path.join(out_dir, "seqlet_confusion.tsv")
 594            data_io.write_seqlet_confusion_df(confusion_df, seqlet_confusion_path)
 595
 596            visualization.plot_hit_vs_seqlet_counts(report_data, out_dir)
 597            visualization.plot_seqlet_confusion_heatmap(
 598                confusion_mat, motif_names_list, out_dir
 599            )
 600
 601    report_path = os.path.join(out_dir, "report.html")
 602    visualization.write_report(
 603        report_df, motif_names_list, report_path, compute_recall, seqlets_df is not None
 604    )
 605
 606
 607def collapse_hits(hits_path: str, out_path: str, overlap_frac: float = 0.2) -> None:
 608    """Collapse overlapping hits by selecting the best hit per overlapping group.
 609
 610    This function processes a set of motif hits and identifies overlapping hits,
 611    keeping only the hit with the highest similarity score within each overlapping
 612    group. This reduces redundancy in hit calls while preserving the most confident
 613    predictions.
 614
 615    Parameters
 616    ----------
 617    hits_path : str
 618        Path to input TSV file containing hit data (hits.tsv or hits_unique.tsv).
 619    out_path : str
 620        Path to output TSV file with additional 'is_primary' column.
 621    overlap_frac : float, default 0.2
 622        Minimum fractional overlap for considering hits as overlapping.
 623        For hits of lengths x and y, minimum overlap = overlap_frac * (x + y) / 2.
 624
 625    Notes
 626    -----
 627    The algorithm uses a sweep line approach with a heap data structure to
 628    efficiently identify overlapping intervals and select the best hit based
 629    on similarity scores.
 630    """
 631    from . import postprocessing
 632
 633    hits_df = data_io.load_hits(hits_path, lazy=False)
 634    hits_collapsed_df = postprocessing.collapse_hits(hits_df, overlap_frac)
 635
 636    data_io.write_hits_processed(
 637        hits_collapsed_df, out_path, schema=data_io.HITS_COLLAPSED_DTYPES
 638    )
 639
 640
 641def intersect_hits(hits_paths: List[str], out_path: str, relaxed: bool = False) -> None:
 642    """Find intersection of hits across multiple Fi-NeMo runs.
 643
 644    This function identifies motif hits that are consistently called across
 645    multiple independent runs, providing a way to assess reproducibility and
 646    identify high-confidence hits.
 647
 648    Parameters
 649    ----------
 650    hits_paths : List[str]
 651        List of paths to input TSV files from different runs.
 652    out_path : str
 653        Path to output TSV file containing intersection results.
 654        Duplicate columns are suffixed with run index.
 655    relaxed : bool, default False
 656        If True, uses relaxed intersection criteria based only on motif names
 657        and untrimmed coordinates. If False, assumes consistent region definitions
 658        and motif trimming across runs.
 659
 660    Notes
 661    -----
 662    The strict intersection mode requires consistent input regions and motif
 663    processing parameters across all runs. The relaxed mode is more permissive
 664    but may not be suitable when genomic coordinates are unavailable.
 665    """
 666    from . import postprocessing
 667
 668    hits_dfs = [data_io.load_hits(hits_path, lazy=False) for hits_path in hits_paths]
 669    hits_df = postprocessing.intersect_hits(hits_dfs, relaxed)
 670
 671    data_io.write_hits_processed(hits_df, out_path, schema=None)
 672
 673
 674def cli() -> None:
 675    """Command-line interface for the Fi-NeMo motif instance calling pipeline.
 676
 677    This function provides the main entry point for all Fi-NeMo operations including:
 678    - Data preprocessing from various formats (bigWig, HDF5, TF-MoDISco)
 679    - Motif hit calling using the Fi-NeMo algorithm
 680    - Report generation and visualization
 681    - Post-processing operations (hit collapsing, intersection)
 682    """
 683    parser = argparse.ArgumentParser()
 684    subparsers = parser.add_subparsers(required=True, dest="cmd")
 685
 686    extract_regions_bw_parser = subparsers.add_parser(
 687        "extract-regions-bw",
 688        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 689        help="Extract sequences and contributions from FASTA and bigwig files.",
 690    )
 691
 692    extract_regions_bw_parser.add_argument(
 693        "-p",
 694        "--peaks",
 695        type=str,
 696        required=True,
 697        help="A peak regions file in ENCODE NarrowPeak format.",
 698    )
 699    extract_regions_bw_parser.add_argument(
 700        "-C",
 701        "--chrom-order",
 702        type=str,
 703        default=None,
 704        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.",
 705    )
 706    extract_regions_bw_parser.add_argument(
 707        "-f",
 708        "--fasta",
 709        type=str,
 710        required=True,
 711        help="A genome FASTA file. If an .fai index file doesn't exist in the same directory, it will be created.",
 712    )
 713    extract_regions_bw_parser.add_argument(
 714        "-b",
 715        "--bigwigs",
 716        type=str,
 717        required=True,
 718        nargs="+",
 719        help="One or more bigwig files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 720    )
 721
 722    extract_regions_bw_parser.add_argument(
 723        "-o",
 724        "--out-path",
 725        type=str,
 726        required=True,
 727        help="The path to the output .npz file.",
 728    )
 729
 730    extract_regions_bw_parser.add_argument(
 731        "-w",
 732        "--region-width",
 733        type=int,
 734        default=inspect.signature(extract_regions_bw)
 735        .parameters["region_width"]
 736        .default,
 737        help="The width of the input region centered around each peak summit.",
 738    )
 739
 740    extract_chrombpnet_regions_h5_parser = subparsers.add_parser(
 741        "extract-regions-chrombpnet-h5",
 742        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 743        help="Extract sequences and contributions from ChromBPNet contributions H5 files.",
 744    )
 745
 746    extract_chrombpnet_regions_h5_parser.add_argument(
 747        "-p",
 748        "--peaks",
 749        type=str,
 750        default=None,
 751        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 752    )
 753    extract_chrombpnet_regions_h5_parser.add_argument(
 754        "-C",
 755        "--chrom-order",
 756        type=str,
 757        default=None,
 758        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.",
 759    )
 760
 761    extract_chrombpnet_regions_h5_parser.add_argument(
 762        "-c",
 763        "--h5s",
 764        type=str,
 765        required=True,
 766        nargs="+",
 767        help="One or more H5 files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 768    )
 769
 770    extract_chrombpnet_regions_h5_parser.add_argument(
 771        "-o",
 772        "--out-path",
 773        type=str,
 774        required=True,
 775        help="The path to the output .npz file.",
 776    )
 777
 778    extract_chrombpnet_regions_h5_parser.add_argument(
 779        "-w",
 780        "--region-width",
 781        type=int,
 782        default=inspect.signature(extract_regions_chrombpnet_h5)
 783        .parameters["region_width"]
 784        .default,
 785        help="The width of the input region centered around each peak summit.",
 786    )
 787
 788    extract_regions_h5_parser = subparsers.add_parser(
 789        "extract-regions-h5",
 790        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 791        help="Extract sequences and contributions from ChromBPNet contributions H5 files. DEPRECATED: Use `extract-regions-chrombpnet-h5` instead.",
 792    )
 793
 794    extract_regions_h5_parser.add_argument(
 795        "-p",
 796        "--peaks",
 797        type=str,
 798        default=None,
 799        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 800    )
 801    extract_regions_h5_parser.add_argument(
 802        "-C",
 803        "--chrom-order",
 804        type=str,
 805        default=None,
 806        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.",
 807    )
 808
 809    extract_regions_h5_parser.add_argument(
 810        "-c",
 811        "--h5s",
 812        type=str,
 813        required=True,
 814        nargs="+",
 815        help="One or more H5 files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 816    )
 817
 818    extract_regions_h5_parser.add_argument(
 819        "-o",
 820        "--out-path",
 821        type=str,
 822        required=True,
 823        help="The path to the output .npz file.",
 824    )
 825
 826    extract_regions_h5_parser.add_argument(
 827        "-w",
 828        "--region-width",
 829        type=int,
 830        default=inspect.signature(extract_regions_chrombpnet_h5)
 831        .parameters["region_width"]
 832        .default,
 833        help="The width of the input region centered around each peak summit.",
 834    )
 835
 836    extract_bpnet_regions_h5_parser = subparsers.add_parser(
 837        "extract-regions-bpnet-h5",
 838        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 839        help="Extract sequences and contributions from BPNet contributions H5 files.",
 840    )
 841
 842    extract_bpnet_regions_h5_parser.add_argument(
 843        "-p",
 844        "--peaks",
 845        type=str,
 846        default=None,
 847        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 848    )
 849    extract_bpnet_regions_h5_parser.add_argument(
 850        "-C",
 851        "--chrom-order",
 852        type=str,
 853        default=None,
 854        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.",
 855    )
 856
 857    extract_bpnet_regions_h5_parser.add_argument(
 858        "-c",
 859        "--h5s",
 860        type=str,
 861        required=True,
 862        nargs="+",
 863        help="One or more H5 files of contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 864    )
 865
 866    extract_bpnet_regions_h5_parser.add_argument(
 867        "-o",
 868        "--out-path",
 869        type=str,
 870        required=True,
 871        help="The path to the output .npz file.",
 872    )
 873
 874    extract_bpnet_regions_h5_parser.add_argument(
 875        "-w",
 876        "--region-width",
 877        type=int,
 878        default=inspect.signature(extract_regions_bpnet_h5)
 879        .parameters["region_width"]
 880        .default,
 881        help="The width of the input region centered around each peak summit.",
 882    )
 883
 884    extract_regions_modisco_fmt_parser = subparsers.add_parser(
 885        "extract-regions-modisco-fmt",
 886        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 887        help="Extract sequences and contributions from tfmodisco-lite input files.",
 888    )
 889
 890    extract_regions_modisco_fmt_parser.add_argument(
 891        "-p",
 892        "--peaks",
 893        type=str,
 894        default=None,
 895        help="A peak regions file in ENCODE NarrowPeak format. If omitted, downstream outputs will lack absolute genomic coordinates.",
 896    )
 897    extract_regions_modisco_fmt_parser.add_argument(
 898        "-C",
 899        "--chrom-order",
 900        type=str,
 901        default=None,
 902        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.",
 903    )
 904
 905    extract_regions_modisco_fmt_parser.add_argument(
 906        "-s",
 907        "--sequences",
 908        type=str,
 909        required=True,
 910        help="A .npy or .npz file containing one-hot encoded sequences.",
 911    )
 912
 913    extract_regions_modisco_fmt_parser.add_argument(
 914        "-a",
 915        "--attributions",
 916        type=str,
 917        required=True,
 918        nargs="+",
 919        help="One or more .npy or .npz files of hypothetical contribution scores, with paths delimited by whitespace. Scores are averaged across files.",
 920    )
 921
 922    extract_regions_modisco_fmt_parser.add_argument(
 923        "-o",
 924        "--out-path",
 925        type=str,
 926        required=True,
 927        help="The path to the output .npz file.",
 928    )
 929
 930    extract_regions_modisco_fmt_parser.add_argument(
 931        "-w",
 932        "--region-width",
 933        type=int,
 934        default=inspect.signature(extract_regions_modisco_fmt)
 935        .parameters["region_width"]
 936        .default,
 937        help="The width of the input region centered around each peak summit.",
 938    )
 939
 940    call_hits_parser = subparsers.add_parser(
 941        "call-hits",
 942        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
 943        help="Call hits on provided sequences, contributions, and motif CWM's.",
 944    )
 945
 946    call_hits_parser.add_argument(
 947        "-M",
 948        "--mode",
 949        type=str,
 950        default=inspect.signature(call_hits).parameters["mode"].default,
 951        choices={"pp", "ph", "hp", "hh"},
 952        help="The type of attributions to use for CWM's and input contribution scores, respectively. 'h' for hypothetical and 'p' for projected.",
 953    )
 954
 955    call_hits_parser.add_argument(
 956        "-r",
 957        "--regions",
 958        type=str,
 959        required=True,
 960        help="A .npz file of input sequences, contributions, and coordinates. Can be generated using `finemo extract-regions-*` subcommands.",
 961    )
 962    call_hits_parser.add_argument(
 963        "-m",
 964        "--modisco-h5",
 965        type=str,
 966        required=True,
 967        help="A tfmodisco-lite output H5 file of motif patterns.",
 968    )
 969
 970    call_hits_parser.add_argument(
 971        "-p",
 972        "--peaks",
 973        type=str,
 974        default=None,
 975        help="DEPRECATED: Please provide this file to a preprocessing `finemo extract-regions-*` subcommand instead.",
 976    )
 977    call_hits_parser.add_argument(
 978        "-C",
 979        "--chrom-order",
 980        type=str,
 981        default=None,
 982        help="DEPRECATED: Please provide this file to a preprocessing `finemo extract-regions-*` subcommand instead.",
 983    )
 984
 985    call_hits_parser.add_argument(
 986        "-I",
 987        "--motifs-include",
 988        type=str,
 989        default=None,
 990        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.",
 991    )
 992    call_hits_parser.add_argument(
 993        "-N",
 994        "--motif-names",
 995        type=str,
 996        default=None,
 997        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.",
 998    )
 999
1000    call_hits_parser.add_argument(
1001        "-o",
1002        "--out-dir",
1003        type=str,
1004        required=True,
1005        help="The path to the output directory.",
1006    )
1007
1008    call_hits_parser.add_argument(
1009        "-t",
1010        "--cwm-trim-threshold",
1011        type=float,
1012        default=inspect.signature(call_hits)
1013        .parameters["cwm_trim_threshold_default"]
1014        .default,
1015        help="The default threshold to determine motif start and end positions within the full CWMs.",
1016    )
1017    call_hits_parser.add_argument(
1018        "-T",
1019        "--cwm-trim-thresholds",
1020        type=str,
1021        default=None,
1022        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.",
1023    )
1024    call_hits_parser.add_argument(
1025        "-R",
1026        "--cwm-trim-coords",
1027        type=str,
1028        default=None,
1029        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.",
1030    )
1031
1032    call_hits_parser.add_argument(
1033        "-l",
1034        "--global-lambda",
1035        type=float,
1036        default=inspect.signature(call_hits).parameters["lambda_default"].default,
1037        help="The default L1 regularization weight determining the sparsity of hits.",
1038    )
1039    call_hits_parser.add_argument(
1040        "-L",
1041        "--motif-lambdas",
1042        type=str,
1043        default=None,
1044        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.",
1045    )
1046    call_hits_parser.add_argument(
1047        "-a",
1048        "--alpha",
1049        type=float,
1050        default=None,
1051        help="DEPRECATED: Please use the `--lambda` argument instead.",
1052    )
1053    call_hits_parser.add_argument(
1054        "-A",
1055        "--motif-alphas",
1056        type=str,
1057        default=None,
1058        help="DEPRECATED: Please use the `--motif-lambdas` argument instead.",
1059    )
1060
1061    call_hits_parser.add_argument(
1062        "-f",
1063        "--no-post-filter",
1064        action="store_true",
1065        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.",
1066    )
1067    call_hits_parser.add_argument(
1068        "-q",
1069        "--sqrt-transform",
1070        action="store_true",
1071        help="Apply a signed square root transform to the input contributions and CWMs before hit calling.",
1072    )
1073    call_hits_parser.add_argument(
1074        "-s",
1075        "--step-size-max",
1076        type=float,
1077        default=inspect.signature(call_hits).parameters["step_size_max"].default,
1078        help="The maximum optimizer step size.",
1079    )
1080    call_hits_parser.add_argument(
1081        "-i",
1082        "--step-size-min",
1083        type=float,
1084        default=inspect.signature(call_hits).parameters["step_size_min"].default,
1085        help="The minimum optimizer step size.",
1086    )
1087    call_hits_parser.add_argument(
1088        "-j",
1089        "--step-adjust",
1090        type=float,
1091        default=inspect.signature(call_hits).parameters["step_adjust"].default,
1092        help="The optimizer step size adjustment factor. If the optimizer diverges, the step size is multiplicatively adjusted by this factor",
1093    )
1094    call_hits_parser.add_argument(
1095        "-c",
1096        "--convergence-tol",
1097        type=float,
1098        default=inspect.signature(call_hits).parameters["convergence_tol"].default,
1099        help="The tolerance for determining convergence. The optimizer exits when the duality gap is less than the tolerance.",
1100    )
1101    call_hits_parser.add_argument(
1102        "-S",
1103        "--max-steps",
1104        type=int,
1105        default=inspect.signature(call_hits).parameters["max_steps"].default,
1106        help="The maximum number of optimization steps.",
1107    )
1108    call_hits_parser.add_argument(
1109        "-b",
1110        "--batch-size",
1111        type=int,
1112        default=inspect.signature(call_hits).parameters["batch_size"].default,
1113        help="The batch size used for optimization.",
1114    )
1115    call_hits_parser.add_argument(
1116        "-d",
1117        "--device",
1118        type=str,
1119        default=None,
1120        help="DEPRECATED: Please use the `CUDA_VISIBLE_DEVICES` environment variable to specify the GPU device.",
1121    )
1122    call_hits_parser.add_argument(
1123        "-J",
1124        "--compile",
1125        action="store_true",
1126        help="JIT-compile the optimizer for faster performance. This may not be supported on older GPUs.",
1127    )
1128
1129    report_parser = subparsers.add_parser(
1130        "report",
1131        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
1132        help="Generate statistics and visualizations from hits and tfmodisco-lite motif data.",
1133    )
1134
1135    report_parser.add_argument(
1136        "-r",
1137        "--regions",
1138        type=str,
1139        required=True,
1140        help="A .npz file containing input sequences, contributions, and coordinates. Must be the same as that used for `finemo call-hits`.",
1141    )
1142    report_parser.add_argument(
1143        "-H",
1144        "--hits",
1145        type=str,
1146        required=True,
1147        help="The output directory generated by the `finemo call-hits` command on the regions specified in `--regions`.",
1148    )
1149    report_parser.add_argument(
1150        "-p",
1151        "--peaks",
1152        type=str,
1153        default=None,
1154        help="DEPRECATED: Please provide this file to a preprocessing `finemo extract-regions-*` subcommand instead.",
1155    )
1156    report_parser.add_argument(
1157        "-m",
1158        "--modisco-h5",
1159        type=str,
1160        default=None,
1161        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.",
1162    )
1163    report_parser.add_argument(
1164        "-I",
1165        "--motifs-include",
1166        type=str,
1167        default=None,
1168        help="DEPRECATED: This information is now inferred from the outputs of `finemo call-hits`.",
1169    )
1170    report_parser.add_argument(
1171        "-N",
1172        "--motif-names",
1173        type=str,
1174        default=None,
1175        help="DEPRECATED: This information is now inferred from the outputs of `finemo call-hits`.",
1176    )
1177
1178    report_parser.add_argument(
1179        "-o",
1180        "--out-dir",
1181        type=str,
1182        required=True,
1183        help="The path to the report output directory.",
1184    )
1185
1186    report_parser.add_argument(
1187        "-W",
1188        "--modisco-region-width",
1189        type=int,
1190        default=inspect.signature(report).parameters["modisco_region_width"].default,
1191        help="The width of the region around each peak summit used by tfmodisco-lite.",
1192    )
1193    report_parser.add_argument(
1194        "-t",
1195        "--cwm-trim-threshold",
1196        type=float,
1197        default=inspect.signature(report).parameters["cwm_trim_threshold"].default,
1198        help="DEPRECATED: This information is now inferred from the outputs of `finemo call-hits`.",
1199    )
1200    report_parser.add_argument(
1201        "-n",
1202        "--no-recall",
1203        action="store_true",
1204        help="Do not compute motif recall metrics.",
1205    )
1206    report_parser.add_argument(
1207        "-s",
1208        "--no-seqlets",
1209        action="store_true",
1210        help="DEPRECATED: Please omit the `--modisco-h5` argument instead.",
1211    )
1212
1213    collapse_hits_parser = subparsers.add_parser(
1214        "collapse-hits",
1215        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
1216        help="Identify best hit by motif similarity among sets of overlapping hits.",
1217    )
1218
1219    collapse_hits_parser.add_argument(
1220        "-i",
1221        "--hits",
1222        type=str,
1223        required=True,
1224        help="The `hits.tsv` or `hits_unique.tsv` file from `call-hits`.",
1225    )
1226    collapse_hits_parser.add_argument(
1227        "-o",
1228        "--out-path",
1229        type=str,
1230        required=True,
1231        help='The path to the output .tsv file with an additional "is_primary" column.',
1232    )
1233    collapse_hits_parser.add_argument(
1234        "-O",
1235        "--overlap-frac",
1236        type=float,
1237        default=inspect.signature(collapse_hits).parameters["overlap_frac"].default,
1238        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.",
1239    )
1240
1241    intersect_hits_parser = subparsers.add_parser(
1242        "intersect-hits",
1243        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
1244        help="Intersect hits across multiple runs.",
1245    )
1246
1247    intersect_hits_parser.add_argument(
1248        "-i",
1249        "--hits",
1250        type=str,
1251        required=True,
1252        nargs="+",
1253        help="One or more hits.tsv or hits_unique.tsv files, with paths delimited by whitespace.",
1254    )
1255    intersect_hits_parser.add_argument(
1256        "-o",
1257        "--out-path",
1258        type=str,
1259        required=True,
1260        help="The path to the output .tsv file. Duplicate columns are suffixed with the positional index of the input file.",
1261    )
1262    intersect_hits_parser.add_argument(
1263        "-r",
1264        "--relaxed",
1265        action="store_true",
1266        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.",
1267    )
1268
1269    args = parser.parse_args()
1270
1271    if args.cmd == "extract-regions-bw":
1272        extract_regions_bw(
1273            args.peaks,
1274            args.chrom_order,
1275            args.fasta,
1276            args.bigwigs,
1277            args.out_path,
1278            args.region_width,
1279        )
1280
1281    elif args.cmd == "extract-regions-chrombpnet-h5":
1282        extract_regions_chrombpnet_h5(
1283            args.peaks, args.chrom_order, args.h5s, args.out_path, args.region_width
1284        )
1285
1286    elif args.cmd == "extract-regions-h5":
1287        print(
1288            "WARNING: The `extract-regions-h5` command is deprecated. Use `extract-regions-chrombpnet-h5` instead."
1289        )
1290        extract_regions_chrombpnet_h5(
1291            args.peaks, args.chrom_order, args.h5s, args.out_path, args.region_width
1292        )
1293
1294    elif args.cmd == "extract-regions-bpnet-h5":
1295        extract_regions_bpnet_h5(
1296            args.peaks, args.chrom_order, args.h5s, args.out_path, args.region_width
1297        )
1298
1299    elif args.cmd == "extract-regions-modisco-fmt":
1300        extract_regions_modisco_fmt(
1301            args.peaks,
1302            args.chrom_order,
1303            args.attributions,
1304            args.sequences,
1305            args.out_path,
1306            args.region_width,
1307        )
1308
1309    elif args.cmd == "call-hits":
1310        if args.alpha is not None:
1311            warnings.warn(
1312                "The `--alpha` flag is deprecated and will be removed in a future version. Please use the `--global-lambda` flag instead."
1313            )
1314            args.global_lambda = args.alpha
1315        if args.motif_alphas is not None:
1316            warnings.warn(
1317                "The `--motif-alphas` flag is deprecated and will be removed in a future version. Please use the `--motif-lambdas` flag instead."
1318            )
1319            args.motif_lambdas = args.motif_alphas
1320
1321        call_hits(
1322            args.regions,
1323            args.peaks,
1324            args.modisco_h5,
1325            args.chrom_order,
1326            args.motifs_include,
1327            args.motif_names,
1328            args.motif_lambdas,
1329            args.out_dir,
1330            args.cwm_trim_coords,
1331            args.cwm_trim_thresholds,
1332            args.cwm_trim_threshold,
1333            args.global_lambda,
1334            args.step_size_max,
1335            args.step_size_min,
1336            args.sqrt_transform,
1337            args.convergence_tol,
1338            args.max_steps,
1339            args.batch_size,
1340            args.step_adjust,
1341            args.device,
1342            args.mode,
1343            args.no_post_filter,
1344            args.compile,
1345        )
1346
1347    elif args.cmd == "report":
1348        report(
1349            args.regions,
1350            args.hits,
1351            args.modisco_h5,
1352            args.peaks,
1353            args.motifs_include,
1354            args.motif_names,
1355            args.out_dir,
1356            args.modisco_region_width,
1357            args.cwm_trim_threshold,
1358            not args.no_recall,
1359            not args.no_seqlets,
1360        )
1361
1362    elif args.cmd == "collapse-hits":
1363        collapse_hits(args.hits, args.out_path, args.overlap_frac)
1364
1365    elif args.cmd == "intersect-hits":
1366        intersect_hits(args.hits, args.out_path, args.relaxed)