seismicrna.clusterscan package

seismicrna.clusterscan.run(input_path: Iterable[str | Path] = Sentinel.UNSET, *, branch: str = '', tmp_pfx: str | Path = './tmp', keep_tmp: bool = False, validate_gaps: bool = False, gap_min_assoc: float = 0.1, min_clusters: int = 1, max_clusters: int = 0, min_em_runs: int = 6, max_em_runs: int = 30, min_em_iter: int = 10, max_em_iter: int = 500, em_thresh: float = 0.37, min_marcd_run: float = 0.016, max_pearson_run: float = 0.9, max_arcd_vs_ens_avg: float = 0.2, max_gini_run: float = 0.667, jackpot: bool = True, jackpot_conf_level: float = 0.95, max_jackpot_quotient: float = 1.1, max_jackpot_sims: int = 12, jackpot_max_data: int = 268435456, max_loglike_vs_best: float = 0.0, min_pearson_vs_best: float = 0.97, max_marcd_vs_best: float = 0.008, try_all_ks: bool = False, write_all_ks: bool = False, cluster_pos_table: bool = True, cluster_abundance_table: bool = True, verify_times: bool = True, brotli_level: int = 10, self_contained: bool = False, num_cpus: int = 4, force: bool = False, seed: int | None = None)

Cluster the domains detected by filterscan.

Parameters:

  • branch (str) – Create a new branch of the workflow with this name [keyword-only, default: ‘’]

  • tmp_pfx (str | pathlib._local.Path) – Write all temporary files to a directory with this prefix [keyword-only, default: ‘./tmp’]

  • keep_tmp (bool) – Keep temporary files after finishing [keyword-only, default: False]

  • validate_gaps (bool) – In clusterscan, test each gap between domains for independence and merge failing gaps; off by default (expensive) [keyword-only, default: False]

  • gap_min_assoc (float) – Minimum deviation of fitted joint cluster proportions from independence to keep a gap unmerged; higher merges fewer gaps [keyword-only, default: 0.1]

  • min_clusters (int) – Start at this many clusters [keyword-only, default: 1]

  • max_clusters (int) – Stop at this many clusters (0 for no limit) [keyword-only, default: 0]

  • min_em_runs (int) – Run EM (successfully) at least this number of times for each K [keyword-only, default: 6]

  • max_em_runs (int) – Run EM (successfully or not) at most this number of times for each K [keyword-only, default: 30]

  • min_em_iter (int) – Run EM for at least this many iterations [keyword-only, default: 10]

  • max_em_iter (int) – Run EM for at most this many iterations [keyword-only, default: 500]

  • em_thresh (float) – Stop EM when the log likelihood increases by less than this threshold [keyword-only, default: 0.37]

  • min_marcd_run (float) – Remove runs with two clusters that differ by less than this MARCD [keyword-only, default: 0.016]

  • max_pearson_run (float) – Remove runs with two clusters more similar than this correlation [keyword-only, default: 0.9]

  • max_arcd_vs_ens_avg (float) – Remove runs where a cluster differs by more than this ARCD from the ensemble average at any position [keyword-only, default: 0.2]

  • max_gini_run (float) – Remove runs where any cluster’s Gini coefficient exceeds this limit [keyword-only, default: 0.667]

  • jackpot (bool) – Calculate the jackpotting quotient to find over-represented reads [keyword-only, default: True]

  • jackpot_conf_level (float) – Confidence level for the jackpotting quotient confidence interval [keyword-only, default: 0.95]

  • max_jackpot_quotient (float) – Remove runs whose jackpotting quotient exceeds this limit [keyword-only, default: 1.1]

  • max_jackpot_sims (int) – Maximum number of simulations to compute the jackpotting quotient [keyword-only, default: 12]

  • jackpot_max_data (int) – Skip calculating the jackpotting quotient if the reads × positions exceeds this limit [keyword-only, default: 268435456]

  • max_loglike_vs_best (float) – Remove Ks with a log likelihood gap larger than this (0 for no limit) [keyword-only, default: 0.0]

  • min_pearson_vs_best (float) – Remove Ks where every run has less than this correlation vs. the best [keyword-only, default: 0.97]

  • max_marcd_vs_best (float) – Remove Ks where every run has more than this MARCD vs. the best [keyword-only, default: 0.008]

  • try_all_ks (bool) – Try all numbers of clusters (Ks), even after finding the best number [keyword-only, default: False]

  • write_all_ks (bool) – Write all numbers of clusters (Ks), rather than only the best number [keyword-only, default: False]

  • cluster_pos_table (bool) – Tabulate relationships per position for cluster data [keyword-only, default: True]

  • cluster_abundance_table (bool) – Tabulate number of reads per cluster for cluster data [keyword-only, default: True]

  • verify_times (bool) – Verify that report files from later steps have later timestamps [keyword-only, default: True]

  • brotli_level (int) – Compress pickle files with this level of Brotli (0 - 11) [keyword-only, default: 10]

  • self_contained (bool) – Write self-contained batch files that skip loading predecessor batches (Filter/Cluster steps); costs more disk space [keyword-only, default: False]

  • num_cpus (int) – Use up to this many CPUs simultaneously [keyword-only, default: 4]

  • force (bool) – Force all tasks to run, overwriting any existing output files [keyword-only, default: False]

  • seed (int | None) – Seed for the random number generator [keyword-only, default: None]

Subpackages

Submodules

Validate the gap between two adjacent filterscan domains by testing whether the domains’ clusters are independent on the reads that span the gap.

filterscan cuts where anti-correlated “bridge” pairs are depleted, but that depletion happens both at a genuine independent gap and across a positively-correlated (co-folding) region – so a spurious cut can fragment one module into two. This test distinguishes the two cases: it reuses the already-computed per-domain cluster mutation-rate profiles (the expensive part) as fixed components and fits only the joint cluster proportions on the spanning reads – a concave, stable, fast reduced EM (no re-clustering, none of the merged-EM instability). It then compares, by BIC, an independent (product) model against a free-association (saturated) model. If association genuinely improves the fit, the gap is spurious and the domains should be merged.

Because a spanning read’s likelihood factorizes into an A-side part (its A-domain positions under the A-profiles) and a B-side part, each read’s component log-likelihood is just A_i(a) + B_i(b) and both models are fit over the same fixed components.

seismicrna.clusterscan.gap.evaluate_gap(merged_uniq_reads: UniqReads, mus_a: pd.DataFrame, mus_b: pd.DataFrame, left_end3: int, right_end5: int, gap_min_assoc: float) → tuple[bool, dict]

Test whether the gap between domain A and domain B is genuine.

Parameters:

  • merged_uniq_reads (UniqReads) – Unique reads over the merged (A ∪ B) region.

  • mus_a (pd.DataFrame) – Cluster mutation rates for domain A and domain B, indexed by position (1-indexed) with one column per cluster.

  • mus_b (pd.DataFrame) – Cluster mutation rates for domain A and domain B, indexed by position (1-indexed) with one column per cluster.

  • left_end3 (int) – Last position of domain A.

  • right_end5 (int) – First position of domain B. May be greater than left_end3 + 1: a gap-mode of “omit” (the default) leaves an uncalled background region between adjacent domains, rather than domain B starting immediately after domain A.

  • gap_min_assoc (float) – Minimum deviation of the fitted joint proportions from independence (max |pi_sat - outer(u, v)|) required to call a gap spurious.

Returns:

(passed, info) – passed is True to keep the gap, False to merge the domains. info holds the fitted statistics for logging.

Return type:

tuple[bool, dict]

class seismicrna.clusterscan.io.ClusterScanFile

Bases: HasRegFilePath, ABC

classmethod get_step()

Step of the workflow.

class seismicrna.clusterscan.io.ClusterScanIO

Bases: ClusterScanFile, RegFileIO, ABC

seismicrna.clusterscan.main.run(input_path: Iterable[str | Path] = Sentinel.UNSET, *, branch: str = '', tmp_pfx: str | Path = './tmp', keep_tmp: bool = False, validate_gaps: bool = False, gap_min_assoc: float = 0.1, min_clusters: int = 1, max_clusters: int = 0, min_em_runs: int = 6, max_em_runs: int = 30, min_em_iter: int = 10, max_em_iter: int = 500, em_thresh: float = 0.37, min_marcd_run: float = 0.016, max_pearson_run: float = 0.9, max_arcd_vs_ens_avg: float = 0.2, max_gini_run: float = 0.667, jackpot: bool = True, jackpot_conf_level: float = 0.95, max_jackpot_quotient: float = 1.1, max_jackpot_sims: int = 12, jackpot_max_data: int = 268435456, max_loglike_vs_best: float = 0.0, min_pearson_vs_best: float = 0.97, max_marcd_vs_best: float = 0.008, try_all_ks: bool = False, write_all_ks: bool = False, cluster_pos_table: bool = True, cluster_abundance_table: bool = True, verify_times: bool = True, brotli_level: int = 10, self_contained: bool = False, num_cpus: int = 4, force: bool = False, seed: int | None = None)

Cluster the domains detected by filterscan.

Parameters:

  • branch (str) – Create a new branch of the workflow with this name [keyword-only, default: ‘’]

  • tmp_pfx (str | pathlib._local.Path) – Write all temporary files to a directory with this prefix [keyword-only, default: ‘./tmp’]

  • keep_tmp (bool) – Keep temporary files after finishing [keyword-only, default: False]

  • validate_gaps (bool) – In clusterscan, test each gap between domains for independence and merge failing gaps; off by default (expensive) [keyword-only, default: False]

  • gap_min_assoc (float) – Minimum deviation of fitted joint cluster proportions from independence to keep a gap unmerged; higher merges fewer gaps [keyword-only, default: 0.1]

  • min_clusters (int) – Start at this many clusters [keyword-only, default: 1]

  • max_clusters (int) – Stop at this many clusters (0 for no limit) [keyword-only, default: 0]

  • min_em_runs (int) – Run EM (successfully) at least this number of times for each K [keyword-only, default: 6]

  • max_em_runs (int) – Run EM (successfully or not) at most this number of times for each K [keyword-only, default: 30]

  • min_em_iter (int) – Run EM for at least this many iterations [keyword-only, default: 10]

  • max_em_iter (int) – Run EM for at most this many iterations [keyword-only, default: 500]

  • em_thresh (float) – Stop EM when the log likelihood increases by less than this threshold [keyword-only, default: 0.37]

  • min_marcd_run (float) – Remove runs with two clusters that differ by less than this MARCD [keyword-only, default: 0.016]

  • max_pearson_run (float) – Remove runs with two clusters more similar than this correlation [keyword-only, default: 0.9]

  • max_arcd_vs_ens_avg (float) – Remove runs where a cluster differs by more than this ARCD from the ensemble average at any position [keyword-only, default: 0.2]

  • max_gini_run (float) – Remove runs where any cluster’s Gini coefficient exceeds this limit [keyword-only, default: 0.667]

  • jackpot (bool) – Calculate the jackpotting quotient to find over-represented reads [keyword-only, default: True]

  • jackpot_conf_level (float) – Confidence level for the jackpotting quotient confidence interval [keyword-only, default: 0.95]

  • max_jackpot_quotient (float) – Remove runs whose jackpotting quotient exceeds this limit [keyword-only, default: 1.1]

  • max_jackpot_sims (int) – Maximum number of simulations to compute the jackpotting quotient [keyword-only, default: 12]

  • jackpot_max_data (int) – Skip calculating the jackpotting quotient if the reads × positions exceeds this limit [keyword-only, default: 268435456]

  • max_loglike_vs_best (float) – Remove Ks with a log likelihood gap larger than this (0 for no limit) [keyword-only, default: 0.0]

  • min_pearson_vs_best (float) – Remove Ks where every run has less than this correlation vs. the best [keyword-only, default: 0.97]

  • max_marcd_vs_best (float) – Remove Ks where every run has more than this MARCD vs. the best [keyword-only, default: 0.008]

  • try_all_ks (bool) – Try all numbers of clusters (Ks), even after finding the best number [keyword-only, default: False]

  • write_all_ks (bool) – Write all numbers of clusters (Ks), rather than only the best number [keyword-only, default: False]

  • cluster_pos_table (bool) – Tabulate relationships per position for cluster data [keyword-only, default: True]

  • cluster_abundance_table (bool) – Tabulate number of reads per cluster for cluster data [keyword-only, default: True]

  • verify_times (bool) – Verify that report files from later steps have later timestamps [keyword-only, default: True]

  • brotli_level (int) – Compress pickle files with this level of Brotli (0 - 11) [keyword-only, default: 10]

  • self_contained (bool) – Write self-contained batch files that skip loading predecessor batches (Filter/Cluster steps); costs more disk space [keyword-only, default: False]

  • num_cpus (int) – Use up to this many CPUs simultaneously [keyword-only, default: 4]

  • force (bool) – Force all tasks to run, overwriting any existing output files [keyword-only, default: False]

  • seed (int | None) – Seed for the random number generator [keyword-only, default: None]

class seismicrna.clusterscan.report.ClusterScanReport(**kwargs: Any | Callable[[Report], Any])

Bases: RegReport, ClusterScanIO

classmethod get_file_seg_type()

Type of the last segment in the path.

classmethod get_result_report_fields()

Result fields of the report.

seismicrna.clusterscan.write.clusterscan(filterscan_report_file: Path, *, branch: str, tmp_pfx: str | Path, keep_tmp: bool, brotli_level: int, force: bool, num_cpus: int, validate_gaps: bool, gap_min_assoc: float, min_clusters: int, max_clusters: int, min_em_runs: int, max_em_runs: int, jackpot: bool, jackpot_conf_level: float, max_jackpot_quotient: float, max_jackpot_sims: int, jackpot_max_data: int, min_em_iter: int, max_em_iter: int, em_thresh: float, min_marcd_run: float, max_pearson_run: float, max_arcd_vs_ens_avg: float, max_gini_run: float, max_loglike_vs_best: float, min_pearson_vs_best: float, max_marcd_vs_best: float, try_all_ks: bool, write_all_ks: bool, cluster_pos_table: bool, cluster_abundance_table: bool, verify_times: bool, self_contained: bool, seed: int | None)

Cluster the domains detected by filterscan for one FilterScanReport. If validate_gaps is set, also validate every gap between adjacent domains: test on the reads spanning the gap whether the two domains’ clusters are independent. If a gap fails (the domains are coupled), merge the two domains, re-cluster the merged region, and repeat until every remaining gap passes. Write a ClusterScanReport recording the final domains, their numbers of clusters, and the original filterscan domains comprising each final domain.