piaso-tools 1.2.5

pip install -U piaso-tools

runSVD on a cytome runs its passes in Rust, and converges by default

The streaming SVD's cost is the pass over the store: one sparse product per pass, twenty-odd passes per solve. Those passes now run in a Rust engine that reads the chunks straight out of the file, in whichever layout the matrix is stored (row-major or column-major, any integer width or float, and the zlib blobs older files carry under the zstd label), blocks each chunk against the dense operand so every random access stays in cache, and lets each thread own its slice of the output. Every entry of every product is one sum in a fixed order, so the result does not depend on the thread count and reruns are bit-identical. Measured on a 37,609-cell, 317,726-peak matrix, one round of the solver went from 27 s to 9.5 s; on a million-cell tile matrix, from 199 s to 135 s. n_threads=0 (the default) uses every core the process is allowed; PIASO_SVD_PROFILE=1 prints, per pass, how long was spent reading, bucketing and multiplying.

The solver changed with it. method="auto", now the default, runs block Lanczos with a stopping rule (tol=1e-3) and reports what it did:

Solver: block Lanczos (at most 12 rounds = 26 passes)
Solver done: block Lanczos, 9 rounds, 20 passes (converged)
Solver time: 43.0 s, of which 41.6 s in 20 passes over the store and 1.4 s in numpy between them

The old power iteration at its default of seven rounds was not converged on a peak matrix; it needs about twenty-two. method="power" with n_iter still runs exactly that many rounds when you want a fixed budget, and method="krylov" with n_iter does the same for Lanczos.

Two more defaults for peak and tile matrices on a cytome. TF-IDF is applied on the fly to raw ATAC or tile counts unless auto_tfidf=False; the pass engine weights each entry as it reads it, so no normalised layer is written. The guard that refuses to normalise an already-normalised layer runs first: inferred TF-IDF steps aside with a warning on a layer that does not look like raw counts, and an explicit auto_tfidf=True raises. And the reads take no file locks: once the write-ahead log has been checkpointed, the engine opens the store as an immutable snapshot, which is what lets it run at full speed from a network file system. On a cluster node reading from a network mount, a twenty-pass solve went from 5 min 13 s to 44.9 s; from the node's local disk the same solve is 29 s. The installation page has a section on running on a cluster.

leiden_local re-embeds every group with runSVD

dr_method="X_svd" is the new default: each group is embedded by runSVD on its own selected features and clustered from there. On a cytome the per-group work streams from disk for any modality, with the TF-IDF runSVD applies to peaks or tiles and the store's highly_variable selection, and each group's subset carries only the one matrix it needs. "X_svd_full" uses every feature rather than the selected ones; "X_pca" keeps the INFOG route. min_cells, svd_method, n_iter and svd_tol pass through, and the function no longer touches scanpy. infog_svd likewise defaults to the infog layer, and where there is no highly_variable column it uses every gene and says so.

Feature selection by name, and two primate genomes

piaso.pp.selectFeatures(ds, features) flags the named features (peak ids, gene symbols, Ensembl ids, whatever the cytome stores) in one call, writing the boolean column runSVD and leiden_local read by default. match="all" and min_matched= turn a partial match into an error.

piaso.data.fetch_genome knows rhesus macaque (Mmul_10) and common marmoset (mCalJac1), with chromosome names as the Cell Ranger ARC references carry them; piaso.data.makeGenomeFiles builds the same files for any assembly from a GTF and a chromosome-sizes file. A genome without a GTF source no longer announces one.

Every reader takes a column-major matrix

Peak and tile matrices stored column-major (the layout the SVD engine prefers; cytome 0.3.3 reads and subsets it as such) are read that way throughout: score, calculateCellMetrics, calculateFeatureMetrics, normalize_log1p, the ligand-receptor readers, the TF-IDF raw-count guard and the feature reads behind plotting all take the columns when that is what the store holds, instead of asking for row chunks a column-major matrix does not have.

A breaking change: getMarkers(as_dict=True) returns the dictionary

It used to return (DataFrame, dict). Almost every use of it wanted the dict — it is the form predictCellTypeByMarker takes — and paid for that by unpacking a pair and discarding the first element. Asking for a dict now gets a dict:

marker_db = piaso.tl.getMarkers(study="AllenWholeMouseBrain_isocortex",
                                as_dict=True)
piaso.tl.predictCellTypeByMarker(adata, marker_gene_set=marker_db)

as_dict='both' returns the old pair. It is deprecated and will be removed in the next release. as_dict=False, the default, is unchanged.

The change announces itself once per session on the first as_dict=True, because a caller that keeps unpacking does not always fail loudly: with exactly two cell types, df, d = getMarkers(as_dict=True) succeeds and binds two cell-type names. Search your code for , ... = piaso.tl.getMarkers( — the tutorials are updated.

Dot plots: the same numbers, two orders of magnitude faster

piaso.pl.dotplot computed each group's mean and expressed fraction with one Python comparison per (feature, group) pair, each scanning every cell. It now builds a one-hot group indicator once and takes two sparse products. The numbers are identical to 5e-16. A 100,000-cell, 60-feature, 30-group plot goes from 68 s to under half a second; at 20,000 cells it is 9.1 s to 0.05 s.

A missing groupby column is an error, not a wrong answer

SQLite has a compatibility misfeature: a double-quoted token that does not resolve to a column is treated as a string literal rather than raising. So SELECT "Leiden" FROM cells on a cytome with no Leiden column returns the string 'Leiden' for every row — and the caller sees one group, named after the typo, with a marker table that looks entirely plausible. That is how a COSG run reported 1 groups for a column that was never there.

Every place a user-supplied column name reaches such a query now looks it up first: cosg.run_cosg_cytome, piaso.tl.runSCALAR, and the new fragment-length functions. The error lists the columns that do exist. A case variant is still accepted, because SQLite column names are case-insensitive and leiden really does name Leiden. piaso.utils.require_cells_column is the shared check.

groups= draws a legend of what is on the plot

plotEmbedding(color='CellTypes', groups=['A','B']) draws A and B in their palette colours and everything else as one grey cloud. The legend was built from the full category list, so on a forty-cluster column it listed forty entries in forty colours, thirty-eight of which appeared nowhere on the figure — while the grey covering most of the points had no entry at all.

It now lists the highlighted groups and one final entry for the grey, labelled by the new na_label (default 'other') and shown only when something was actually left out. Highlighted groups keep the colour they have on an unfiltered plot.

color= can be the values, not only a column name

Anything computed on the fly — a regulon activity vector, a per-cell score never written back, a mask — had to be assigned into obs first. Passing the array reached color in adata.obs.columns and came back as TypeError: unhashable type: 'numpy.ndarray'.

plotEmbedding and plot_embeddings_split now accept a numpy array, a pandas Series or a Categorical as one value per cell. Numeric values get a colorbar, everything else a category legend; a Categorical keeps its own order, and a named Series titles the panel. A wrong length names both numbers. A plain list/tuple still means one panel per colour.

A Series is aligned by name when its index resolves to the cell names — anything that came out of a groupby, a sort or a .loc colours the right cells instead of being silently permuted. A RangeIndex, or an index sharing no labels with the cells, is used positionally; an index that resolves some of its labels and not others is refused, naming the count, because neither reading is safe. Duplicate cell names make alignment ambiguous, so those fall back to position with a warning.

predictCellTypeByMarker can weight its markers

A marker list treats a textbook marker and a weak, broadly expressed one identically. marker_gene_weights= passes per-gene weights through to piaso.tl.score(gene_weights=): a dict keyed by cell type, a DataFrame read column by column, or a sequence in set order.

The weights are used as given, with no normalisation per cell type. score() divides by median(weights) * n_genes and builds its control sets from the same weights, so rescaling one cell type's weights cannot change which cell type wins — measured at 9e-16 across a 1000× rescale. Normalising would be a step that provably changes nothing.

What a rescale cannot survive is a median of zero: the divisor becomes zero and every score becomes NaN without a word. score() now rejects that, and non-finite weights, naming the gene set. Mismatched lengths, a missing cell type and an unknown one are each refused with both numbers.

Smaller things

• predictCellTypeByGDR fits on the shared label column. With a query that lacks the reference's label column, the classifier was fitted on a column that only the reference rows carried and failed with a KeyError; it now fits on the combined column both sides share.

• Twenty-two public functions have complete docstrings, checked by a test against their signatures: every parameter present, every default the one the code uses, the cytome-only parameters and the deprecated aliases documented where they live. The remaining functions are on the same test's backlog.

• leiden_local stops warning about its own work. It called neighbors without passing the modality on, so an ATAC run warned on every call that it had found the embedding PIASO itself had just written under another modality. The argument is forwarded; graphs and labels are unchanged, and a genuine cross-modality guess still warns.

• A per-cell column will not land on someone else's. SQLite column names are case-insensitive, and cytome resolves a new spelling to an existing column with only a warning — so writing tss_score into a cytome already carrying a TSS_score column from another tool replaced it. piaso.utils.write_cells_column is the shared guard: the same spelling overwrites (re-running a QC step replaces its own column), a spelling that differs only in case is refused unless overwrite=True.

• The TF-IDF cache no longer grows with the dataset. compute_tfidf_stats kept its per-cell depth and per-feature IDF vectors in ds.metadata, which is a JSON document rewritten whole on every save — large enough on a real ATAC cytome for cytome to warn about its size. The depth is now a cells column, the IDF stays the feature column it was already written to, and metadata keeps the scale factor and the two column names. At 20K cells by 30K peaks that is 98 bytes instead of 1.0 MB. Cytomes written the old way are read unchanged.