convert: write SCope looms with --scope (#61)

docs/COMMANDS.md

@@ -131,7 +131,7 @@ adata convert cistopic_obj.pkl -o atac.h5ad
 | `cellbender` | `.h5` | CellBender `remove-background` 0.2 and later | no |
 | `10x-mtx` | a folder (`--to 10x-mtx`) | Cell Ranger 2 and 3+, STARsolo (velocyto too) | Cell Ranger 3+ layout, gzipped |
 | `parse` | a folder | Parse Biosciences split-pipe: a run, a sample or a `DGE_*` folder | no |
-| `loom` | `.loom` | loompy, velocyto, kallisto\|bustools, SCope (pySCENIC, SCENIC+, pycistopic); spec 2 and 3 | spec 3.0.0, as loompy 3 writes it |
+| `loom` | `.loom` | loompy, velocyto, kallisto\|bustools, SCope (pySCENIC, SCENIC+, pycistopic); spec 2 and 3 | spec 3.0.0, as loompy 3 writes it; a SCope loom with `--scope` |
 | `text` | `.csv`, `.tsv`, `.tab`, `.txt`, `.data`, each also `.gz` | a dense matrix, as anndata's `read_csv` and `read_text` read it | no |
 | `hdf5` | `.hdf5`, `.he5`, `.hdf`; a plain `.h5` with `--from hdf5` | one 2-D dataset, as anndata's `read_hdf` reads it | no |
 | `cistopic` | `.pkl`, `.pickle`, gzipped too | a pycisTopic `CistopicObject` or `CistopicImputedFeatures`, without pycisTopic (needs `pyadata-cli[pickle]`) | no |
@@ -449,8 +449,56 @@ where genes go, so var is named by region, and SCENIC+'s `linkedGene` is a
 var column. The compound attributes are read a block of cells at a time,
 so memory is one block across all regulons, never the table.
 
-Writing a loom writes a plain one: an AnnData's regulons, embeddings and
-clusterings are not yet written in SCope's form (#61).
+#### Writing SCope looms
+
+`--scope` writes a loom the SCope viewer reads: an AnnData's regulons,
+embeddings and clusterings in SCope's form, with a `MetaData` naming them.
+Without it a loom export is plain, and one of an AnnData that came from a
+SCope loom says that `--scope` would keep its regulons.
+
+```bash
+adata convert scenic.h5ad -o scenic.loom --scope
+adata convert scenic.h5ad -o scenic.loom --scope --clusterings leiden,cell_type
+adata convert scenic.h5ad -o scenic.loom --scope --default-embedding X_umap
+```
+
+| AnnData | SCope loom |
+|---|---|
+| a dataframe in obsm or varm (`RegulonsAUC`, `Regulons`, `ClusterMarkers_<id>`) | a compound `col_attrs` or `row_attrs` of its numeric columns |
+| `obsm` arrays of two columns | `Embeddings_X` and `Embeddings_Y`, a field per id; the default (id -1) as `Embedding` too |
+| the categorical obs columns a SCope import recorded, and those `--clusterings` names | `Clusterings`, int codes per clustering (and not a plain attribute too), with `ClusterID` for the first; their clusters described in `MetaData` |
+| other obs columns | plain attributes, as ever, listed in `MetaData` as metrics (numbers) and annotations (up to 245 values) |
+| `uns["scope"]` | `MetaData`'s regulon thresholds (from the `regulon_thresholds` table, as edited), and the recorded ids and names |
+
+| Flag | Meaning |
+|---|---|
+| `--scope` | Write a SCope loom |
+| `--clusterings` | Categorical obs columns to write as clusterings, comma separated |
+| `--default-embedding` | The obsm key SCope shows first (default: the recorded default, else `X_umap`, else `X_tsne`, else the first) |
+
+An AnnData that came from a SCope loom keeps its embeddings' and
+clusterings' ids and names, recorded in `uns["scope"]`, so a SCope loom
+converts to h5ad and back with the same attributes. As SCope expects,
+clusterings and their clusters are numbered 0..n-1 in order, regulon names
+lose their spaces (`TP53 (5g)` becomes `TP53_(5g)`, alike in the AUC and
+the Regulons, with a warning), and cluster markers are written only as
+`ClusterMarkers_<id>` with their metrics (`ClusterMarkers_<id>_avg_logFC` and the like, which
+`clusterMarkerMetrics` then lists, under the names a SCope import recorded), their fields renumbered with the clusters; a
+clustering `--clusterings` names owns `ClusterMarkers_<its id>`. `ClusterID` is always written afresh from the first
+clustering, in place of an obs column of that name, so the two cannot
+disagree. `MetaData` is plain JSON, which SCope reads as it
+reads compressed. A clustering that is not categorical is refused, as is a
+table with more columns than HDF5 can hold in one compound attribute (its
+type must fit 64 KiB, a limit pySCENIC's own export shares), before any
+output exists. Attributes are written a block of cells at a time. A loom's
+attributes share one flat namespace per axis, with or without `--scope`:
+where an obs column, an obsm entry or one of SCope's parts ask for one name,
+the first keeps it (loom's and SCope's own attributes, then columns, SCope's
+tables, arrays) and the others are left out with a warning.
+
+SCope opens looms without loompy's validation, which rejects compound
+attributes; so do `loompy.connect(..., validate=False)` and the pySCENIC
+tutorial's recipe, which CI uses to read our SCope looms back.
 
 ### Text tables
 

docs/TESTING.md

@@ -36,7 +36,7 @@ separate jobs.
 | `test_cistopic.py` | pycisTopic pickles read through the allowlist: a CistopicObject's matrices, frames, model and projections, an imputed-features pickle, a pickle naming `os.system` refused unrun, the missing-extra hint |
 | `test_pandas_frame.py` | Each pandas dtype written as anndata writes it, on every backend |
 | `test_concat_axis_elements.py` | `concat` on obsm dataframes and varm, against `anndata.concat` for each join and merge strategy |
-| `test_scope.py` | SCope looms as pySCENIC's `export2loom` writes them: regulons as dataframes, embeddings and clusterings named from `MetaData`, plain or compressed `MetaData`, region looms, R's logicals and factors, renamed fields, small blocks |
+| `test_scope.py` | SCope looms, read as pySCENIC's `export2loom` writes them and written with `--scope` (round trips, plain AnnData, refusals, markers): regulons as dataframes, embeddings and clusterings named from `MetaData`, plain or compressed `MetaData`, region looms, R's logicals and factors, renamed fields, small blocks |
 | `test_table.py` | The streaming CSV reader against pandas: column types, quoting, row subsets, blocks, refusals |
 | `test_mtx_codec.py`, `test_csr_build.py` | The streaming Matrix Market reader and writer against scipy, and the bucket sort they and `recast` share |
 | `test_external_readers.py` | 10x conversions, HDF5 and MTX, Parse runs and looms, checked by scanpy (Parse's own recipe for Parse), anndata with loompy (the pySCENIC tutorial's recipe for SCope looms), anndata's `read_csv`, `read_text` and `read_hdf`, a pandas 1.5 pycisTopic pickle read back there, and PyTables in their own uv environments (`integration`) |

src/adata/cli.py

@@ -537,6 +537,9 @@ def concat(
     "key": "--key",
     "obs_names_key": "--obs-names-key",
     "var_names_key": "--var-names-key",
+    "scope": "--scope",
+    "clusterings": "--clusterings",
+    "default_embedding": "--default-embedding",
 }
 
 #: Flags that only make sense for format conversion; with recast's flags,
@@ -738,6 +741,24 @@ def convert(
         help="HDF5 input: the dataset naming the cells (default row_names, "
         "col_names with --transpose)",
     ),
+    scope: bool = typer.Option(
+        False,
+        "--scope",
+        help="loom output: write a SCope loom: regulons, embeddings, clusterings "
+        "and MetaData, as the SCope viewer reads them",
+    ),
+    clusterings: Optional[str] = typer.Option(
+        None,
+        "--clusterings",
+        help="SCope output: categorical obs columns written as clusterings, "
+        "comma separated (besides those a SCope import recorded)",
+    ),
+    default_embedding: Optional[str] = typer.Option(
+        None,
+        "--default-embedding",
+        help="SCope output: the obsm key SCope shows first (default: the "
+        "recorded default, else X_umap, else X_tsne, else the first)",
+    ),
     var_names_key: Optional[str] = typer.Option(
         None,
         "--var-names-key",
@@ -773,7 +794,7 @@ def convert(
     is read as anndata.read_loom reads it; a SCope loom (pySCENIC, SCENIC+,
     pycistopic) also gives its regulons as dataframes in obsm and varm, its
     embeddings and clusterings named from its MetaData, and the MetaData in
-    uns["scope"]. Features are named by gene id, unique but in probe-level
+    uns["scope"]; --scope writes them back as SCope reads them. Features are named by gene id, unique but in probe-level
     (Flex) matrices; --var-names gene_symbols names them as scanpy does.
 
     Nothing is written until every check has passed, and the output is moved
@@ -800,6 +821,7 @@ def convert(
         adata convert velocyto.loom -o velocyto.h5ad
         adata convert pyscenic_output.loom -o scenic.h5ad
         adata convert counts.h5ad -o counts.loom
+        adata convert scenic.h5ad -o scenic.loom --scope --clusterings leiden
         adata convert GSE123_counts.txt.gz -o counts.h5ad --transpose
         adata convert matrix.hdf5 -o matrix.h5ad --key data/counts
         adata convert cistopic_obj.pkl -o atac.h5ad
@@ -900,6 +922,9 @@ def convert(
                 key=key,
                 obs_names_key=obs_names_key,
                 var_names_key=var_names_key,
+                scope=scope,
+                clusterings=clusterings,
+                default_embedding=default_embedding,
             ),
         )
     except Exception as e:

src/adata/conversion/base.py

@@ -59,6 +59,11 @@ class Options:
     obs_names_key: Optional[str] = None
     var_names_key: Optional[str] = None
 
+    # Writing a SCope loom.
+    scope: bool = False
+    clusterings: Optional[str] = None
+    default_embedding: Optional[str] = None
+
 
 @dataclass
 class ConversionPlan: