Malva platform: unified single-cell sequence search
Malva is a scalable single-cell and spatial transcriptomics search platform that provides a unified interface for searching millions of cells by RNA sequence. The platform combines five independent, modular components connected through a central orchestrator (Extended Data Fig. 1b). Users can submit queries through the public RESTful API, including gene symbols, DNA or RNA sequences, and natural-language requests such as “cells expressing markers of neurodegeneration.”
The central orchestrator routes each request to the appropriate service. Natural-language queries are first processed by an integrated language model, which converts biological concepts into relevant gene sets. Sequence queries are decomposed into k-mers and searched against the distributed Malva Index. Frequently requested genes are accelerated through pre-computed results. Malva Index returns cell identifiers together with expression values or sequence-coverage profiles. A metadata service then adds standardized cell-level annotations, such as cell type, and sample-level annotations, including tissue, disease and developmental stage. Users can therefore filter search results or perform targeted analyses within defined sample subsets.
Malva can continuously expand its searchable corpus. A crawler monitors public repositories, including GEO, SRA, the Human Cell Atlas (HCA) Data Portal and ENA, for newly available datasets. These datasets are processed asynchronously into independent Malva Index instances and progressively merged into the main corpus. Load balancing across multiple workers supports reliable performance when hundreds of users submit queries simultaneously. Results are returned in standard formats, including JSON and h5ad, and contain complete metadata for downstream analysis with tools such as scanpy. The server uses an open-source software stack comprising Flask v3.1.0, numpy v2.2.6, pandas v2.3.2, SQLite v3.49.1 and polars v1.27.1.
Malva data sources, dataset crawling and availability
Malva Index incorporates public single-cell RNA-sequencing and spatial transcriptomics datasets from major repositories, including the Human Cell Atlas, GEO, SRA, ENA, ArrayExpress and CNGBdb. Dataset ingestion is managed by the malvacrawl service. This service maintains a registry of candidate datasets identified primarily through the HCA Data Portal, supplemented by automated parsing of GEO and ENA metadata feeds and manual ingestion of selected accessions from other resources.
Each candidate dataset receives an internal UUID and is checked against existing records using GEO or ENA accession numbers to prevent duplicate ingestion. Because public metadata are often incomplete or inconsistent—for example, sequencing technology and barcode format may be missing—malvacrawl applies a two-stage annotation workflow. An instruction-tuned Llama 3.1 language model first parses dataset descriptions and extracts candidate labels for technology, organism and tissue. These labels are then manually reviewed and mapped to controlled ontology terms with text2term v4.5.0 (see ‘Metadata standardization’).
For raw data acquisition, malvacrawl attempts to download FASTQ files using SRA prefetch, HTTP or FTP, in that order of preference, with sra-tools v3.4.1. Archived BAM files are converted to FASTQ using bamtofastq v1.4.1. Barcode and sequence parsing during preprocessing use dnaio v1.2.4. Barcoded technologies, particularly 10x Genomics assays, generally require paired-end sequencing. If only one FASTQ file is available, or if the dataset is archived only as BAM files, the BAM files are downloaded and converted using standard 10x conversion utilities.
Technology type is determined from deposited metadata, the structure of the R1 read—including barcode length and whitelist match rate—or manual inspection when the technology remains ambiguous. Datasets that are incompatible with Malva’s whitelist-based barcode system (see ‘Barcode preprocessing’) are deferred for later integration. For plate-based, non-barcoded assays such as Smart-seq2, each FASTQ file represents a single cell. These datasets are ingested in bulk mode as per-cell indices and later combined into sample-level indices while preserving cellular resolution.
Malva Index currently spans approximately 60 million cells from thousands of experiments. It includes diverse technologies such as 10x Chromium, Visium, Open-ST, Stereo-seq and Smart-seq2, covering multiple organs, species and biological conditions. The platform operates on on-premise infrastructure and applies query limits to maintain consistent performance for all users.
Querying genes, sequences and biological concepts
Malva maintains local databases that translate user requests into searchable nucleotide sequences. Gene symbols and Ensembl identifiers are mapped to genomic coordinates on standard reference genomes, including hg38 for human and mm10 for mouse, using default repeat masking55. For expression searches, Malva retrieves complete transcript sequences from Ensembl. Coverage queries that display read distribution along a gene use genomic coordinates to define the target region and report k-mer matches at each position.
Specialized databases support semantic biological searches. Malva integrates Llama 3.1 through Ollama as a natural-language translation layer between user prompts and gene sets56. The model can optionally use a Retrieval-Augmented Generation (RAG) database constructed from MSigDB57, which contains curated gene sets from Reactome, Gene Ontology (GO) and other pathway resources.
For a request such as “inflammatory response genes in macrophages,” the system retrieves relevant gene sets from MSigDB and uses the language model to refine them according to the query context. This allows researchers to explore biological processes without manually compiling gene lists. The resulting gene sets are converted into nucleotide sequences and searched across Malva Index using the same workflow as any other sequence query.
Pre-computed gene expression matrices for faster searches
To accelerate frequently used gene-level searches, Malva pre-computes and stores expression matrices for commonly queried gene sets. For each species represented in Malva Index, comprehensive transcript catalogues are compiled from Ensembl. These catalogues include protein-coding transcripts, non-coding RNAs and 3′ UTR sequences. The reference sequences are queried against Malva Index to generate pseudocount matrices for every indexed dataset.
The pre-computation workflow aggregates k-mer counts across all isoforms of each gene, producing a cell-by-gene matrix comparable to a conventional single-cell expression matrix. These matrices are stored in a distributed database indexed by gene symbol, Ensembl ID and dataset identifier. Common gene queries can therefore be answered directly from pre-computed values, reducing response times from seconds to milliseconds.
The database is updated incrementally as new datasets are added to Malva Index, ensuring that gene-level results reflect the complete searchable corpus. For custom or non-standard gene symbols, Malva falls back to real-time k-mer searching, preserving flexibility while optimizing performance for common queries.
Metadata standardization across single-cell studies
Cross-study aggregation requires consistent metadata. Malva harmonizes heterogeneous annotations from the HCA Data Portal, GEO, SRA, ArrayExpress and CNGB. Source metadata are supplied in different formats and use inconsistent terminology; the same tissue, disease, developmental stage or experimental condition may be described using multiple names.
Malva uses an automated metadata standardization pipeline based on text2term v4.5.058. This biomedical ontology-mapping tool uses semantic similarity algorithms to connect free-text terms with established ontology identifiers. The pipeline first normalizes null values and common terminology variants. Multi-value fields containing several biological conditions are separated and mapped individually to maximize annotation coverage.
Terms are then matched using TF–IDF (term frequency–inverse document frequency) semantic matching and hierarchical fallback strategies. Diseases are mapped to MONDO, anatomical structures to UBERON, phenotypic traits to PATO and developmental stages to age-related ontologies. Organs are additionally connected to hierarchical parent systems by traversing UBERON relationships. This multi-level harmonization maps biologically equivalent terms to consistent identifiers, enabling reliable cross-study aggregation and systematic biological searches while retaining both detailed and hierarchical relationships.
Automated cell-type classification
Malva includes automated cell-type annotations to support aggregation and interpretation of search results. Author-provided labels are frequently inconsistent between studies, and many datasets lack cell-type annotations entirely. Malva addresses this limitation with a two-stage de novo annotation pipeline.
First, a marker-based scoring system assigns preliminary labels. Curated marker gene sets for major cell types were compiled according to Cell Ontology classifications, with each cell type represented by 10–20 literature-supported markers. Cells are clustered using principal component analysis (PCA), k-nearest neighbours and the Leiden algorithm implemented through leidenalg v0.10.2. Differentially expressed genes are identified for each cluster, and enrichment scores are calculated for each cell-type signature. Cluster assignments consider both the number of expressed markers and their specificity.
Second, a language model refines the cluster-level annotations using biological context, as previously described59. The model receives the top differentially expressed genes together with the preliminary cell-type annotation, tissue of origin, developmental stage and disease state. It generates refined cell-type labels and evaluates whether each label is biologically plausible for the associated tissue context. Implausible assignments are flagged for correction, and final annotations are mapped to a controlled vocabulary based on the Cell Ontology.
Validation of the natural-language query interface
We evaluated the accuracy of Malva’s natural-language query translation using 10,000 synthetic test cases generated with gpt-oss:120b. This test-case generator was independent of the production model, Llama 3.1 8B. Queries were distributed across eight categories: simple gene queries (15%), gene queries with filters (20%), marker-gene requests (15%), sequence searches (10%), database identifier lookups (10%), pathway queries (10%), unsupported queries (7%) and over-inference traps (15%).
Over-inference traps tested whether the model incorrectly introduced biological filters that were not explicitly requested—for example, adding a “T cell” filter to an FOXP3 query when the user had not specified a cell type. We measured query-type accuracy, gene-list F1 score, filter precision and filter recall. The complete set of evaluation queries, expected outputs and prompts used for production-model instruction tuning are publicly available (Data availability).
In production, categorical metadata filters—such as tissue or disease terms—are resolved through a hybrid matching system. The input is first compared with a fuzzy dictionary containing all unique metadata values, using Levenshtein distance and a configurable threshold. If no confident match is found, Malva performs vector similarity searching against a sentence-embedding database built with the same model used for ontology mapping. The nearest ontology term is returned when its cosine similarity exceeds the defined threshold.
Malva Index architecture for scalable sequence search
Malva Index stores single-cell sequencing data as searchable collections of k-mers—fixed-length nucleotide sequences—linked to their cellular origins. Reads are decomposed into non-overlapping k-mers of length k, with a default value of k = 24. A final overlapping k-mer is added when necessary to ensure complete sequence coverage. For a sequence of length l, k-mers are extracted at positions 0, k, 2k, and so on, together with one k-mer beginning at position l − k when required.
Non-overlapping k-mer sampling, rather than minimizer-based sampling, was selected deliberately. Although minimizers guarantee shared k-mers between overlapping sequences, Malva’s windowed query strategy provides full sensitivity for queries containing at least 2k nucleotides (see ‘Sequence querying and pseudoquantification’). Non-overlapping sampling reduces indexing time and memory requirements, provides uniform detection probability and simplifies index merging through a union operation.
Each k-mer is encoded as a 64-bit integer using two bits per nucleotide: A = 00, C = 01, G = 10 and T = 11. Cellular origin is represented by a composite identifier encoding both the cell barcode and dataset source. Upper bits identify the dataset and lower bits encode the cell barcode. For example, allocating 8 bits to the dataset ID and 24 bits to the cell barcode supports up to 256 datasets with approximately 16 million cells per dataset. This hierarchical identifier allows every query result to be linked to both a specific cell and its source dataset.
The index uses a two-tier structure conceptually similar to an inverted web-search index. A sequence array contains all unique k-mers identified across datasets in numerical order. A location array stores the composite cell-dataset identifiers in which each k-mer occurs. A pointer array connects the two by recording the start and end positions of each k-mer within the location array. Finding all cells containing a particular k-mer therefore requires one binary search in the sequence array followed by a read from the corresponding location range.
Datasets are processed independently before being merged into the unified index. FASTQ files are decompressed and parsed in parallel, extracting barcodes and sequences at rates exceeding 2 million reads per core per second. Reads are streamed through memory-efficient buffers, and the resulting k-mers are lexicographically sorted and deduplicated while retaining the cells in which they occur.
Dataset-specific indices are merged hierarchically. Smaller indices are combined into progressively larger indices using k-way merge algorithms that require limited memory. This design supports processing datasets that exceed available RAM and provides near-linear scaling with data volume. Detailed descriptions of index chunk generation, chunk merging and query scoring—including pseudocode suitable for reimplementation—are provided in Supplementary Methods.
For indices larger than available memory, k-mers are divided into disk-resident buckets. A lookup table kept in memory maps each k-mer prefix to its bucket offset and size. A lookup requires one prefix-table access followed by a binary search within the bucket, typically involving fewer than 20 comparisons for billion-scale indices. Prefix tables and suffix arrays are accessed through memory-mapped file input–output, while cell blocks are loaded on demand.
The final index is stored in a compressed binary format designed for disk-based retrieval. Arrays are divided into fixed blocks of 512 entries and independently compressed using delta encoding followed by Blosc compression. This produces a 10–20× reduction compared with raw FASTQ files (Supplementary Methods). On an index containing N k-mers organized into pages of P k-mers, searches require approximately logpN page accesses—typically three to four random disk reads for billion-scale indices. A least recently used (LRU) cache retains frequently accessed pages in memory.
Barcode preprocessing and cell calling
Cell barcodes are extracted using technology-specific configurations and matched exactly against whitelist sequences, including 10x Chromium v2 and v3 barcode sets. All barcode–k-mer associations are retained during indexing. Barcodes with Hamming distance greater than zero from the whitelist are discarded, and error correction is not applied.
Technologies without vendor-supplied whitelists, including split-seq and sci-seq variants, require barcode resolution before indexing. Tools such as flexiplex44, matchbox60 and splitcode61 can convert non-standard barcode layouts into whitelist-compatible R1/R2 FASTQ pairs. Malva can then index these files using its standard workflow. Current Malva datasets were generated using protocols with fixed-length barcode subsequences, including 10x Chromium, Visium and Stereo-seq, or Smart-seq2, for which one index is generated per cell and later merged using the malva combine command.
Before cell-type annotation, Malva applies a UMI-based filtering method to retain barcodes likely to represent true cells. This method is analogous to the order-of-magnitude (OrdMag) algorithm used by Cell Ranger62. Barcodes are ranked by total UMI count, and the expected number of cells is estimated by minimizing the squared relative difference between a candidate cell count x and the number of barcodes exceeding a UMI threshold derived from x.
For each candidate x, the 99th percentile UMI count among the top x barcodes is divided by 10 to define a threshold. The number of barcodes above that threshold is represented by f(x). The optimal x minimizes (f(x) − x)2/x. The final UMI threshold is then calculated from the 99th percentile of the top barcodes divided by 10, and barcodes exceeding this threshold are retained as cells.
Sequence querying and k-mer pseudoquantification
Malva searches nucleotide sequences by decomposing each query into k-mers and looking for matches in Malva Index. Queries longer than k are evaluated with a sliding-window strategy so that all possible k-mer alignments are considered. For a query sequence q and window size w, set to 64 by default, Malva extracts windows of length w and calculates the fraction of matching k-mers in each window.
The match score for window W and cell C is defined as the number of matching k-mers divided by the total number of possible k-mers: \(\mathrm{score}(W,C)=\frac{|M(W,C)|}{w-k+1}\), where M(W,C) is the set of k-mers in W that occur in cell C. Low-complexity k-mers are masked by default with dustmasker v1.0.055. Users can optionally exclude repetitive k-mers that occur frequently across the reference transcriptome to reduce false positives.
A window is considered a positive match when its score exceeds threshold τ, which defaults to 0.65. For sequences longer than w, a cell is classified as positive when any individual window passes this threshold. Because a single nucleotide substitution affects only one k-mer, the method tolerates sequencing errors and mutations while preserving signal from the remaining matching k-mers.
Malva aggregates matches across all windows and ranks cells according to the number of significant windows. Under typical settings (w = 64, k = 24 and τ = 0.65), the method detects sequences with approximately 15% divergence while maintaining low false-positive rates. Query performance is improved by batching k-mer lookups, reducing disk access and exploiting the sorted index for cache-efficient retrieval. Full details of batching, index access and caching are provided in Supplementary Methods.
Sequence coverage quantification
Malva preserves transcript-level spatial information by reporting k-mer matches independently for each sliding window. A sliding average, with a default window of 24 bp, smooths the resulting coverage profile and produces visualizations comparable to alignment-based methods. Low-complexity k-mers are masked with dustmasker v1.0.055 to prevent artificial peaks. Without masking, shared repetitive sequences can create position-specific spikes of 100–1,000×. This positional coverage information supports the detection of alternative polyadenylation sites, splice-junction usage and other transcript-structure features.
Malva benchmarking and performance evaluation
We benchmarked Malva against established sequence-search and single-cell indexing tools, evaluating indexing speed, memory use and query accuracy. Comparisons included kallisto-bustools v0.28.222, MetaGraph v0.4.1 using the RowDiff/Multi-BRWT variant27, BLAST+ v2.16.0 using MEGABLAST63, a STAR v2.7.11b-based workflow23,64, Fulgor v4.0.0 using the meta-diff variant31 and REINDEER v1.4.7-1-g0812ad128. Sequence-search comparisons included MetaGraph, BLAST+, Fulgor and REINDEER.
Benchmark datasets varied according to each method’s architectural requirements. MetaGraph and BLAST were evaluated with the Stereo-seq CIRSTA mouse liver atlas from CNGBdb (CNP0003447), which contains approximately 60 billion reads, 7 TB of compressed raw data, 20 samples and approximately 4 billion unique spatial barcodes65. Because REINDEER and Fulgor require each cellular barcode to be supplied as a separate input file, processing millions of barcodes is impractical. REINDEER also requires bcalm66 to be run independently for each sample.
These tools were therefore benchmarked on a composite dataset containing six Visium samples: GSM5115019, GSM5115020, GSM5115021, GSM5115022, GSM5115023 and GSM5115024. The dataset included approximately 1.5 billion reads, 180 GB of compressed raw data and 30,000 unique spatial locations. It represented spatial resolutions ranging from 0.6 µm in Open-ST to 55 µm in Visium, as well as different sequencing depths and total data volumes.
For large-scale testing, we randomly selected 1% of the datasets in Malva Index, representing approximately 1.5 TB of raw data, and reindexed them with Malva and Fulgor. These methods showed the strongest runtime and memory performance during initial benchmarks. All analyses used identical settings. Each benchmark was run independently with a maximum allocation of 96 threads and 512 GB of memory on a multi-socket server containing Intel Xeon Platinum 8454H CPUs, 256 total threads and 4 TB of RAM. Resource use was monitored with snakemake v7.32.468.
Validation of biological signal preservation
We assessed whether Malva k-mer pseudocounts preserve biologically meaningful information comparable to conventional expression pipelines. Representative datasets included high- and low-resolution spatial transcriptomics, healthy and disease samples, and both human and mouse single-cell data. These comprised an Open-ST human tumour section69 (SRA accession SRX23008310), 10x v3 embryonic E13.5 mouse midbrain70 (SRA accession SRR11947658), 10x v3 human PBMCs downloaded from https://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v2/pbmc_1k_v3_fastqs.tar, and Visium mouse brain data downloaded from https://s3-us-west-2.amazonaws.com/10x.files/samples/spatial-exp/1.3.0/Visium_Adult_Mouse_Brain/Visium_Adult_Mouse_Brain_fastqs.tar.
For every dataset, Malva pseudocounts were compared with the deposited count matrices. We calculated gene-level and cell-level total counts and measured correlations between the two quantification approaches. To evaluate downstream performance, we clustered all benchmark datasets using species-specific Ensembl transcript sequences after masking low-complexity and repetitive k-mers.
Standard scanpy processing included normalization, log transformation, highly variable gene selection, PCA, nearest-neighbour graph construction and Leiden clustering with leidenalg v0.10.2. Cluster agreement between Malva and conventional expression matrices was assessed using Normalized Mutual Information. Reference clusters were also transferred to Malva objects by cell barcode. Marker genes were identified with Wilcoxon rank-sum tests comparing each cluster with all remaining cells, and concordance was measured using fold-change comparisons.
Trajectory analysis used scanpy v1.11.4 diffusion pseudotime (DPT) on the E13.5 mouse midbrain dataset71,72. Analysis was restricted to neuronal cells, with a representative radial glia cell selected as the root. Pseudotime was calculated independently from the Malva pseudocount matrix and the reference UMI-based matrix. Results were merged by cell barcode to quantify agreement between the inferred trajectories.
Simulation of sequencing errors
We tested Malva’s resistance to sequencing errors by introducing synthetic substitutions into Open-ST human tumour, 10x v3 human PBMC and Visium mouse brain datasets. Additional error rates ranged from 0.1% to 10% and were cumulative with errors already present in short-read sequencing data. At each rate, nucleotide positions were selected randomly within reads, and the original base was replaced with equal probability by one of the three alternative bases.
Malva indices were rebuilt from the altered reads, followed by the standard evaluation workflow. We measured per-cell and per-gene correlation with reference counts, clustering agreement using NMI and marker-gene concordance. For variant-detection experiments, we also measured sensitivity and specificity for cell-line mixture probes across error rates using both overlapping and non-overlapping k-mer indexing.
Validation of sequence-based variant detection
We evaluated Malva’s ability to detect sequence variants by comparison with Flexiplex v1.02.3, a specialized sequence-search tool44. The benchmark used a seven-cell-line mixture dataset (SRR10971813) containing known variants: a BCAS4–BCAS3 fusion in MCF7 cells, adenovirus E1A sequences in HEK293T cells and mitochondrial SNV rs878887783 in T47D cells.
Four target sequences were searched: the BCAS4–BCAS3 SNV variant (CCGATCCTGGGGCCGAGGTACCTTTGACAGGAGC), adenovirus 5 E1A (TTTGGACTTGAGCTGTAAACGCCCCAGGCCATAA) and a probe for rs878887783 (ATTACTCCTGCCATCATGACCCCTGGCCATAATATGATTTATCT). For the fusion variant, the reference sequence was excluded because MCF7 cells were not expected to express it. For the mitochondrial variant, we used a shortened probe because the full Flexiplex probe produced zero counts under the published supplementary protocol; the modified probe generated comparable detection numbers.
Flexiplex was run with default settings in grep-based detection mode. Cellular barcodes were recovered by matching Flexiplex read identifiers to the corresponding R1 FASTQ file and extracting the first 16 nucleotides. Malva used a 24-nucleotide sliding window and a 100% k-mer presence threshold. Malva directly returned cell barcodes and k-mer match counts. Wall-clock time included index construction for Malva, performed once, and query processing for each probe. True positives, false positives and false negatives were calculated by comparing detected cells with expected cell-line identities from the deposited annotated objects.
We further assessed clinically relevant somatic mutation detection using a lung adenocarcinoma dataset, PRJCA00173173, containing matched tumour and normal samples with characterized EGFR mutation status. Four tumours carried EGFR exon 19 deletions, and one carried the exon 21 L858R mutation (c.2573T>G). Ground truth was established by aligning reads to GRCh38 with STAR v2.7.11b and examining coverage across the EGFR locus.
Malva probes targeted exons 18–21 as positive controls for EGFR expression, the wild-type exon 21 sequence and a 45-bp probe centred on L858R (GTCAAGATCACAGATTTTGGGCGGGCCAAACTGCTGGGTGCGGAA). A matched wild-type probe was also used (GTCAAGATCACAGATTTTGGGCTGGCCAAACTGCTGGGTGCGGAA). Probes were searched against sample-specific Malva indices using w = 24 and τ = 1.0 for variant-specific detection.
Specialized sequence searches across Malva Index
Large-scale arbitrary nucleotide sequence search. To screen for exogenous sequences, we compiled viral genomes from NCBI RefSeq. The panel included SARS-CoV-2 (NC_045512.2), hepatitis C virus (NC_004102.1), measles virus (NC_001498.1), West Nile virus (NC_001563.2) and Zika virus (NC_012532.1) as negative controls; HERV-K consensus sequences; PhiX174 (NC_001422.1) as a positive control; and the lentiviral vector pBA439 (AddGene #85967).
These sequences were searched against the complete Malva Index using default settings of k = 24, a 64-nucleotide window and a score threshold of 0.65. Tissue-level and study-level coverage profiles were generated. We also screened for Mycoplasma, a common and frequently undetected cell-culture contaminant. The analysis included Mycoplasma hyorhinis, Mycoplasma hominis and Mycoplasma fermentans, which are frequently associated with human and mammalian cell cultures, as well as Mycoplasma agassizii, a turtle respiratory pathogen not typically associated with human or standard mammalian cell-culture workflows.
Germline variant probe design. To screen for germline variants, we generated probes targeting common SNPs in protein-coding regions. Phase 3 data from the 1000 Genomes Project were downloaded for human chromosomes 1–22. We retained biallelic SNPs with MAF ≥ 1% located in exons or UTRs of protein-coding genes defined by GENCODE v19. For each variant, genomic sequence surrounding the position was extracted from the GRCh37 reference genome, and paired 64-nucleotide probes were generated containing either the reference or alternative allele.
ClinVar classifications were used to categorize variant pathogenicity as pathogenic, likely pathogenic, variant of uncertain significance (VUS), likely benign, benign or unknown. Ancestral allele information from 1000 Genomes alignments was incorporated to determine whether the reference or alternative allele represented the ancestral or derived state. Probes with extreme GC content below 30% or above 70% were removed. The final collection contained 296,985 variant pairs across all chromosomes.
Variant probes were queried against the complete Malva Index through the standard API, returning cell-level match counts for reference and alternative alleles. Derived allele frequencies were calculated for evolutionary analyses, and observed frequencies in healthy and disease samples were compared with expected population frequencies. Ancestral-state annotations were used to test for systematic deviations that could indicate selection.
Isoform usage. Cell-type-specific isoform usage was evaluated with the Tabula Muris dataset, which provides full-length transcript coverage for exon-level analysis. Mouse genome annotations from GRCm39.113 were used to query individual genes against sample-specific Malva indices. Overlapping annotated exons were merged into consensus gene models, and mean exon coverage was calculated across cell types.
A nested marker-analysis strategy was applied. Marker genes were first identified through one-versus-all differential expression testing with the Wilcoxon rank-sum test and Benjamini–Hochberg correction, using log2 fold change ≥1.0 and adjusted P value ≤ 0.05. For marker genes, exon-level differential analysis identified cell-type-specific marker exons within broadly expressed genes. This enabled detection of cell-type-specific exon usage even when the corresponding host gene was expressed across several cell types.
We also examined 3′ UTR length regulation using the MOSTA Stereo-seq mouse embryonic development atlas, spanning E9.5–E15.5. The same reads were aligned with STAR v2.7.11b using default parameters to generate comparison coverage tracks. Malva coverage analysis focused on the Add2 locus and its proximal and distal polyadenylation sites, previously reported in bulk and single-cell RNA-sequencing studies. Proximal-to-distal 3′ UTR usage ratios were calculated and mapped spatially to identify tissue-specific isoform preferences during embryonic development.
Detection of circRNAs. To detect the circular RNA CDR1as, we designed a 36-nucleotide probe centred on the backsplice junction. The probe contained 18 nucleotides upstream of the splice donor site and 18 nucleotides downstream of the splice acceptor site, thereby targeting the junction that distinguishes circular from linear transcripts. Malva Index searches used w = 24 and τ = 0.65.
Benchmarking alternative polyadenylation detection against polyASite
Malva’s ability to detect alternative polyadenylation was evaluated against polyASite v3.051, a database of polyadenylation sites identified from Human Cell Atlas samples using 3′ end sequencing protocols. For each annotated polyadenylation site, we designed two probe types: polyA probes containing 44 bp of upstream genomic sequence followed by a 20-bp polyA tract (20*A), and control probes containing 64 bp of genomic sequence around the same position without the polyA tract.
Both probe sets were searched against the HCA samples used to construct polyASite. Spearman correlations between Malva pseudocounts and polyASite usage scores were calculated separately for each probe and tissue type. Because Malva pseudocounts were not depth-normalized, tissue-specific correlations were used to evaluate relative usage patterns rather than absolute expression levels. Probe sequences are publicly available (Data availability).
Large-scale somatic mutation detection and comparison with scTML
Malva’s mutation-detection performance was benchmarked at scale using the scTML database53, which contains uniformly processed single-cell RNA-sequencing data from 280 tumour samples spanning 16 cancer types, together with matched variant calls. The scTML workflow used GATK HaplotypeCaller to identify variants from aligned reads, followed by quality and read-support filtering.
We selected SNVs in exons or UTRs and excluded intronic variants because of their lower sequence complexity and greater repeat content. For each variant, a 45-bp probe centred on the variant position was generated and searched against the corresponding sample-specific Malva index. Cells were classified as variant-positive when they met the k-mer match criteria (w = 24, τ = 1.0). Sensitivity was calculated as the proportion of scTML-positive cells detected by Malva. Specificity was measured using matched normal cells or variant-negative tumour cells.
Reference-free single-cell analysis using k-mer composition
Conventional single-cell analysis depends on mapping reads to reference genomes and quantifying annotated genes. This approach can miss unannotated transcripts and is less suitable for non-model organisms. Malva provides a reference-free analysis workflow that clusters cells according to their complete k-mer composition, revealing sequence relationships that may not be visible in gene-based analyses.
Scalable k-mer bucketing
Directly analysing complete k-mer repertoires is computationally challenging because even moderate datasets can contain billions of unique sequences. Malva addresses this problem with a bucketing strategy that groups related k-mers while retaining sequence-level information.
The bucketing algorithm uses MinHash to assign k-mers to a fixed number of groups according to their w-mer composition, where w < k. For each k-mer, all contiguous w-mers are extracted. Typically, w = 16 when k = 24. Two independent hash functions, MurmurHash3 and xxHash64, are calculated for every w-mer. The minimum hash value from each function is retained, and the two minima are combined into a composite signature. The bucket assignment is determined by taking this signature modulo the predefined number of buckets, which defaults to 100,000.
This method increases the likelihood that k-mers with similar w-mer composition are assigned to the same or nearby buckets. The resulting cell-by-bucket matrix reduces the feature space from billions of unique k-mers to a manageable number of features while preserving sequence relationships. Assignments are deterministic and consistent across datasets and species. k-mers are processed in chunks, and bucket counts are accumulated per cell using sparse matrix operations, enabling analysis of datasets containing hundreds of millions of cells.
Assessing preservation of local cellular structure
We evaluated whether k-mer bucket representations preserve biologically meaningful neighbourhoods by comparing sequence-based embeddings with PCA embeddings generated from gene-count matrices. For every cell in the overlapping dataset, we identified k = 15 nearest neighbours in both representations.
Wasserstein distances were calculated between the distances from each cell to its nearest neighbours in sequence space and the distances from that cell in sequence space to neighbours identified in gene-count space. The resulting distribution was compared with a random null model using ROC-AUC analysis. Local Structure Agreement scores substantially above 0.5 indicate better-than-random preservation of local cellular relationships.
Cell clustering and discovery of marker sequences
We analysed the cell-by-bucket matrix using a standard scanpy workflow. Bucket counts were normalized to account for sequencing-depth differences, and highly variable buckets were selected according to dispersion relative to the mean. PCA was performed using 30 components, followed by nearest-neighbour graph construction and Leiden clustering with leidenalg v0.10.2.
Similarity between clusters generated from k-mer profiles and gene-count matrices was quantified with Normalized Mutual Information:
$$\mathrm{NMI}({Y}_{1},{Y}_{2})=\frac{2\times I({Y}_{1}; {Y}_{2})}{[H({Y}_{1})+H({Y}_{2})]}$$
Here, Y1 and Y2 represent the two clusterings, H is the entropy function and I(Y1; Y2) = H(Y1) − H(Y1|Y2) is the mutual information function.
To identify the biological features underlying each cluster, differentially abundant buckets were detected using a t-test with Benjamini–Hochberg correction. Features were selected using an adjusted P value < 0.05 and log2 fold change > 1. The k-mers assigned to significant buckets were then retrieved and assembled de novo with SPAdes v4.2.0 in RNA mode74.
Contigs were annotated in two stages. First, they were aligned to host references with minimap2 v2.1.1-r341, using spliced presets for genome mapping and ungapped mapping for transcriptome alignment against GRCh38. Transcript-level mapping used Ensembl cDNA and lncRNA catalogues. Contigs were classified as genic when they overlapped annotated transcripts, intergenic when they mapped uniquely outside annotated features or unmapped when no confident alignment was found.
For loci containing multiple contigs, per-cell Malva coverage profiles were compared across the locus to distinguish isoform-specific from locus-wide signal. Contigs without confident host assignments, including those with rRNA signatures, were analysed with Kraken2 v2.1.5 using a standard RefSeq PlusPF database containing bacterial, archaeal, viral, human and UniVec/adapter sequences. Expression was recalculated using both annotated and unannotated contigs, and cell-level totals and feature abundances were compared.
We also tested reference-free Malva analysis in an organism with limited genomic resources using publicly available Placozoa single-cell RNA-sequencing data75 from SRA project PRJNA982113. Reads were indexed with the standard Malva configuration of k = 24 and embedded using the sequence-bucketing workflow. Dimensionality reduction and clustering followed the scanpy procedure described above.
Reporting summary
Additional information about the research design is available in the Nature Portfolio Reporting Summary linked to this article.
Source: www.nature.com