seismic cluster

Purpose

seismic cluster fits a mixture model to the filtered reads, finding subpopulations (clusters) with distinct mutation patterns. Use it when you expect your RNA preparation to contain multiple structural conformations or when DMS-MaPseq reveals heterogeneous reactivity.

Inputs

Filter output directories or report files: One or more paths to seismic filter output. See Specify Input Files.

Outputs

All outputs go into {out}/{sample}/cluster/{ref}/{reg}/.

cluster-batch-{num}.brickle: Per-cluster read assignments for one batch. See Brickle: Compressed Python Objects.
cluster-report.json: Summary of settings and the best number of clusters found. See Cluster Report.

Quick example

Cluster filter output, searching up to 3 clusters:

seismic cluster -k 3 out/sample-1/filter/ref-1

Options

Number of clusters

--max-clusters N (-k): Maximum number of clusters to try (default 0 = no limit; set this).
--min-clusters N: Minimum number of clusters to try (default 1).
--try-all-ks/--stop-best-k: Try all values of K even after finding the apparent best (default off).
--write-all-ks/--write-best-k: Write results for all K values, not just the best (default off).

EM algorithm

--min-em-runs N (-e): Minimum successful EM runs per K (default 6).
--max-em-runs N (-E): Maximum EM attempts per K (default 30).
--min-em-iter N: Minimum EM iterations per run (default 10).
--max-em-iter N: Maximum EM iterations per run (default 500).
--em-thresh F: Stop EM when log-likelihood improvement falls below F (default 0.37).

Run quality filters (used to discard low-quality EM solutions)

--max-pearson-run F: Discard runs where two clusters correlate above F (default 0.9).
--min-marcd-run F: Discard runs where two clusters are less than F apart (default 0.016).
--max-gini-run F: Discard runs where any cluster has a Gini coefficient above F (default 0.667).

K selection filters

--min-pearson-vs-best F: Remove K values where every run correlates below F with the best K (default 0.97).
--max-marcd-vs-best F: Remove K values where every run is more than F from the best K (default 0.008).

Jackpotting

--jackpot/--no-jackpot: Check for over-represented reads caused by PCR jackpotting (default on).
--max-jackpot-quotient F: Discard runs whose jackpotting quotient exceeds F (default 1.1).

Branches

--branch NAME (-b): Write outputs to {out}/{sample}/cluster_{NAME}/. See Branches.

Performance

--num-cpus N — multiprocessing; see Parallelize Tasks. --force — overwrite existing outputs.

The auto-generated Command Line Reference lists every option with its current default.

Common unexpected results

Clustering always finds 1 cluster: The data may not have enough reads or positional coverage to separate clusters. Check that --min-ninfo-pos in filter is not too strict, or that you have enough reads.
Clustering is very slow: Reduce --max-em-runs or --max-em-iter. Also reduce --num-cpus if memory is the bottleneck.