Single-Cell RNA-seq Analysis

A single-cell count matrix records how many RNA molecules from each gene were observed in each cell. Most entries are zero, so BioLang keeps 10x matrices sparse through filtering, normalization, variable-gene selection, and PCA.

Load 10x Data

sc.load reads the Matrix Market file, gene features, and cell barcodes. The result keeps raw counts in both matrix and layers.counts, with cell metadata in obs and gene metadata in var.

import "singlecell" as sc

let cells = sc.load("filtered_feature_bc_matrix")
println(str(cells.n_cells) + " cells x " + str(cells.n_genes) + " genes")
println(str(nnz(cells.matrix)) + " non-zero counts")

AnnData exchange

BioLang reads and writes AnnData Zarr stores natively and preserves a sparse X matrix plus observation and variable index names. Arbitrary metadata columns and auxiliary layers are not yet copied. Direct HDF5 .h5ad access requires Python/anndata or a container conversion to Zarr.

let cells = read_anndata("study.zarr")
write_anndata("study_biolang.zarr", cells)

Quality Control

Low-complexity cells may be empty droplets. Very high gene counts may indicate two cells captured together. A high mitochondrial fraction often signals a damaged cell. Thresholds depend on tissue and protocol, so inspect the distributions before fixing them.

let metrics = cell_qc(cells.matrix, cells.genes)
println("median genes per cell: " + str(median(col(metrics, "n_genes"))))
println("median mitochondrial %: " + str(median(col(metrics, "pct_mito"))))

let filtered = cells
    |> sc.filter_genes(3)
    |> sc.filter_cells(200, 5000, 20.0)

println(str(filtered.n_cells) + " cells retained")

Normalize and Select Genes

Library-size normalization makes cells with different sequencing depths comparable. log1p reduces the influence of very large counts. Variable-gene selection focuses the later analysis on genes that distinguish cells rather than genes that are nearly constant.

let prepared = filtered
    |> sc.normalize(10000.0)
    |> sc.variable_genes(2000)

println(str(len(prepared.hvg_genes)) + " variable genes selected")

PCA, Neighbors, and Leiden Clusters

PCA compresses thousands of genes into a small set of components. Neighbors are computed from those PCA scores. Leiden then finds communities in the stored neighbor graph; it does not cluster the raw expression matrix or rebuild the graph. The current exact neighbor search takes quadratic time in the number of cells, so very large atlases should be sampled or sent to a suitable remote backend.

let reduced = prepared |> sc.run_pca(30)
let graphed = reduced |> sc.neighbors(15)
let result = graphed |> sc.cluster_leiden(15, 0.5)

println("PCA components: " + str(result.pca.n_components))
println("neighbor edges: " + str(len(result.knn)))
println("clusters: " + str(len(unique(result.clusters))))

Resolution sensitivity

Cluster IDs are labels, not biological identities. Higher resolution usually creates more clusters. Compare several values and check markers, sample composition, and stability rather than selecting a number only because the plot looks tidy.

for resolution in [0.2, 0.5, 0.8, 1.0] {
    let labels = leiden_graph(graphed.knn, graphed.n_cells, resolution)
    println("resolution " + str(resolution) + ": " +
            str(len(unique(labels))) + " clusters")
}

Markers and Annotation

A marker is a gene whose expression differs between groups. Marker genes support a cell-type interpretation, but no single marker is definitive. Use multiple genes and biological context.

log2fc is positive when the gene is higher in the first cluster, matching Seurat's FindMarkers(ident.1, ident.2). The pct_a and pct_b columns give the fraction of cells expressing the gene in each cluster, which is the specificity evidence a mean hides.

let markers = sc.marker_table(result, 0, 1)
let strongest = markers
    |> filter(|row| row.padj < 0.05)
    |> sort(|a, b| if abs(a.log2fc) > abs(b.log2fc) { -1 } else { 1 })
    |> take(10)

for marker in strongest {
    println(marker.gene + " log2FC=" + str(marker.log2fc) +
            " FDR=" + str(marker.padj))
}

Plots

UMAP is a view of local relationships, not proof of discrete cell types. BioLang computes it from PCA scores and can color points using the stored Leiden labels. The plotting API returns SVG text, so the same result can be rendered inline in Desktop and notebooks or written as a publication-ready vector file.

write_text("umap.svg", sc.plot_umap(result, "PBMC clusters"))
write_text("pca.svg", sc.plot_pca(result, "PBMC PCA"))
write_text("cd3d.svg", sc.plot_feature(result, "CD3D", "CD3D expression"))
write_text("qc.svg", sc.plot_qc_dashboard(cells, "Cell quality"))

Diagnostic and comparison plots

The advanced plotting module covers the common views needed to inspect quality, embeddings, differential expression, donor consistency, cell composition, cluster stability, and enrichment. Split feature plots and group heatmaps help distinguish a condition effect from a change in cell abundance.

# conditions and cell_types contain one label per cell.
write_text("condition.svg",
    sc.plot_embedding(result, conditions, "Condition"))
write_text("ifit1_by_condition.svg",
    sc.plot_feature_split(result, "IFIT1", conditions, "IFIT1"))
write_text("group_means.svg",
    sc.plot_group_heatmap(result, ["IFIT1", "ISG15", "EPCAM"],
                          conditions, "Mean expression"))

Donor-Aware Advanced Analysis

Cells from one donor are repeated observations, not independent biological replicates. Aggregate raw counts by donor, condition, and cell type before donor-level differential testing. This avoids treating thousands of cells as thousands of independent patients.

let profiles = sc.pseudobulk_profiles(
    result, donors, conditions, cell_types, 10
)

let de = sc.paired_pseudobulk_de(
    result, donors, conditions, cell_types, "T cell",
    "control", "treated", 10
)

write_text("volcano.svg", sc.plot_volcano(de, "T-cell response"))
write_text("ma.svg", sc.plot_ma(de, "T-cell mean-abundance"))
write_text("donors.svg",
    sc.plot_donor_expression(
        result, "IFIT1", donors, conditions, cell_types, "T cell",
        "control", "treated", 10
    ))
write_text("pseudobulk_pca.svg",
    sc.plot_pseudobulk_pca(
        result, donors, conditions, cell_types, "T cell", 10
    ))

paired_pseudobulk_de is an exploratory paired log2-CPM t-test. Use its effect sizes and donor plots for exploration. For confirmatory claims, export the raw pseudobulk counts and fit a negative-binomial model with DESeq2, edgeR, or an equivalent validated method.

Composition and clustering diagnostics

let composition = sc.composition_table(
    donors, conditions, cell_types
)
let composition_test = sc.paired_composition_test(
    donors, conditions, cell_types, "control", "treated", 20
)
write_text("composition.svg",
    sc.plot_composition(
        donors, conditions, cell_types, "Cell-type composition"
    ))

let stability = sc.cluster_stability(result, [0.2, 0.5, 0.8, 1.0], 15)
write_text("stability.svg",
    sc.plot_cluster_stability(stability, "Resolution stability"))

let diagnostics = sc.cluster_diagnostics(result)
write_text("silhouette.svg",
    sc.plot_silhouette(diagnostics, "Cluster separation"))

Run the complete example

The packaged example creates a deterministic multi-donor data set, executes every advanced analysis, and writes SVG and CSV results. It is the runnable source for the variables used in the excerpts above.

bl examples singlecell --copy singlecell-examples
cd singlecell-examples
bl run advanced_analysis.bl

One-Call Exploration

standard is useful for a first look. It prints the explicit pipeline and identifies which parameters are defaults. For a report or publication, keep the expanded pipeline and record every threshold.

let result = sc.standard(
    cells,
    resolution: 0.5,
    n_hvg: 2000,
    k: 15,
    min_genes: 200,
    max_genes: 5000,
    max_pct_mito: 20.0
)

Cross-Validation

The installed package includes seeded Scanpy and Seurat workflows under examples/validation. Compare cells by barcode and partitions with adjusted Rand index. Do not compare raw cluster numbers or PCA signs, because both can change without changing the biological partition. Copy the package examples into an independent working directory; no BioLang repository checkout is required.

bl examples singlecell --copy singlecell-examples
cd singlecell-examples
python make_demo_10x.py --output validation_data
bl run validation/biolang_reference.bl
python validation/scanpy_reference.py validation_data validation_scanpy_labels.csv
python validation/compare_labels.py validation_scanpy_labels.csv validation_biolang_labels.csv

# Cross-check the advanced paired test against SciPy and R.
python validation/advanced_reference.py
Rscript validation/advanced_reference.R

See the Single-Cell RNA-seq with BioLang book for interpretation, experimental design, batch effects, reproducibility, and the complete API reference.