MainRNA is of fundamental importance for life. Knowledge of RNA sequences, RNA expression, splicing, isoforms, structure and modifications are central for understanding and targeting cellular processes. In the past decade, new sequencing technologies have made possible RNA sequencing (RNA-seq) in single, dissociated cells and subsequently for cells within whole tissue slices. Since cells are the basic units of life, these technologies have transformed life sciences, including medical and clinical research and applications, and annually generate petabytes of sequence information from hundreds of millions of cells1,6,7,8,9. The Human Cell Atlas consortium systematically quantifies RNA in healthy human cells2,3,10, for instance, while researchers also apply these technologies to diseases11. Moreover, vast amounts of data are generated for model systems, including organoids, and will be applied to investigate biodiversity on Earth as a whole12,13.It is therefore of utmost importance to have tools that can both integrate these fragmented datasets into a unified resource and make them efficiently interrogable. Thirty-five years ago, BLAST pioneered this principle for sequence databases, transforming biological research by making it possible to locate a nucleotide sequence of interest in a vast collection of genomic sequences—critically, not only with accuracy but also with unprecedented efficiency14. These design principles made it broadly usable and ultimately indispensable across the life sciences. Yet, no existing computational tool—such as cell atlases (CELLxGENE, STOmicsDB, Single Cell Expression Atlas and others15,16,17,18,19,20,21)—can empower researchers to efficiently search for specific RNA sequences across the recently created vast single-cell, single-nuclei or spatial transcriptomics (SC/ST) data. At best, these platforms provide pre-computed gene counts that reflect the aggregated number of RNA sequences mapped to a reference. Therefore, critical information is missing: about transcripts that do not map to the reference (for example, transcripts generated by genetic lesions in cancer cells or transcripts that do not map to the reference haplotypes), RNA sequence editing, mutations, pathogens present in the cells and more.Standard single-cell and single-nuclei pipelines remain reference-centric: they process sequencing reads by aligning them to a user-provided reference22,23,24,25,26. Applying these approaches at atlas scale requires downloading and reprocessing petabytes of raw files, which is beyond the computational resources available to most laboratories. By contrast, bulk, sample-centric methods that index short sub-sequences (k-mers) directly from the data treat the dataset itself as the searchable unit rather than forcing every read through a reference. They enable ultra-rapid search across large sequencing collections (MetaGraph, REINDEER, Sequence Bloom Tree, Mantis, Fulgor, Themisto and Logan, among others27,28,29,30,31,32,33). These approaches have already proved to be very useful: several studies have discovered new plastic-degrading enzymes33, novel RNAs that shape tumour identity34,35,36, sequence variation across large cancer cohorts37, disease-associated sequences invisible to alignment-based detection38, functionally related long non-coding RNAs (lncRNAs) that share k-mer patterns despite lacking sequence homology39, and even new forms of life40. However, these tools were designed for collections of bulk samples (hundreds to hundreds of thousands of genomes or experiments). Their construction typically requires a compacted de Bruijn graph built from the union of all input data, an operation that must be at least partially repeated when new data are added. Annotation is stored either as a label-wide matrix or as inverted lists whose size increases with total label count; in both cases, efficient queries generally require the index to reside in memory, which at single-cell scale (tens to hundreds of millions of individual cells) would demand terabytes of RAM. These properties make existing tools impractical for this regime.Here we introduce Malva (https://malva.mdc-berlin.de), a computational platform that unifies SC/ST data and makes them efficiently searchable at nucleotide resolution. Malva comprises a sequence search algorithm, a continuously expanding searchable corpus (Malva Index) and a public application programming interface (API); just as the success of BLAST rested on speed and accuracy, Malva makes sequence search and analyses practical at atlas scale.Malva search can start from nucleotide sequences, gene identifiers or natural language prompts. This enables analyses on demand across millions of cells in seconds, from viral genomes and germline variants to isoform usage. Beyond search, Malva can perform reference-free cell-type clustering41,42,43 and de novo assembly of cell-type-specific sequences.Malva does not replace reference-based quantification pipelines, which remain suitable for gene-level analysis; rather, it extends them by enabling sequence-level queries that are otherwise inaccessible at very large scale, which is also critical in scenarios where references are not readily available13. To our knowledge, Malva is the first platform for large-scale reference-free sequence search and analysis at single-cell resolution.Fast sequence search at single-cell resolutionEnabling sequence queries across all public single-cell data, beyond pre-computed gene counts, requires addressing three fundamental challenges. First, preserving nucleotide-level resolution across massive datasets means that traditional alignment pipelines, which discard sequence information, are insufficient. Second, comparing data across studies demands harmonization and standardized metadata that bridge different technologies and processing methods. Third, even with aggressive compression, these data are too large to download locally, requiring decentralized access via an API over the Internet.To address all three challenges, we developed Malva (Fig. 1). Malva is a platform that comprises a new search algorithm, a continuously expanding searchable corpus (Malva Index), and a public API. Malva Index is built on thousands of public datasets and is accessible via the API for both interactive and programmatic access.Fig. 1: Malva enables instantaneous sequence-based queries across single-cell atlases.Full size imageConventional single-cell portals lack sequence-level resolution and existing sequence search tools lack single-cell resolution. Malva integrates these tools, providing real-time results at atlas scale. Researchers can query Malva Index using raw nucleotide sequences, gene symbols or natural language descriptions. Malva resolves these inputs to their underlying sequences and returns all matching cells together with rich metadata at both the cell and sample level (for example, cell types, disease status, sex, age and study identifiers). Results can be fed into computational systems for automated analysis (for example, neural networks) or analysed interactively, for example as expression distributions across cell types (A), coverage profiles along the queried sequence queries (B) or spatial maps when coordinates are available (C).From a user’s perspective, Malva transforms how researchers explore transcriptomic data. Traditional workflows of a researcher interested in circular RNA (circRNA) expression, for instance, require downloading entire atlases, running specialized circRNA detection pipelines and aggregating results across samples—a process that requires days of computation and specialized expertise. With Malva, searching for a circRNA across millions of cells takes seconds. In general, users can instantly query along four dimensions: the feature that they seek (a raw sequence, gene or pathway), where to search (filters such as tissue, disease state or protocol), how to measure it (base-level coverage or aggregated pseudocounts) and scale (individual cells, whole libraries or aggregated groups). These capabilities are demonstrated through three complementary search modes.First, researchers can provide a nucleotide sequence (for example, ranging from a short 24-mer splice probe to a full viral genome) and immediately retrieve the cells that contain it. Second, they can specify a genomic interval and view live coverage tracks coloured by lineage or tissue. Third, they can type natural language queries, such as “cells expressing high levels of mitochondrial genes in tumor samples”, which an integrated agent translates to structured searches (Table 1 and Supplementary Fig. 1).Table 1 Example questions supported by Malva and what they returnFull size tableIndexing and query performanceMalva’s search algorithm borrows the k-mer indexing principle from bulk sequence search and adapts it to single-cell data: each read is decomposed into k-mers (sub-sequences of length k) stored together with the cell barcode of origin, enabling constant-time lookups without alignment (Extended Data Fig. 1 and Methods). To overcome the scalability limitations of existing k-mer indices at single-cell scale (Supplementary Note 1), Malva makes three design choices: non-overlapping k-mer sampling eliminates de Bruijn graph construction; per-k-mer sparse inverted lists store only the expressing cell identifiers, decoupling query cost from total index size; and independent per-sample indexing with incremental merging yields an on-disk layout that supports memory-efficient queries and continuous corpus expansion (Supplementary Methods).To begin, we assessed indexing performance on a Stereo-seq mouse liver atlas containing about 61 billion reads across 20 sections, each with more than 3 billion reads and hundreds of millions of spatial barcodes. Malva completed indexing in 24 central processing unit (CPU) hours, between 4 and 25 times faster than other preprocessing algorithms (that is, k-mer indexing and reference mapping), with a peak memory footprint of 8 GB, whereas all other methods were benchmarked without memory constraints and required substantially more resources (Extended Data Fig. 2a and Supplementary Table 1). Indexing speed saturates gzip decompression of FASTQ files, an unavoidable input–output bottleneck; memory usage scales linearly with the number of unique cell barcodes rather than with sequencing depth, and index sizes scale with sequencing depth per cell rather than cell count or total unique k-mers (Extended Data Fig. 2b,c and Methods).Next, we validated that query performance enabled real-time analysis. When searching for a sequence, Malva tiles it into k-mers and counts exact matches per cell (Supplementary Fig. 2c). These raw match totals, termed pseudocounts, emerge directly from k-mer matches without alignment (guidelines for parameter choice provided in Supplementary Tables 2 and 3 and Supplementary Note 2). Search latency depends only on disk access: 70 ms for a single k-mer, 0.9 s for a 1 kb transcript and approximately 1 min for 1,000 transcripts on a single CPU core, all computed on demand rather than from pre-computed gene matrices (Extended Data Fig. 2d). Other methods had worse time or memory performance, and traditional pipelines would have required read realignment for every new search query (Supplementary Note 1). We additionally selected a large collection of single-cell datasets (about 1.5 TB of raw data) and reindexed them with Malva and Fulgor, the two methods with the best runtime and memory usage; Malva achieved faster indexing with a smaller memory footprint while maintaining comparable query throughput from disk (Supplementary Table 1).Third, we validated quantification accuracy of Malva on datasets covering different SC/ST technologies, tissue types and disease states (Extended Data Fig. 3 and Methods). Across datasets, pseudocounts showed strong correlation with reference counts at both cellular (Extended Data Fig. 3a) and gene (Extended Data Fig. 3b) levels. We note that pseudocounts are not corrected for unique molecular identifiers (UMIs) and therefore correlate with read counts rather than molecule counts; users should consider these as semi-quantitative measures (Discussion). Systematic offset observed across single-cell counts can be corrected through library-specific normalization based on sequencing saturation (Extended Data Fig. 4 and Supplementary Note 3), analogously to approaches for bulk transcriptomics37. This correction addresses first-moment differences (that is, mean count inflation due to PCR amplification) but does not account for gene- or sequence-specific amplification biases.Having established quantitative concordance, we assessed whether Malva preserves biological signatures that distinguish cell types. Using the original cell-type annotations for each dataset, differential expression analysis showed that cell-type-specific enrichment across genes correlated strongly between Malva and reference methods across all datasets (Extended Data Fig. 3c). Moreover, clusters from pseudocounts preserved spatial organization in tissue sections (Extended Data Fig. 3d) and cell-type separation in single-cell data (Extended Data Fig. 3e). Additionally, developmental trajectory analysis on the embryonic day 13.5 (E13.5) mouse midbrain showed that pseudotime inference using pseudocounts closely recapitulated trajectories derived from reference counts (Extended Data Fig. 3f).We further benchmarked Malva’s sequence search accuracy by its ability to identify cell lines through variant detection. The benchmark contains a mixture of seven cell lines and was originally developed to test Flexiplex, a sequence search tool44. Malva indexed the more than 200 million reads 4 times faster than Flexiplex’s total search time, then resolved each query in milliseconds (Extended Data Fig. 5a). Across all target sequences, Malva maintained both high sensitivity (ability to detect true positives) and specificity (ability to avoid false positives) (Extended Data Fig. 5b–d). Because Malva indexes non-overlapping k-mers, sensitivity for short probes (length = k) depends on alignment with indexed positions, but the sliding window query algorithm guarantees detection for probes longer than k (Extended Data Fig. 5e,f and Supplementary Note 2). Systematic evaluation across k-mer sizes and probe lengths confirmed that our default parameters operate within the optimal detection regime (Supplementary Fig. 2); small deviations from perfect sensitivity reflect methodological differences with Flexiplex, which permits Hamming distance 1 mismatches that exact k-mer matching excludes.We additionally report that Malva accurately detected tumour driver mutations on a small cohort of lung cancer patients and matched normal samples (Extended Data Fig. 6), and it is robust to sequencing errors; we benchmarked this by simulating error rates from 0.1% to 10% across benchmark datasets. Sensitivity, specificity, pseudocount correlations and downstream analysis results (for example, cell-type clustering) remained stable up to error rates of 2% and degraded gradually up to very high (10%) error rates (Supplementary Fig. 3).In summary, Malva achieved faster indexing and querying times than existing methods while maintaining a low memory footprint, and showed high concordance with traditional pipelines for expression quantification.Capturing cellular diversity at scaleHaving validated technical performance, we constructed Malva Index by compiling more than 140 TB of public SC/ST data into a searchable resource, in about 70 h of computational wall time (around 9,700 h of CPU time). The build of the human Index at the time of analysis (v2025.04a) encodes more than 100 billion total unique 24-mers from around 51 million cells (592 studies, 7,966 samples) and an additional approximately 10 million cells from mouse samples.We achieved this by implementing a crawler that continuously collects FASTQ files mainly from the Human Cell Atlas (HCA) Data Portal that is also compatible with the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) and Gene Expression Omnibus (GEO), European Nucleotide Archive (ENA) and CNGBdb data portals, extracts barcodes or spatial coordinates across major protocols, indexes and merges each read into Malva Index, a single harmonized container (Fig. 2a and Methods). As of May 2025, Malva Index covered around 80% of the Human Cell Atlas, or around 20% of human, single-cell runs from the SRA. Malva Index has less than 10% the storage footprint of source files, a critical reduction given the continuous expansion of these publicly archived data (Fig. 2b). Malva also maintains individual sample-level indices, enabling local queries when global access is unnecessary. Cell barcodes are matched exactly against technology-specific whitelists and filtered to ‘true’ cells via an automated approach (Methods); barcode error correction and doublet removal are omitted to maintain high processing throughput at atlas scale (Discussion).Fig. 2: Malva Index.Full size imagea, Public sequence archives (NCBI GEO, SRA, ENA, CNGBdb and the HCA data portal) are continuously crawled for single-cell RNA-seq (scRNA-seq), single-nucleus RNA-seq (snRNA-seq) and spatial RNA-seq datasets, which are indexed and integrated into Malva Index. b, Cumulative expansion of the human Malva Index (v2025.04a at the time of analysis, spanning around 60 million cells), with sequencing volume stratified by organ system and technology type. c, Disease representation in Malva Index, shown as a word cloud in which word size and position in the vertical axis reflect the number of indexed cells per condition, with words at the top corresponding to the largest number. Autoimmune, neurodegenerative and oncological disorders are among the most prevalent. COPD, chronic obstructive pulmonary disease; SCC, squamous cell carcinoma; SLE, systemic lupus erythematosus. d, Key metrics and metadata (part of the query results output) of the human Malva Index, organized by anatomical system and organ. Panels show cell numbers, total and unique annotated cell types, donor age distribution, sex balance, unique k-mer counts and technology contributions. Bars distinguish healthy (blue) and disease (red) samples.Malva Index spans embryonic, adult and pathological human material from all major organ systems, with a diverse representation of cell types, donor ages, male–female ratio, disease status and underlying technology, represented in a curated database of sample- and cell-level annotations (Fig. 2c,d and Supplementary Fig. 4a–c).Notably, we validated that querying remains ultra-rapid on the complete Index. Despite containing tens of billions of k-mers across millions of cells, Malva Index resolves short sequence probes (such as 24-mers) in milliseconds, and returns full-length cDNAs in seconds and catalogues of thousands of transcripts in minutes (Extended Data Fig. 7a). This sub-linear scaling in query time ensures scalability and strengthens future proofness.Together, Malva Index enables interactive exploration across a continuously expanding, vast corpus of SC/ST data. Users can interrogate Malva Index via the public API, which is freely accessible with ORCID accounts linked to academic users.Enabling queries beyond current portalsCurrent atlas-scale single-cell portals excel at interrogating counts from predefined reference genes. Thus, interrogation of unexpected (for example, bacterial, viral, fungal or synthetic) sequences, mutations, isoforms or regulatory variants is not possible without downloading and re-analysing petabytes of data. Malva enables these via instantaneous sequence search (Fig. 3 and Extended Data Fig. 7).Fig. 3: Malva enables sequence-level queries beyond current data portals.Full size imagea, Sequence search across Malva Index detects viral, retroviral and contaminant sequences. Columns show genomes or probes (for example, HERV-K, PhiX spike-in, lentiviral vector (pLentiv.), Mycoplasma species, SARS-CoV-2, HCV1 (hepatitis C virus probe), measles (ME), West Nile virus (WNV) and Zika virus (ZV)), rows represent human tissues. Dot size indicates the fraction of positive cells; colour denotes normalized expression. Mycoplasma species expected in human cell culture (Myco+) and negative-control species (Myco–) are included. hERVs, human endogenous retroviruses. b, Malva recapitulates the known distribution of human germline genotypes at single-cell resolution. Systematic screening of 296,985 common SNPs (64-bp probes) from the 1000 Genomes Project. Filtering pipeline from all variants (approximately 40 million) to the selected set (approximately 300,000 coding, biallelic SNPs with MAF ≥ 1%). Scatter plots (mean ± s.e.) show correlation between Malva-observed allele frequencies and population expectations in healthy and disease samples; random mutation controls are generated by altering a non-SNP position (Methods). Inset, distribution of deviations from the expected alternative allele frequency when subset to derived alleles, confirming no systematic bias in health or disease. c, Malva recapitulates isoform switches that remain hidden in conventional gene expression analysis. Exon usage at the Ptprc locus in Tabula Muris4. Coverage-like queries reveal cell-type-specific inclusion of exon A (encoding the CD45RA isoform in B cells) versus exclusion (encoding CD45RO-like isoform in other immune cells). Rows represent cell types, with dot size and colour reflecting exon inclusion. DC, dendritic cell; HC, haematopoietic cell; HSC, haematopoietic stem cell; NK, natural killer. d, Malva analysis explains 3′ UTR usage patterns in its spatial context—that is, 3′ UTR length regulation of Add2 during spatiotemporal mouse embryogenesis. Left, coverage profiles across stages E9.5–E15.5 show proximal (cyan) and distal (magenta) polyadenylation site (PAS) usage; Malva-derived profiles (top) closely match STAR alignments (bottom). Right: Stereo-seq spatial maps (MOSTA atlas50) highlight distal UTR usage in developing neural tissue versus proximal usage in liver. Chr., chromosome.Screening for arbitrary sequences across Malva IndexUpon providing a nucleotide sequence, Malva returns the cells that express it, regardless of whether it is part of a predefined reference (for example, human genome) or not. To demonstrate this, we compiled a small panel of viral genomes, common laboratory contaminants, endogenous retroviruses and experimental vectors, and interrogated Malva Index to identify cells and samples containing these (Fig. 3a).First, we assessed sensitivity via positive controls (Fig. 3a, left): for instance, PhiX sequencing spike-in, routinely added during library preparation, was detected across most datasets. Similarly, signal for Mycoplasma (a well-documented contaminant in cell culture) was detected in 0.5–2% of samples, mostly at the 23S ribosomal RNA (rRNA) locus (Extended Data Fig. 7c). This is consistent with prior reports of Mycoplasma contamination in RNA-seq datasets45. Similarly, human endogenous retroviral signal (specifically, human endogenous retrovirus K (HERV-K)) was detected, as expected, in embryonic samples where these elements are known to be active. Of note, we observed tissue-specific expression in adult samples—such as healthy kidney and heart—cell lines and diseases including cancer and neurological conditions, supporting reports of HERV-K reactivation in adult tissues and pathological states46,47. We also tested the retrieval of cells carrying specific synthetic constructs. In particular, we queried a lentiviral vector that is used for clonal tracing in a rhabdomyosarcoma study48. Top hits returned by Malva, across the whole vector body (Extended Data Fig. 7d), correspond to cells from that study.Then, we assessed specificity via negative controls (Fig. 3a, right): we did not detect signals for pathogenic viruses that are not represented in the current Malva Index (for example, SARS-CoV-2, hepatitis C, measles, West Nile and Zika viruses). However, false negatives remain possible if they are below the detection threshold. We only observed hits of SARS-CoV-2-positive cells from lung samples that were labelled as healthy. Notably, these correspond to a cell line experiment (derived from healthy donors) in which cells were infected in vitro (Extended Data Fig. 7b), highlighting how Malva overcomes annotation biases.Screening for germline variation across single-cell atlasesWe next used Malva Index to locate cells with specific mutation signatures. We designed sequence probes targeting 296,985 common germline single nucleotide polymorphisms (SNPs) from the 1000 Genomes Project, focusing on biallelic variants in protein-coding regions with minor allele frequency (MAF) ≥ 1% (Fig. 3b, top and Methods); this ensured that we focused on variants that reached equilibrium rather than rare, deleterious mutations.Population-based allele frequencies correlated strongly to Malva-observed frequencies across healthy and disease samples (Fig. 3b). Thus, single-cell RNA data in Malva Index recapitulated the underlying population genetic structure despite heterogeneous sequencing technologies and depths. We additionally observed that common coding variants in our panel behave neutrally regardless of disease status, as there was no deviation between observed and expected derived allele frequencies (Fig. 3b, inset). Finally, we validated that allele frequencies at the single-cell level remained around 50% irrespective of population frequency (Extended Data Fig. 7e), supporting how most common coding SNPs segregate at intermediate frequencies in outbred populations, and that individuals are frequently heterozygous.Together, variant querying recapitulated population genetic structure and aligned with expectations that variants reaching ≥1% frequency in coding regions represent evolutionarily stable, predominantly neutral alleles at a single-cell level.Quantifying isoform usage without mappingGene-centric atlases collapse many isoforms into single counts, obscuring cell-type-specific regulatory switches. We leveraged Tabula Muris4 (profiled at full-length coverage), to systematically compare expression across gene loci as a proxy for isoform usage. We identified 1,138 loci with cell-type-specific differential coverage patterns (Extended Data Fig. 7f): for example, Ptprc (also known as Cd45) in B cells predominantly included exon A (encoding CD45RA isoform) while other immune cells mostly excluded it (instead encoding CD45RO), recapitulating this well-known isoform variation across lymphocyte subsets (Fig. 3c). This pattern emerging solely from coverage-like analysis supports the ability of Malva to interrogate exon and isoform usage.Moreover, we focused on isoforms at the 3′ untranslated region (UTR): these are fundamental for regulation of transcript and protein localization and stability and, more importantly, the SC/ST data in the Malva Index are mostly 3′-biased, therefore suitable for exploring this landscape. As an example, we interrogated the 3′ UTR of Add2, which is known to have tissue- and stage-specific patterns in the mouse embryo49, against a high-resolution spatial atlas of mouse organogenesis50. Malva validated how distal 3′ UTR has the strongest signal in the developing neural tissue, whereas the proximal form accumulates at the developing liver. Malva’s coverage-like profiles closely matched those generated by mapping and alignment with a standard tool (Fig. 3d), confirming that the bulk-level lengthening of Add2 across developmental stages is likely to represent a shift in cell-type proportions (neuronal versus liver progenitors). More generally, systematic benchmarking against polyASite v3.051 confirmed that polyA-specific probes correlate more strongly with database usage scores than control probes, validating probe design for site-specific detection (Extended Data Fig. 8).Querying circRNAs across cell atlasescircRNAs, produced by back-splicing, remain invisible in gene expression matrices provided by most current atlases. Researchers currently rely on specialized circRNA databases that cover limited samples and often do not cover single-cell or spatial data. Malva enables interactive analysis of circRNA expression by searching for probes specific to their backsplice junctions. For instance, we searched for CDR1as, a conserved circRNA that is highly expressed in the brain52. Malva confirmed the expected cell-type specificity of CDR1as at single-cell resolution (Extended Data Fig. 7g): 95% of positive cells were excitatory neurons, with lower detection in inhibitory interneurons and very low detection across other tissues. This specificity, previously known from bulk studies, emerges from droplet-based data that were not specifically designed to capture circRNAs.Querying cancer somatic mutations at scaleWe next tested whether Malva could detect somatic mutations directly from single-cell RNA-seq data at pan-cancer scale. We used the scTML database53, which provides matched mutation calls across 280 tumour samples from 16 cancer types (Extended Data Fig. 9a), and designed allele-specific probes centred on expressed SNVs in exonic and UTR regions (Methods). Malva recovered mutation-positive cells with high specificity and sensitivity across all cancer types (Fig. 4a), despite differences in tumour type, sample size, sequencing depth and mutation frequency. Sensitivity remained high even for variants detected in few cells (Fig. 4b and Extended Data Fig. 9b), and across variants at different features types (Fig. 4c). Thus, this benchmark shows that Malva can recover cancer-associated mutations from large single-cell cohorts without realignment or running specialized variant calling algorithms.Fig. 4: Malva detects somatic mutations across single-cell cancer datasets.Full size imagea, Sensitivity versus specificity across 16 cancer types (280 samples) compared to quantification in the scTML database53. Each point represents one study and cancer type. Points show median specificity and sensitivity for each study and cancer type; error bars indicate the interquartile range. CRC, colorectal cancer; HNSC, head and neck squamous cell carcinoma; KIRC, kidney renal clear cell carcinoma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; OV, ovarian cancer; PAAD, pancreatic adenocarcinoma; STAD, stomach adenocarcinoma; THCA, thyroid carcinoma; UVM, uveal melanoma. b, Sensitivity as a function of mutation frequency (number of cells carrying the variant), aggregated by cancer type. c, Sensitivity stratified by mutation genomic context. UTR3, UTR5, Exonic and Splicing labels indicate the functional annotation of the genomic region of a genetic variant relative to RefSeq genes. Labels joined with semicolons indicate that mutations fall within a region with overlapping features. 3′ UTR variants show the highest sensitivity, followed by exonic variants, with splicing-proximal variants showing moderate sensitivity due to junction-spanning probe constraints. Mean sensitivity across cancer types (approximately 0.94) is similar to mutation detection rates in bulk RNA-seq, where k-mer indices such as Transipedia (running on REINDEER) report recall values of 0.875–0.945 at high-precision thresholds37. Cancer types in b,c are represented by the same colours as in a. In a, PAAD is an outlier due to mouse-derived samples in the original database; therefore, PAAD is excluded in b,c.Together, these analyses demonstrate how, rather than limited to predefined gene sets, Malva allows researchers to interrogate millions of cells across diverse functional layers of RNA biology—expression, mutations and isoforms—via a single sequence search interface.From cell-by-gene to cell-by-sequenceSearching with Malva enables retrieval of cells containing specific biological signals. However, this does not enable cell–cell comparison, which is critical for single-cell clustering, trajectory inference and marker discovery. These operations are carried in a reduced (latent) space computed from cell-by-feature matrices, commonly, cell-by-gene matrices (Fig. 5a).Fig. 5: Sequence-based transcriptomic analysis without predefined gene annotations.Full size imagea, Conceptual view of the single-cell information funnel. Biological (‘physical’) cells are converted to ‘digital’ cells via sequencing, which can be represented either as cell-by-gene or cell-by-sequence matrices. These representations enable reconstruction of latent spaces underlying cell types and trajectories. PC, principal component. b, Latent space agreement comparing sequence-based (k-mer bucket) versus gene count principal component embeddings across healthy human cells. Values greater than 0.5 indicate better-than-random preservation of local neighbourhoods. n = 29,240,986 cells examined over 4,932 sequencing runs. Box plots show the median (centre line) and interquartile range (box edges), and whiskers extend to 1.5 times the interquartile range. AUC, area under the curve. c, Example comparisons of clustering based on gene features versus sequence features in human samples. Panels show cases with low (NMI = 0.31) and high (NMI = 0.65) concordance, corresponding to tissues highlighted in b. Bsl, basal epithelial cell; InfF., inflammatory fibroblast; Lum, luminal epithelial cell; Macro., marophages; Mes., mesenchymal progenitor cell; Msl., muscle fibre precursor; Myo., myofibroblast; MyoF., myofibroblast; N-ker., non-keratinized epithelial cell; Strat., stratified epithelial cell. d, Sequence-based clustering of non-model organisms (Hoilungia and Trichoplax). Clustering structure based on k-mer features closely resembles gene-based annotations. Epid, epithelial dorsal cells; Epigl, epithelial gland-like cells; Epiu, epithelial and unknown cluster cells; Epiv, ventral epithelial cells; Fibre, fibre cells; Gland, glandular cells; Lipo, lipophil cells; Meio, meiotic cells; Pept, peptidergic cells. e, Spatial clustering of a head and neck tumour dataset using gene count-based (left) versus sequence-based (Malva, right) embeddings. KS1 and KS2 are two k-mer-specific clusters. CAF, cancer-associated fibroblasts; Neutr. recr., neutrophil-recruiting population. f, Cluster-specific de novo assembly of marker sequences. Contigs assembled from discriminative k-mers reveal both annotated human features (genes and intergenic regions) and unaligned sequences. The dot plot shows the fraction of cells per cluster expressing each assembled contig (dot size) and the average log-normalized expression rescaled per contig so that the highest mean expression across clusters is set to 1 (dot colour). g, Representative marker sequences assembled from KS1 and KS2 clusters. Examples include a truncated 3′ UTR isoform of KLK10 (human) and bacterial 23S rRNA sequences from Gemella morbillorum. h, Spatial expression maps of selected marker sequences from g, shown as kernel density estimates. Distinct patterns confirm biological specificity: the truncated KLK10 isoform localizes to tumour regions, whereas microbial rRNA contigs map to extra-tissue areas, co-occurring with oral bacteria such as G. morbillorum and C. gingivalis.Malva Index encodes billions of sequences into a large sequence-by-cell matrix. By transposing it, we could in principle perform sequence-based, cell–cell comparison. In practice, these cell-by-sequence matrices are too large to handle (millions of cells times billions of features). Thus, we collapse sequence k-mers into buckets of sequence similarity (Methods). This ensures that cell-by-feature matrices are small enough such that dimensionality reduction is tractable.We tested this idea on approximately 30 million cells from around 5,000 samples from healthy individuals and found that the principal components derived from cell-by-bucket matrices were comparable to those derived from cell-by-gene matrices (Fig. 5b). We quantified this via two metrics: (1) neighbourhood preservation; and (2) clustering similarity (Methods). Across samples and tissue types, the degree of concordance between sequence-based and gene-based embeddings was variable (Fig. 5c and Extended Data Fig. 10a); in some tissues, cells that are neighbours in cell-by-gene space are also neighbours in cell-by-bucket space (high normalized mutual information (NMI)); in others, the two representations captured distinct structures (low NMI). This variation reflects the complementary nature of sequence-based analysis: although it can recover gene-driven neighbourhoods, it may also reveal additional structure. Therefore, these embeddings are an orthogonal layer of information to gene-based analysis.This approach becomes particularly useful in settings where high-quality reference genomes are lacking, sparse or incomplete; for example, for non-model organisms such as the phylum Placozoa (Fig. 5d), showcasing the utility of sequence-based representations for large-scale biodiversity studies13. We observed, however, that samples coming from single-nuclei technologies yield consistently lower values across both metrics. We attributed this to the over-representation of low-complexity intronic (nuclei-enriched) RNA (Extended Data Fig. 10b), which is equivalent to yields from low sequence diversity (for example, using 4-mers and below) (Extended Data Fig. 10c and Supplementary Note 4).Reference-free cell clusteringCell clusters are often characterized by their marker genes. However, in bucket-based clustering, features (groups of similar sequences) are not directly explainable. We assessed cell-by-bucket representations on a spatial transcriptomics sample (high-resolution profiling of a human tumour), in which tissue architecture provides an orthogonal readout for interpretability; that is, the degree (or lack of) of spatial organization can inform biological relevance. Indeed, gene-based and bucket-based representations captured similar major tissue compartments (tumour and stromal), plus two k-mer-specific clusters (KS1 and KS2) with distinct spatial organization, one within the tissue boundaries and another outside the tissue boundaries (Fig. 5e). We hypothesized that there might be sequence signatures that are specific to these two clusters, whose signal is absent from the original reference-based gene representation.De novo detection of marker sequencesAcross all bucket-based clusters, we extracted their specific collections of k-mers and used them as the input for de novo assembly (Methods). This yielded thousands of contiguous unique sequences (contigs), which we termed ‘marker sequences’ (Extended Data Fig. 10d). Upon mapping these to the human genome, around 70% of these contigs were annotated to known genes (mostly within 1 kb of the transcription end site (TES)), around 29% were annotated to intergenic regions, and approximately 1% remained unmapped (Extended Data Fig. 10d–f). Also, by narrowing the k-mer search space to those informative of cluster identity, downstream assembly (with SPAdes) operated more efficiently (requiring about 1 h versus about 27 h of wall time).We then examined the cluster-specific enrichment of all assembled contigs. All clusters (except KS1 and KS2) were enriched in contigs that mapped to human genes. Some contigs mapped to the same gene bodies but across different locations, and a subset of these were significantly detected as markers (Fig. 5f). One example is KLK10_676, a contig that maps to KLK10, a gene encoding a tumour-enriched serine protease. This was the only KLK10 contig that was enriched in keratin pearl structures (Fig. 5g,h). Coverage-like profiles showed that this feature may correspond to an unannotated 3′ UTR end; in fact, a polyA tract at its 3′ end further supports it as a TES. Conversely, KS1 and KS2 were indeed enriched in intergenic or unknown contigs (Fig. 5f). We first identified that KS2 was enriched in sequences that mapped to human rRNA (Fig. 5g,h), which are often removed during reference preprocessing, and were thus missing from the original cell-by-gene matrix. The KS1 cluster was also enriched in contigs that are compatible with bacterial 23S rRNA from several species that are known to occur in the oral mucosa (Fig. 5g,h). This confirms that this cluster separated from a genuine microbial signal that is only captured in the sequence-based analysis.Together, these results show that Malva can unbiasedly detect signals from SC/ST datasets in a reference-agnostic manner, providing insights that remain invisible to gene-centric analysis.DiscussionWe built Malva to make raw single-cell and spatial transcriptomic sequences searchable at (sub)cellular resolution. We indexed public datasets across millions of cells and show that sequences can be queried directly and instantly, without first reducing the data to genes, transcripts or predefined annotations. This makes it possible to search for and compare between splice junctions, mutations, viral sequences, synthetic constructs and non-reference RNAs across many studies, in seconds. This is particularly relevant for artificial intelligence-assisted analysis: neural networks learn from biological data but typically cannot return to the underlying experiments during inference. Malva provides a retrieval layer that allows models and agents to query observed sequences, retrieve cellular evidence and refine their predictions against public data.Many biological questions begin with such queries. A mutation can define a cancer clone, a splice junction can define a cell state, a viral RNA can define an infected cell, and a barcode can define a perturbation. These signals are not recoverable from processed count matrices. Malva keeps them accessible and links them to cells, tissues, disease states and metadata.Beyond individual searches, Malva can classify cells from sequence composition, assemble cell-type-specific marker sequences and may serve as a real-time, sequence data source for building artificial intelligence models. It therefore provides a practical discovery layer for asking sequence-defined questions in single-cell biology.LimitationsMalva currently supports queries of at least 24 nucleotides and uses exact k-mer matching, which enables fast searches of high specificity and sensitivity, with sensitivity being reduced only when the underlying data contain very high error rates, which are not typical in single-cell datasets sequenced with standard technology. Malva reports k-mer-derived pseudocounts rather than absolute molecule counts, and results should be followed by task-specific quality control, quantification or validation when needed. The current Malva Index also reflects the public data available today. Many datasets are short-read and 3′-biased, and spatial metadata are often incomplete or difficult to standardize. Finally, sequence search in human transcriptomic data can reveal donor-specific genotype information; this risk should be treated similarly to the risks associated with raw genome and exome sequencing data (Supplementary Discussion).ConclusionAs single-cell and spatial datasets expand, searchable sequence indexes may become a natural way to reuse public data. Future releases of Malva Index could include more organisms, spatial sequencing- and imaging-based assays, long-read transcriptomics, chromatin accessibility, proteomics and RNA modifications, among others. More generally, Malva can enable public single-cell archives to become live molecular resources, queried by sequence and linked back to cells, tissues and biological context. Ultimately, Malva may help to provide a coordinate system for mapping one’s own data onto global datasets to understand health and disease mechanisms and trajectories54.MethodsThe Malva platformMalva provides a unified interface to search millions of cells by any RNA sequence. Malva integrates five modular components that operate independently but communicate through a central orchestrator (Extended Data Fig. 1b). When a user submits a query via the public RESTful API—a gene symbol, DNA sequence, or natural language request such as “cells expressing markers of neurodegeneration”—it is routed by a central orchestrator. Natural language queries first pass through an integrated language model that translates biological concepts into gene sets. The search engine then decomposes sequences into k-mers and queries our distributed index, leveraging pre-computed results for common genes. Malva Index returns cell identifiers with expression levels or coverage profiles, which are then enriched by a metadata service that appends standardized cell-level (for example, cell-type annotation) and sample-level (for example, tissue, disease and developmental stage) annotations. This enables users to filter results or perform targeted searches on specific sample subsets. The platform can also expand datasets available in Malva Index: a crawler monitors public repositories (GEO, SRA, HCA Data Portal and ENA) for new datasets, which are asynchronously processed into separate Malva Index instances and progressively merged into the searchable corpus. Load balancing across multiple workers ensures consistent performance even as hundreds of users query simultaneously. Results return in standard formats (JSON, h5ad) ready for downstream analysis in familiar tools such as scanpy, with full metadata to enable further filtering and aggregation. The server code operates on a standard, open-source stack of tools (Flask v3.1.0, numpy v2.2.6, pandas v2.3.2, SQLite v3.49.1 and polars v1.27.1).Data sources, crawling and availabilityMalva Index incorporates public single-cell and spatial transcriptomic datasets from major repositories (Human Cell Atlas, GEO, SRA, ENA, ArrayExpress and CNGBdb). Data are ingested into Malva Index via the malvacrawl service. The service maintains a registry of datasets identified primarily from the Human Cell Atlas Data Portal, supplemented by automated parsing of GEO and ENA metadata feeds, and manual ingestion of selected identifiers from other platforms (such as ArrayExpress or CNGB). Each candidate dataset is assigned an internal UUID (universally unique identifier) and checked against existing entries by GEO or ENA accession to prevent duplication. Because deposited metadata are frequently incomplete or inconsistent (for example, the sequencing technology or barcode format is often not specified), malvacrawl utilizes a two-stage annotation strategy. First, a large language model (instruction-tuned Llama 3.1) performs an initial parse of the dataset description, extracting candidate technology, organism and tissue labels. Second, these candidates are manually curated and mapped to a controlled ontology using text2term v4.5.0 (see ‘Metadata standardization’). For raw data acquisition, malvacrawl attempts to download FASTQ files via SRA prefetch, HTTP, or FTP in order of preference, using sra-tools (v3.4.1). For datasets archived as BAM files, we converted them back to FASTQ using bamtofastq (v1.4.1). Barcode and sequence parsing during preprocessing relied on dnaio (v1.2.4). For barcoded technologies (primarily 10x genomics), paired-end data are expected; therefore, if only a single FASTQ file is available, or only BAM files are archived, these BAM files are downloaded and converted back to FASTQ using standard 10x conversion utilities. The technology is identified either from deposited metadata, from the structure of the R1 read (barcode length and whitelist match rate), or, when ambiguous, via manual inspection. Datasets that are not compatible with Malva’s whitelist-based barcode system (see ‘Barcode preprocessing’) are deferred for future integration. For non-barcoded, plate-based assays (for example, Smart-seq2), each FASTQ file comprises one cell, and datasets are ingested in ‘bulk’ mode as per-cell indices, and later merged into a sample-level index that preserves cellular resolution.Currently spanning ~60 million cells from thousands of experiments, Index includes diverse technologies (10x Chromium, Visium, Open-ST, Stereo-seq, Smart-seq2) across multiple organs, species, and conditions. The platform runs on-premise infrastructure with query limits to ensure consistent performance for all users.QueryingMalva maintains local databases to translate user queries into searchable sequences. When users query gene symbols or Ensembl IDs, these are mapped to genomic coordinates on standard reference genomes (hg38 for human, mm10 for mouse) with default repeat masking55. For expression queries, we retrieve full transcript sequences from Ensembl. For coverage queries that visualize read distribution along genes, we use the genomic coordinates to define the region of interest and return k-mer matches at each position.Several specialized databases support specific semantic biological queries. To support these via natural language, we integrated Llama 3.1 via Ollama as a translation layer between user prompts and gene sets56. The language model is optionally augmented with a Retrieval-Augmented Generation (RAG) database built from MSigDB57, which includes curated gene sets from Reactome, Gene Ontology (GO) and other pathway databases. When users submit queries such as “inflammatory response genes in macrophages”, the system retrieves relevant gene sets from MSigDB and uses the language model to refine the selection based on the specific query context. This approach provides a convenient interface for exploring biological concepts without requiring users to manually compile gene lists. The resulting gene sets are converted to sequences and searched across Malva Index like any other query.Pre-computed gene expression matricesTo accelerate common gene-level queries, we pre-compute and store gene expression matrices for frequently accessed gene sets. For each species in Malva Index, we compile comprehensive transcript catalogues from Ensembl, including all protein-coding transcripts, non-coding RNAs, and 3′ UTR sequences. We then query Malva Index with these reference sequences to generate pseudocount matrices for every indexed dataset.The pre-computation process aggregates k-mer counts across all isoforms of each gene, producing a cell-by-gene matrix analogous to standard single-cell expression matrices. These matrices are stored in a distributed database indexed by gene symbol, Ensembl ID, and dataset identifier. When users query common genes through the API, the system retrieves pre-computed values directly from the database rather than searching the k-mer index, reducing query latency from seconds to milliseconds. The database is updated incrementally as new datasets are added to Index, ensuring that gene-level queries always reflect the complete corpus of indexed data. For custom non-standard gene symbols, the system falls back to real-time k-mer queries, maintaining flexibility while optimizing performance for common use cases.Metadata standardizationCritical to the ability of Malva to aggregate results across studies is the harmonization of heterogeneous metadata from diverse sources (HCA Data Portal, GEO, SRA, ArrayExpress and CNGB). Raw metadata arrives in various formats with inconsistent naming conventions, where identical biological entities may be described using different terminologies across repositories. Similarly, developmental stages, disease states, and experimental descriptors vary widely between data sources.We implemented an automated standardization pipeline based on text2term v4.5.058, a specialized biomedical ontology mapping tool that utilizes semantic similarity algorithms to map free-text metadata terms to established ontology standards. The pipeline first preprocesses metadata by standardizing null values and normalizing common terminology variants. Multi-value fields containing multiple biological conditions are decomposed and mapped individually to ensure comprehensive coverage. The system then applies TF–IDF (term frequency–inverse document frequency)-based semantic matching with hierarchical fallback strategies to map terms to domain-appropriate ontologies: diseases to MONDO, anatomical structures to UBERON, phenotypic traits to PATO, and developmental stages to age-related ontologies. To enable system-level biological analysis, organs are additionally mapped to their hierarchical parent anatomical systems through ontology traversal within UBERON. This multi-layered semantic harmonization approach ensures that functionally equivalent terms from different repositories are mapped to identical standardized identifiers, enabling robust cross-study aggregation and systematic biological queries while preserving both granular and hierarchical biological relationships.Cell-type classificationTo enable aggregation of query results at the cell-type level, we implemented an automated annotation pipeline for all indexed cells. Author-provided cell-type labels are often heterogeneous across studies or entirely absent in many submissions. To address this, we perform de novo annotation using a two-stage approach.First, we apply a marker-based scoring system. We compiled curated marker gene sets for major cell types following Cell Ontology classifications, with each cell type defined by 10–20 literature-validated markers. We then perform clustering using standard methods (principal component analysis (PCA), k-nearest neighbours, Leiden algorithm via leidenalg v0.10.2) and identify differentially expressed genes for each cluster. Across clusters, we calculate enrichment scores for each type signature; thus, clusters are assigned cell types based on the differential expression of marker genes, considering both the number of expressed markers and their specificity to the cluster.Second, we refine these cluster-based annotations using a language model that integrates biological context, as previously described59. The model receives the top differentially expressed genes for each cluster along with metadata (marker-based cell-type annotation, tissue origin, developmental stage, disease state) and produces refined cell-type assignments. To ensure biological accuracy, the model evaluates whether each annotation is plausible given the tissue context, flagging unlikely assignments for correction. We mapped the output cell-type annotations to a predefined, controlled vocabulary based on the cell ontology standard.Natural language query interface validationTo systematically evaluate the accuracy of natural language query translation, we generated 10,000 synthetic test cases using gpt-oss:120b as a test case generator, ensuring independence from our production model (Llama 3.1 8B). Test cases 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 were specifically designed to test whether the model inappropriately adds biological context as filters (for example, adding ‘T cell’ filter when querying FOXP3 without explicit cell-type specification).We evaluated four metrics: query type accuracy (correct classification into the eight categories), gene list F1 (harmonic mean of precision and recall for extracted gene names), filter precision (fraction of returned filters that were explicitly stated in the query), and filter recall (fraction of explicitly stated filters that were captured). The complete set of 10,000 evaluation queries, together with expected outputs and the prompts used during instruction-tuning of the production model, are publicly available (Data availability).In the production environment, categorical metadata filtering (for example, resolving a user’s tissue or disease term to standardized ontology values) is performed via a hybrid approach. First, the input term is matched against a fuzzy dictionary built from all unique values in the metadata store, using Levenshtein distance with a configurable threshold. If no confident match is found, the system falls back to vector similarity search against a sentence-embedding database (built with the same model used for ontology mapping), returning the nearest ontology term above a cosine similarity threshold.Malva Index architectureMalva Index organizes single-cell sequencing data as a searchable collection of k-mers (fixed-length nucleotide sequences) linked to their cellular origins. Each sequencing read is decomposed into non-overlapping k-mers of length k (default k = 24), with a final overlapping k-mer to ensure complete sequence coverage. Specifically, for a sequence of length l, we extract k-mers at positions 0, k, 2k, ..., plus one k-mer starting at position l − k if needed. The use of non-overlapping k-mers rather than minimizer-based sampling is a purposeful design choice. We acknowledge that minimizers were developed to guarantee shared k-mers between overlapping sequences (essential for assembly); however, Malva’s windowed query strategy already ensures full sensitivity for queries with ≥2k nucleotides (see ‘Sequence querying and pseudoquantification’). Importantly, non-overlapping sampling greatly reduces indexing time and memory cost, provides uniform detection probability and simplifies index merging to a union operation.Each k-mer is encoded as a 64-bit integer using 2-bit nucleotide encoding (A = 00, C = 01, G = 10, T = 11), enabling efficient storage and comparison. The cellular origin is tracked through a composite identifier that encodes both the cell barcode and dataset source. We allocate the upper bits for dataset identification and lower bits for cell-specific barcodes; for example, using 8 bits for dataset ID and 24 bits for cell barcodes allows indexing up to 256 datasets with ~16 million cells each. This hierarchical encoding enables queries to return both the specific cell and its source dataset.Malva Index uses a two-tiered structure optimized for sequence search conceptually similar to inverted indices for web search. The sequence array stores all unique k-mers found across datasets, sorted numerically. The location array contains the composite cell-dataset identifiers where each k-mer appears. These are linked via the pointer array mapping each k-mer to its locations, storing start and end positions within the location array. This design enables rapid lookup: finding all cells containing a specific k-mer requires just one binary search in the sequence array followed by a single read from the location array.Individual datasets are processed independently before being merged into the unified index. FASTQ files are decompressed and parsed in parallel, extracting cell barcodes and sequences at rates exceeding 2 million reads per core per second. For each dataset, reads are streamed through memory-efficient buffers, extracting k-mers and their associated barcodes. The k-mers are lexicographically sorted and deduplicated, keeping track of all cells where each appears. These dataset-specific indices are then merged hierarchically: smaller indices combine into progressively larger ones using k-way merge algorithms that require minimal memory. This approach enables processing of datasets far larger than available RAM while maintaining near-linear scaling with data size. Detailed algorithmic descriptions of the index chunk generation, chunk merging and query scoring procedures, including pseudocode sufficient for reimplementation, are provided in Supplementary Methods.To handle indices exceeding available memory, k-mers are partitioned into buckets, each stored contiguously on disk (Supplementary Methods). A lookup structure, always resident in memory, maps k-mer prefixes directly to the disk offset and size of its bucket. Locating a k-mer requires one prefix table lookup followed by a binary search within the bucket, typically fewer than 20 comparisons for billion-scale indices. The prefix table and suffix arrays are accessed via memory-mapped file input–output; cell data blocks are read on demand for each queried bucket.The final merged index is stored in a compressed binary format optimized for disk-based retrieval. Arrays are partitioned into fixed blocks (512 entries) that are independently compressed using delta-encoding followed by Blosc compression, achieving 10–20× reduction compared to raw FASTQ files (Supplementary Methods). The hierarchical organization enables efficient searches: locating a k-mer on an index with N k-mers in pages of P k-mers requires in the order of logpN page accesses, typically 3–4 disk random reads for billion-scale indices. A small least recently used (LRU) cache maintains frequently accessed pages in memory.Barcode preprocessingCell barcodes are extracted according to technology-specific configurations and matched via exact lookup against whitelists (for example, 10x Chromium v2/v3 barcode sets). All barcode–k-mer associations are stored during indexing. Barcodes with Hamming distance >0 from the whitelist are discarded; error correction is not performed. For technologies that do not provide a vendor-supplied barcode whitelist, such as split-seq or sci-seq variants, barcodes must be resolved as a preprocessing step before indexing. Tools such as flexiplex44, matchbox60, or splitcode61 can convert non-standard barcode layouts into whitelist-compatible R1/R2 FASTQ pairs, after which Malva indexes the data in the standard way. The datasets currently in Malva Index were generated using protocols for which barcodes can be extracted as fixed-length sub-sequences (for example, 10x Chromium, Visium, Stereo-seq), or, in the case of Smart-seq2, where one index per cell is generated and subsequently merged using the malva combine command.Prior to cell-type annotation, we apply a UMI-based threshold method to retain barcodes that are likely to correspond to true cells, analogous to the order of magnitude (OrdMag) algorithm used by Cell Ranger62: barcodes are ranked by descending total UMI counts, and an 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. Specifically, for each candidate x, we compute the 99th percentile UMI count among the top x barcodes, divide this value by 10 to obtain a threshold, and count barcodes exceeding this threshold. The optimal x minimizes (f(x) − x)2/x, where f(x) is the count of barcodes above the threshold. Once the expected cell number is determined, the final UMI threshold is set as the 99th percentile of the top barcodes divided by 10, and all barcodes with UMI counts above this threshold are retained as cells.Sequence querying and pseudoquantificationMalva identifies query sequences by decomposing them into k-mers and searching Malva Index for matches. For queries longer than k, we utilize a sliding window approach to ensure all possible k-mer alignments are considered. Given a query sequence q and window size w (default w = 64), we extract all windows of length w and evaluate the fraction of matching k-mers in each window.The match score for a window W and cell C is defined as the ratio of matching k-mers to total possible k-mers: \(\mathrm{score}(W,C)=\frac{|M(W,C)|}{w-k+1}\), where M(W,C) represents the set of k-mers in W that appear in cell C. By default, we mask low-complexity k-mers using dustmasker v1.0.0 with default parameters prior to scoring55. Optionally, users can further reduce false positives from repetitive sequences by excluding k-mers that appear frequently across the reference transcriptome; this is achieved by counting k-mer occurrences across all transcripts and ignoring those exceeding a user-specified threshold during query, though this option is not enabled by default. Windows with scores exceeding threshold τ (default τ = 0.65) are considered positive matches. For sequences longer than w, presence in a cell is determined as soon as any single window passes the threshold.This approach tolerates sequencing errors and mutations: a single nucleotide change affects only one k-mer, allowing the remaining k-mers to drive detection. The complete query process aggregates matches across all windows, returning cells ranked by the number of significant windows. For typical use cases (w = 64, k = 24, τ = 0.65), the method reliably detects sequences with up to ~15% divergence while maintaining low false positive rates. Query optimization includes batching k-mer lookups to minimize disk access and exploiting the sorted index structure for cache efficiency. Complete specifications of the batching strategy, index access patterns, and caching parameters are provided in Supplementary Methods.Sequence coverage quantificationFor coverage analysis, we preserve spatial resolution along transcripts by reporting k-mer matches for each sliding window independently. We apply a sliding average (default window: 24 bp) to smooth the coverage profile, producing visualizations comparable to those from alignment-based tools. As with pseudoquantification, we mask low-complexity k-mers using dustmasker v1.0.055 during coverage calculation to prevent artifactual peaks. Without masking, certain positions can show 100–1,000× spikes due to sequences shared across many genes. This positional information enables detection of alternative polyadenylation sites, splice junction usage, among others.BenchmarkingWe evaluated Malva’s performance against established tools for both indexing speed and query accuracy. For indexing benchmarks, we compared Malva with kallisto-bustools v0.28.222, MetaGraph v0.4.1 (RowDiff/Multi-BRWT variant, indexing only)27, BLAST+ v2.16.0 (MEGABLAST variant, indexing only)63 and a STAR (v2.7.11b)-based workflow23,64, Fulgor v4.0.0 (meta-diff variant)31, and REINDEER v1.4.7-1-g0812ad128, on various datasets ranging sequencing depths and number of unique labels (spatial locations). For sequence search capabilities, we compared against MetaGraph, BLAST+, Fulgor and REINDEER.For benchmarking of computational complexity, dataset selection varied by method due to architectural constraints: for MetaGraph and BLAST, we used the Stereo-seq CIRSTA mouse liver atlas (CNGBdb accession: CNP0003447), containing ~60 billion reads (7 TB of raw, compressed data) across 20 samples with ~4 billion unique spatial barcodes65. REINDEER and Fulgor require treating each cellular barcode as a separate input file, making analysis of millions of cellular barcodes impractical due to filesystem overhead and preprocessing requirements (REINDEER additionally requires running bcalm66 on each sample independently). We therefore benchmarked these methods on a composite dataset from 6 Visium samples (accession numbers: GSM5115019, GSM5115020, GSM5115021, GSM5115022, GSM5115023, GSM5115024) totalling ~1.5 billion reads (~180 GB of raw, compressed data) across ~30,000 unique spatial locations67, providing comparable read depth while maintaining feasible file counts. These datasets tested the ability to handle varying spatial resolutions (0.6 µm in Open-ST to 55 µm in Visium), sequencing depths, and total data volumes. For large-scale benchmarking, we additionally selected 1% of datasets at random from Malva Index (~1.5 TB of raw data) and reindexed them with Malva and Fulgor, the two methods with best runtime and memory usage during indexing and querying, using identical settings as in other experiments. Each benchmark was run independently and was allocated at most 96 threads and 512 GB of memory on a multi-socket server with Intel Xeon Platinum 8454H CPUs (256 threads total) and 4 TB RAM, with resource monitoring via snakemake v7.32.468.Validation of biological signal preservationWe evaluated whether Malva’s k-mer-based pseudocounts produce biologically meaningful results comparable to traditional pipelines. We selected representative datasets spanning high- and lower-resolution spatial transcriptomics, healthy and disease samples: an Open-ST tumour section (high spatial resolution, small scale69; SRA accession SRX23008310), 10x v3 embryonic (E13.5) mouse midbrain (mouse single-cell, small scale70; SRA accession SRR11947658), 10x v3 human PBMCs (human single-cell, small scale; downloaded from https://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v2/pbmc_1k_v3_fastqs.tar), and Visium mouse brain (low spatial resolution, small scale, 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). Across datasets, we compared Malva pseudocounts to the original deposited count matrices. For each, we computed per-cell gene counts and calculated correlations between total counts per gene and per cell.To assess whether pseudocounts support downstream analyses, we performed clustering on all benchmarked datasets using species-specific Ensembl transcript sequences after masking low-complexity and repetitive k-mers. We applied standard scanpy processing (normalization, log transformation, highly variable gene selection, PCA, nearest neighbour graph construction, and Leiden clustering via leidenalg v0.10.2) and compared cluster assignments between Malva and traditional pipelines using Normalized Mutual Information. Additionally, the clusters from the reference objects were merged to the Malva objects by cell barcode, and used to assess the similarity marker gene detection between Malva and original counts, by applying a Wilcoxon rank-sum tests (one cluster versus the rest), then measured the concordance of the Fold Changes for each gene in each comparison.For trajectory analysis, we applied the scanpy v1.11.4 implementation of diffusion pseudotime (DPT) on the single-cell E13.5 mouse midbrain dataset71,72. Analyses were restricted to the neuronal subset, and we selected a representative radial glia cell as the root for pseudotime inference. Pseudotime ordering was computed independently on both the Malva pseudocount matrix and the reference UMI-based quantification. To enable direct comparison, we merged results by cell barcode and quantified similarity between pseudotime trajectories.Simulation of sequencing errorsTo evaluate Malva’s robustness to sequencing errors, we introduced synthetic errors into the benchmark datasets (Open-ST human tumour, 10x v3 human PBMCs, and Visium mouse brain) at rates ranging from 0.1% to 10%. These synthetic errors are cumulative to the existing errors of the data induced from the current short-read sequencing platforms. For each error rate, we randomly selected positions within reads and substituted nucleotides with uniform probability across the three alternative bases. We then rebuilt Malva indices from the error-containing reads and repeated our standard benchmarking analyses: per-cell and per-gene correlation with reference counts, clustering agreement (NMI), and marker gene detection concordance. For variant detection benchmarks, we additionally assessed sensitivity and specificity of the cell line mixture probes across error rates using both overlapping and non-overlapping k-mer indexing strategies.Sequence-based variant detection validationTo validate Malva’s accuracy for sequence-based variant detection, we compared against Flexiplex (v1.02.3), a specialized sequence search tool44, using a seven-cell-line mixture dataset (SRR10971813). The dataset contains known sequence variants: MCF7 cells expressing a BCAS4–BCAS3 fusion gene, HEK293T cells containing adenovirus E1A sequences, and T47D cells carrying mitochondrial SNV rs878887783.We searched for four target sequences: BCAS4–BCAS3 SNV variant (CCGATCCTGGGGCCGAGGTACCTTTGACAGGAGC), Adenovirus 5 E1A (TTTGGACTTGAGCTGTAAACGCCCCAGGCCATAA) and a rs878887783 probe (ATTACTCCTGCCATCATGACCCCTGGCCATAATATGATTTATCT). Compared to the Flexiplex benchmarking set-up, for the BCAS4–BCAS3 variant, we excluded the reference (CCGAGCCTGGGGCCGAGGTACCTTTGACAGGAGC), as no MCF7 cells are expected to express it. For the mitochondrial variant, we used a truncated version of the original probe as the full sequence from the Flexiplex paper yielded zero counts when following their supplementary protocol, though our modified probe produced similar detection numbers to their reported results.Flexiplex was run with default parameters using grep-based detection mode. For output processing, we extracted cellular barcodes by matching read identifiers from Flexiplex results to the corresponding R1 FASTQ file and extracting the first 16 nucleotides as cellular barcodes. For Malva, we used a sliding window size of 24 nucleotides with 100% k-mer presence threshold. Malva results directly return cell barcodes and k-mer match counts. Wall-clock time was measured for both workflows including index construction (only in the case of Malva, done once) and query processing (for each probe). Sequence search performance metrics (true positives, false positives and false negatives) were calculated by comparing detected cells against expected cell line identities for each target sequence, from the annotated objects deposited by the authors.To further validate Malva’s ability to detect clinically relevant somatic mutations, we analysed a lung adenocarcinoma dataset (PRJCA00173173) containing matched tumour–normal pairs with characterized EGFR mutation status. From the archived datasets, four tumours carried EGFR exon 19 deletions and one carried the L858R point mutation (c.2573T>G) in exon 21. We first established ground truth by aligning reads with STAR (v2.7.11b) to GRCh38 and visualizing coverage across the EGFR locus. We then designed Malva probes targeting: (1) exons 18–21 as positive controls for EGFR expression; (2) the wild-type exon 21 sequence; and (3) a 45-bp probe centred on the L858R variant (GTCAAGATCACAGATTTTGGGCGGGCCAAACTGCTGGGTGCGGAA), with a pair-matched wild-type probe (GTCAAGATCACAGATTTTGGGCTGGCCAAACTGCTGGGTGCGGAA). Probes were queried against sample-specific Malva indices with parameters w = 24, τ = 1.0 for variant-specific detection.Specialized sequence queries across Malva IndexLarge-scale arbitrary nucleotide sequence search. To systematically screen for exogenous sequences, we compiled viral genomes from NCBI RefSeq: 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). We queried these sequences against the complete Malva Index with default settings (k = 24, window size = 64, threshold = 0.65) and generated tissue- and study-level coverage profiles. We additionally screened Mycoplasma sequences because these bacteria are frequent, often undetected cell culture contaminants. Literature and guidance documents consistently report a small set accounting for most laboratory contamination events including Mycoplasma hyorhinis, Mycoplasma hominis, and Mycoplasma fermentans in human and other mammalian cell cultures, and also included Mycoplasma agassizii, a respiratory pathogen of turtles, not associated with human or standard mammalian cell culture workflows.Germline variant query probe design. To systematically screen for germline genetic variants across Malva Index, we generated sequence probes targeting common SNPs in protein-coding regions. We downloaded phase 3 data from the 1000 Genomes Project across human chromosomes 1–22 and filtered for biallelic SNPs with MAF ≥ 1% located within exons or UTRs of protein-coding genes as defined by GENCODE v19 annotations. For each qualifying variant, we extracted the genomic context from the GRCh37 reference genome and constructed 64-nucleotide probes centred on the variant position, generating paired probes containing either the reference or alternative allele.We annotated variant pathogenicity using ClinVar classifications, categorizing variants as ‘pathogenic’, ‘likely pathogenic’, ‘variant of uncertain significance’ (VUS), ‘likely benign’, ‘benign’ or ‘unknown’ based on clinical significance scores. To enable evolutionary analysis, we incorporated ancestral allele information from the 1000 Genomes Project ancestral sequence alignments, determining whether the reference or alternative allele represented the ancestral versus derived state for each variant. Quality control filters removed probes with extreme GC content (70%) to ensure reliable detection. The final probe set contained 296,985 variant pairs spanning all chromosomes.We queried these variant probes against the complete Malva Index using the standard API, retrieving cell-level match counts for both reference and alternative alleles. For evolutionary analysis, we calculated derived allele frequencies and compared observed frequencies in healthy versus disease samples to expected population frequencies, testing for systematic deviations that might indicate selection effects using the ancestral allele annotations.Isoform usage. To identify cell-type-specific isoform usage, we analysed the Tabula Muris dataset which provides full-length transcript coverage suitable for exon-level quantification. We used mouse genome annotations (GRCm39.113) and queried individual genes against sample-specific Malva indices with standard coverage analysis parameters. For each gene, we merged overlapping annotated exons to create consensus gene models, then quantified average coverage per exon across cell types.We implemented a nested marker analysis approach: first identifying marker genes using one-versus-all differential expression testing with the Wilcoxon rank-sum test and Benjamin-Hochberg multiple testing correction (log2 fold change ≥1.0, adjusted P value ≤ 0.05). For genes identified as cell-type markers, we performed exon-level differential analysis within marker cell types to identify marker exons: exonic regions showing cell-type-specific expression patterns within broadly expressed genes. This approach enabled detection of exon specificity in specific cell types despite the host gene being a marker gene shared by more cell types.Additionally, we analysed 3′ UTR length regulation using the MOSTA Stereo-seq mouse embryonic development atlas spanning E9.5–E15.5. For validation, we aligned the same sequencing reads using STAR (v2.7.11b) with default parameters to generate coverage tracks for comparison with Malva results. We performed standard Malva coverage-like analysis targeting the Add2 locus, focusing on proximal and distal polyadenylation sites previously identified in bulk and single-cell RNA-seq studies. We quantified proximal versus distal 3′ UTR usage by computing coverage ratios in defined UTR regions, then mapped these expression patterns spatially to identify tissue-specific UTR isoform preferences during embryonic development.Detection of circRNAs. To interrogate the expression of the circRNA CDR1as across Malva Index, we designed a 36-nucleotide probe centred on its backsplice junction. Therefore, the probe sequence spans 18 nucleotides upstream of the splice donor site and 18 nucleotides downstream of the splice acceptor site, capturing the unique junction sequence that distinguishes circular from linear transcripts. We used w = 24, τ = 0.65 for querying Malva Index.Benchmarking polyadenylation site detection against polyASiteTo evaluate Malva’s ability to detect alternative polyadenylation events, we leveraged polyASite v3.051, a database of polyadenylation sites identified from a large set of samples from the Human Cell Atlas (also ingested in Malva Index) using 3′ end sequencing protocols. For each annotated polyA site, we designed two probe types: (1) polyA probes consisting of 44 bp of upstream genomic sequence followed by a 20-bp polyA tract (20*A); and (2) control probes consisting of 64 bp of genomic sequence around the same position, but without the polyA signal. We queried both probe sets against the same HCA samples used to build polyASite and computed Spearman correlations between Malva pseudocounts and polyASite usage scores within each tissue type. Correlations were computed per-probe per-tissue rather than globally, as pseudocounts were not depth-normalized; this approach assesses the concordance of relative usage patterns rather than absolute quantification. Probe sequences are publicly deposited (Data availability).Large-scale somatic mutation detection and benchmark against scTMLTo systematically benchmark Malva’s mutation detection capabilities at scale, we utilized the scTML database53, which provides uniformly processed single-cell RNA-seq data from 280 tumour samples across 16 cancer types with matched variant calls. The scTML pipeline utilized GATK HaplotypeCaller for variant detection from aligned reads, with subsequent filtering for quality and read support. We extracted SNVs located in exonic regions or UTRs (excluding intronic variants due to lower sequence complexity and repeat content) and constructed 45-bp probes centred on each variant position. For each probe, we queried the corresponding sample’s Malva index and classified cells as variant-positive based on k-mer matches (w = 24, τ = 1.0). Sensitivity was computed as the fraction of scTML-positive cells detected by Malva; specificity was computed from cells in matched normal samples or variant-negative tumour cells.Reference-free analysis using k-mer compositionTraditional single-cell analysis relies on mapping reads to reference genomes and counting annotated genes. However, this approach cannot capture the full diversity of RNA sequences, particularly in non-model organisms or when studying unannotated transcripts. We developed a reference-free analysis pipeline that clusters cells based on their complete k-mer repertoire, revealing relationships invisible to gene-based methods.Bucketing approach for scalable k-mer analysisDirect analysis of k-mer composition faces a fundamental challenge: even modest datasets contain billions of unique k-mers, making cell-by-k-mer matrices computationally intractable. To address this, we implemented a bucketing strategy that groups similar k-mers while preserving sequence relationships.The bucketing algorithm uses a MinHash-based approach to assign k-mers to a fixed number of buckets based on their w-mer composition, where w 500) are used to compute the PCR bias factor (P), defined as the mean R/U ratio across cells. Malva pseudocounts (X) are then divided by P to yield corrected values (\(\hat{X}\)), mitigating inflation caused by PCR amplification. (b) Per-cell correlation between corrected Malva pseudocounts and reference UMI counts across four datasets: Open-ST human tumor, Visium mouse brain, single-cell human PBMC, and single-cell embryonic mouse brain. Red dashed lines show regression fits for uncorrected pseudocounts, highlighting improvement after correction. (c) Per-gene correlations between corrected Malva pseudocounts and reference quantification for the same datasets.Extended Data Fig. 5 Benchmarking Malva for sequence-based variant detection in mixed cell populations.(a) Query processing times for four representative targets: BCAS4–BCAS3 fusion, adenovirus 5 E1A, and mitochondrial SNV rs878887783. (b) Detection accuracy comparison between Malva and Flexiplex. Bars show true positives (solid), false positives (hatched), and false negatives (gray) for expected cell lines. Numbers indicate cell counts per category; precision values (True Positives over all Positives) are reported above each bar. (c) UMAP visualization of Malva detections across a mixture of seven cell lines. Target sequences localize to the expected backgrounds: BCAS4–BCAS3 variants in MCF7 (brown), adenovirus E1A in HEK293 (red), and rs878887783 in T47D (blue). Gray cells show no target sequence detection. (d) Same as (c), but for Flexiplex detections. (e) Precision (PPV) and (f) sensitivity (TPR, computed over the total detections by flexiplex) as a function of probe length. Solid red: computed on a Malva Index generated from overlapping k-mers; dashed red (blue) gradient: Malva Index from overlapping (non-overlapping) k-mers at 0, 0.1, 1, 2, 5, and 10% simulated error rates (light to dark). 0% error rate corresponds to the baseline, intrinsic error rate from the original FASTQ files (Illumina sequencing). None of the benchmark sequences were affected by low-complexity masking.Extended Data Fig. 6 Validation of EGFR mutation detection in lung adenocarcinoma.(a) Sample overview showing five matched tumor-normal pairs from PRJCA001731 (73), including patient sex, age, and sequencing depth. Four tumors carry EGFR exon 19 deletions; one carries the L858R point mutation. (b) STAR alignment coverage across the EGFR locus (exons 18-21) in tumor (T1-T5) and normal (N1-N5) samples. Yellow shading indicates regions with detected mismatches corresponding to known mutations. (c) Malva probe detection matrix showing pseudocounts for probes targeting exons 18-21, wild-type exon 21 position L858, and the L858R variant. Malva correctly identifies mutation-positive samples while maintaining specificity in normal tissue.Extended Data Fig. 7 Malva enables sequence queries beyond predefined gene annotations.(a) Query performance across different sequence types and lengths, ranging from short sequences (24 bp) to large genomic regions (~2 Mbp). Bubble size reflects the number of cells queried (100,000 to 100 million). (b) Dot plot showing the percentage of SARS-CoV-2–positive cells per cell type and sample (dot size) from77, with expression values z-score normalized within each cell type to highlight relative enrichment (dot color). (c) Coverage-like analysis of Mycoplasma genome across positive samples from Index shows predominant localization to 23S rRNA loci. (d) Coverage-like analysis of lentiviral vector pBA439 across positive samples from Index. (e) Per-sample heterozygosity analysis showing the fraction of cells expressing alternative alleles (mean ± s.e.), centered around 0.5 for individual samples regardless of population allele frequency. (f) Cell-type-specific differential exon usage across the Tabula Muris atlas. Representative examples highlight known regulatory switches, such as Ptprc exon inclusion in B cells and alternative splicing of Cdca2 in neuronal populations. (g) Tissue-specific detection of the circRNA CDR1as across the Human Cell Atlas and Tabula muris, compared with expression of marker genes SLC17A7 (excitatory neurons) and DLX2 (inhibitory neurons).Extended Data Fig. 8 Large-scale benchmarking of polyadenylation site detection.(a) Workflow schematic. PolyA sites from polyASite v3.051 were used to design two probe types: polyA probes (genomic sequence + polyA tract) and control probes (genomic sequence only). Both were queried against overlapping HCA samples in Malva Index (Methods). (b) Example analysis at the SPP1 locus showing polyASite annotations (top), Malva polyA probe signal (middle), and Malva control probe signal (bottom) across tissues. PolyA probes recapitulate site-specific usage patterns while control probes show broader coverage-like profiles. (c) Spearman correlation between Malva pseudocounts and polyASite usage scores across ~1,000 randomly selected genes with annotated polyA sites (totaling ~1.5 M probe-tissue pairs across ~600 samples). Point size depicts the relative number of polyA sites per gene; point color shows the difference in correlation (polyA minus control), showing consistent improvement with polyA-specific probe design.Extended Data Fig. 9 Pan-cancer dataset composition and mutation detection completeness.(a) scTML dataset composition showing sample counts per cancer type. (b) Detection completeness heatmap showing sensitivity across cancer types and mutation frequency bins. OS, osteosarcoma; GBM, glioblastoma; HNSC, head and neck squamous cell carcinoma; LUAD, lung adenocarcinoma; KIRC, kidney renal clear cell carcinoma; PAAD, pancreatic adenocarcinoma; OV, ovarian cancer; STAD, stomach adenocarcinoma; BLCA, bladder carcinoma; CHOL, cholangiocarcinoma; PRAD, prostate adenocarcinoma; THCA, thyroid carcinoma; LIHC, liver hepatocellular carcinoma; OSCC, oral squamous cell carcinoma; CRC, colorectal cancer; UVM, uveal melanoma.Extended Data Fig. 10 De novo sequence assembly with cluster-specific k-mer subsets.(a) Schematic of the workflow for identifying marker k-mers from Malva-derived clusters. Cluster-enriched k-mers are extracted and used either for transcript quantification (gene-based mode) or for de novo assembly (k-mer–based mode). (b) Comparison of sequence-based versus gene-based clustering across healthy human cells, as shown in Fig. 5b. Latent structure agreement (AUC) and Adjusted Mutual Information (AMI) quantify how well sequence-derived embeddings preserve cluster structure relative to gene counts. Lower agreement in single-nucleus data reflects low-complexity intronic sequences. N = 29,240,986 cells examined over 4,932 sequencing runs; boxplots show the median, interquartile range, and whiskers extending to 1.5x the interquartile range. (c) Parameter scan for k-mer bucketing. Varying k-mer size and bucket definitions shows the bounds of informative sequence length required to recover meaningful latent structures. Color scale is clipped to 0.4 as the cutoff for low-medium quality clustering agreement. (d) Efficiency of cluster-specific assembly. Gene detection rate is shown for transcripts reconstructed with SPAdes using either all reads or cluster-specific Malva k-mer subsets. (e) Annotation versus assembly differences. Histogram shows genomic distance between assembled contigs and annotated transcript 3′ ends (TES). Inset: proportion of contigs mapping within different distances from the TES. (f) Mapping of contigs with minimap2 to the human genome, and subsequent classification of unmapped contigs with kraken2. Taxonomic assignment reveals human intergenic sequences, bacterial species and a residual pool of sequences with no confident match.Supplementary informationSupplementary Information (download PDF )This file contains Supplementary Methods, Supplementary Discussion, Supplementary Notes 1–4, Supplementary Figs. 1–4, Supplementary Tables 1–4 and references. The Supplementary Methods detail the data structures and algorithms used for indexing and querying with Malva. The supplementary figures contain benchmarking and evaluation experiments of Malva indexes and the query engine, and the metadata ontologies. The supplementary tables indicate the different types of query, and benchmark experiments of the Malva algorithmReporting Summary (download PDF )Peer Review file (download PDF )Supplementary Code (download ZIP )source code of the Malva indexing and querying algorithmsRights and permissionsOpen Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.Reprints and permissionsAbout this article