Joint analysis of expression and variation at single cell resolution by scVar

Joint analysis of expression and variation at single cell resolution by scVar | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Joint analysis of expression and variation at single cell resolution by scVar View ORCID Profile Wei Zhao , View ORCID Profile Fei Yang , View ORCID Profile Wenting Zong , Xinchang Zheng , View ORCID Profile Yiming Bao doi: https://doi.org/10.1101/2025.10.12.681753 Wei Zhao 1 National Genomics Data Center, China National Center for Bioinformation , Beijing 100101, China 2 Beijing Institute of Genomics, Chinese Academy of Sciences , Beijing 100101, China 3 University of Chinese Academy of Sciences , Beijing 100049, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Wei Zhao Fei Yang 1 National Genomics Data Center, China National Center for Bioinformation , Beijing 100101, China 2 Beijing Institute of Genomics, Chinese Academy of Sciences , Beijing 100101, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Fei Yang Wenting Zong 1 National Genomics Data Center, China National Center for Bioinformation , Beijing 100101, China 2 Beijing Institute of Genomics, Chinese Academy of Sciences , Beijing 100101, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Wenting Zong Xinchang Zheng 1 National Genomics Data Center, China National Center for Bioinformation , Beijing 100101, China 2 Beijing Institute of Genomics, Chinese Academy of Sciences , Beijing 100101, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: baoym{at}big.ac.cn Yiming Bao 1 National Genomics Data Center, China National Center for Bioinformation , Beijing 100101, China 2 Beijing Institute of Genomics, Chinese Academy of Sciences , Beijing 100101, China 3 University of Chinese Academy of Sciences , Beijing 100049, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yiming Bao For correspondence: baoym{at}big.ac.cn Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Tumor heterogeneity represents a critical determinant affecting both cancer diagnosis and therapeutic efficacy. Traditional bulk sequencing approaches are limited in their ability to resolve genomic alterations at the resolution of individual cells. In this study, we present scVar, an integrated analytical framework designed for mutation profiling using single-cell transcriptomic data. scVar enables sensitive and robust detection of single-nucleotide variants, particularly facilitating the identification of low-frequency variants. The framework incorporates customizable filtering parameters including variant allele frequency, and provides comprehensive functional and clinical annotations for identified mutations. Furthermore, scVar includes downstream analytical modules that facilitate the joint investigation of transcriptomic and mutational profiles within single cells. Benchmarking on simulated datasets and matched tumor samples demonstrates that scVar consistently outperforms existing methods, especially in detecting variants of low abundance, and shows strong concordance with whole-exome sequencing data. Application of scVar to non-small cell lung cancer samples effectively characterized ITH across spatially distinct tumor regions and diverse cellular populations. Collectively, scVar offers an integrated platform for concurrent analysis of somatic mutations and transcriptomic data derived from single-cell RNA sequencing, providing a powerful tool for elucidating tumor evolution and the complex interplay between genomic heterogeneity and the tumor microenvironment. Background Cancer continues to represent a significant global health challenge, accounting for approximately 20 million new cases and 10 million deaths in 2022, as reported by the GLOBOCAN 2022 database ( https://gco.iarc.fr/en ). A major obstacle in the advancement of cancer therapeutics is the complexity of tumor heterogeneity, a multifactorial phenomenon that compromises diagnostic accuracy, therapeutic efficacy, and prognostic predictions[ 1 ]. Tumor heterogeneity is characterized by inter-tumor heterogeneity, which denotes the differences among tumors of the same type across different patients, and intra-tumor heterogeneity, which refers to the coexistence of multiple cancer cell subpopulations with distinct genotypes, phenotypes, and functions within a single tumor of the same patient. Intertumoral heterogeneity (ITH), compromises precise cancer diagnosis and effective therapy[ 2 , 3 ]. ITH has been identified as a prognostic factor for reduced survival across multiple cancer types[ 4 ] and a key driver of tumor evolution and therapeutic resistance[ 5 – 7 ]. Advances in high-throughput sequencing platforms, including whole-exome sequencing (WES) and whole-genome sequencing (WGS), have enabled quantitative characterization of tumor heterogeneity. These approaches facilitate the reconstruction of tumor clonal architectures and provide insights into evolutionary trajectories and mechanisms underlying drug resistance[ 8 – 10 ]. However, traditional bulk sequencing approaches lack single-cell resolution, limiting the accurate depiction of subclonal structures and constraining our understanding of the spatial and temporal dynamics of tumor heterogeneity. The development of single-cell sequencing has made it possible to extend the study of genes and their functional alterations to the single-cell level[ 11 ]. Single-cell genome sequencing is the best way to identify somatic mutations in single cells. However, since each cell contains only two copies of DNA, single-cell DNA sequencing (scDNA-seq) faces many challenges, such as amplification bias and low coverage 9 . These can reduce the quality of scDNA-seq and increase the difficulty of mutation identification. Alternatively, somatic mutations can be detected directly from high-throughput single-cell assays, such as single-cell RNA sequencing (scRNA-seq). Notably, 10X Genomics scRNA-seq is widely adopted in transcriptome analyses due to its high throughput, cost efficiency, high mRNA detection sensitivity, and robust capacity to capture rare cell populations[ 12 – 14 ]. Furthermore, this technology has advanced our understanding of tumor biology by elucidating mechanisms underlying tumor heterogeneity, including clonal evolution dynamics[ 15 , 16 ], interactions within the immune microenvironment[ 17 , 18 ], and transcriptional plasticity among cellular subpopulations[ 19 ]. Importantly, scRNA-seq datasets contain latent genomic information, including somatic mutations, which remain underexploited due to methodological constraints[ 20 ]. Existing computational tools for mutation detection from 10X Genomics scRNA-seq data, including cellSNP[ 21 ], Monopogen[ 22 ], and SComatic[ 23 ], exhibit certain limitations. These include restricted accuracy stemming from the resolution of cell type annotations, insufficient sensitivity for detecting low-frequency mutations, and the lack of comprehensive software platforms that integrate single-cell transcriptomic data with mutation information. To overcome these challenges, we developed scVar, a framework that extends mutation detection capabilities from bulk transcriptomic data to the single-cell level. scVar facilitates de novo identification of single-nucleotide variants (SNVs) from 10X Genomics scRNA-seq data without the necessity for paired normal samples or predefined cell-type annotations. The framework integrates seamlessly into standard 10X Genomics workflows and supports comprehensive single-cell transcriptome analysis, SNV detection, as well as downstream comprehensive analysis and visualization. This enables users to correlate transcriptional heterogeneity with genetic variation at single-cell resolution and to investigate the mutation landscape in details. Results Overview of scVar scVar facilitates the comprehensive identification and analysis of somatic single nucleotide mutations from scRNA-seq data through a multi-layered analytical framework ( Fig. 1 ). The workflow initiates with preprocessing and fundamental analysis of 10X scRNA-seq data, encompassing alignment, noise reduction, clustering, and cell-type annotation. The core stage of scVar involves mutation detection and annotation. During this phase, high-quality binary alignment map (BAM) files are reconstructed using cell barcodes, and candidate mutations are identified through GATK[ 24 – 26 ], optimized with parameters specifically adapted for mutation detection in 10X scRNA-seq datasets (Supplementary Note and Supplementary Fig. 1). Subsequently, genotype assignment is conducted at the single-cell level through parallelized computation, accompanied by noise reduction employing stringent filtering criteria. To reduce false positives, scVar applies rigorous filters, including the requirement that mutations be observed in at least two cells and the exclusion of common single nucleotide polymorphisms (SNPs) based on population frequency databases. Functional annotation of identified mutations is performed using the Ensembl Variant Effect Predictor (VEP)[ 27 ] and Annovar[ 28 ], supplemented by clinical relevance assessments referencing Cancer Hotspots[ 29 , 30 ], COSMIC[ 31 ], ClinVar[ 32 ], and other databases. Following this, scVar integrates mutation data with single-cell transcriptomic gene expression profiles into a MuData format file, thereby enabling streamlined and integrated downstream analyses. The analytical pipeline of scVar comprises three principal modules. The first module, mutation profiling and heterogeneity assessment, involves characterization of global and cell population-specific mutation signatures, evaluation of tumor mutational burden (TMB), assessment of clonal diversity using entropy and Simpson’s diversity index, and clonal clustering based on binary sparse mutation matrix. The second module focuses on functional annotation, encompassing the identification of cell population-specific mutations, enrichment analyses of oncogenic pathways and Gene Ontology (GO) terms, as well as pseudotime analysis to characterize the developmental states of mutated cells within the cellular differentiation trajectories. The third module provides statistical summaries such as estimation of cell type proportions, evaluation of alignment quality metrics, and visualization of the most recurrent mutations by mapping mutated cells onto Uniform Manifold Approximation and Projection (UMAP) embeddings with enhanced annotation. Download figure Open in new tab Figure 1. Overview of scVar. Schematic overview of the analytical workflow, illustrating the initial stages of data preprocessing, sequence alignment, quality control, fundamental single-cell transcriptomic analyses, and cell type annotation, followed by the principal steps of mutation processing, including somatic variant identification, genotyping, variant filtering, and functional annotation, and culminating in the integrated analytical modules implemented within the scVar framework. In summary, scVar presents a comprehensive and integrated framework for the detection and analysis of SNVs derived from scRNA-seq data. It facilitates robust downstream analyses aimed at elucidating genomic heterogeneity, mutation dynamics, and cell-type-specific functional drivers. Furthermore, scVar can be seamlessly integrated into the standard 10X Genomics pipeline, thereby furnishing additional insights at the somatic mutation level. Performance evaluation of somatic mutation detection using simulated data To establish a robust evaluation framework for assessing somatic mutation detection performance in 10X Genomics scRNA-seq data, we designed and implemented a comprehensive in silico spike-in experiment using five matched lung tissue samples. These baseline datasets were preprocessed by globally reverting all detectable variants to reference alleles. High-confidence candidate loci were selected with at least 10X coverage. This was followed by in silico simulated mutations at controlled variant allele frequency (VAF) ranging from 10% to 50%, achieved by read-level modifications. Quantitative benchmarking against alternative single-cell variant callers—including cellSNP, Monopogen, and SComatic—revealed that scVar consistently outperformed these methods across all evaluated VAF thresholds. Specifically, scVar achieved the highest recall rates and demonstrated robust F1-scores, surpassing the performance of the second-best method, Monopogen. This performance advantage was especially notable at lower VAF thresholds, where scVar exhibited enhanced sensitivity in detecting rare somatic variants, thereby confirming its efficacy in resolving rare mutations from sparse transcriptomic data ( Fig. 2a ). Download figure Open in new tab Figure 2. Comparison of the performance of scVar against other mutation detection methods. a. Performance of cellSNP, Monopogen, SComatic, scVar for the detection of somatic mutations in the simulation data. b. Performance of cellSNP, Monovar, SComatic, scVar for the detection of somatic mutations in the scRNA-seq data from cSCC. Significance with respect to scVar using the Student’s t-test ( ***P < 0.0001). c. Comparison between the mutational spectra of the mutations identified in cSCC samples using WES and scRNA-seq data across the benchmarked algorithms. The cosine similarities between the mutational spectra derived from scRNA-seq and WES data are illustrated. d. Mutational rates for all samples using the somatic SNVs detected by scVar in cSCC and matched normal skin scRNA-seq data sets. e. Correlation between the mutation rate estimated using the mutations identified in WES and the mutations found in epithelial cells detected by scVar in the scRNA-seq data. The correlation was assessed using a linear regression model. Only genomic regions with sufficient sequencing depth available in both the WES and scRNA-seq datasets were included in this analysis. Mb, megabase. Validation of scVar using scRNA-seq and WES data We further compared scVar with other methods using real data. Specifically, we employed matched WES and scRNA-seq data from eight cutaneous squamous cell carcinoma (cSCC) tumors[ 33 ]. Mutations characterized by high coverage in the WES data and concurrently detected in the corresponding scRNA-seq datasets were designated as the gold standard benchmark. This benchmark encompassed a total of 14,606 mutations located within regions of high sequencing coverage. These mutations, located in high-coverage areas, showed a positive correlation in VAF between WES and scRNA-seq (R = 0.52, P = 3.94×10^(-132), Supplementary Figs. 2a, 2b). Recall was operationally defined as the proportion of true mutations accurately identified by each computational tool, whereas precision was calculated based on mutations uniquely reported by the software but absent in both WES and scRNA-seq BAM files (see Methods). Benchmarking outcomes revealed that scVar significantly outperformed existing approaches. Compared to Monopogen, the second-best performing method, scVar achieved markedly higher recall rates (0.61–0.90 compared to 0.54–0.68; P = 0.00079, two-sided t-test) and F1 scores (0.75–0.92 versus 0.69–0.80; P = 0.0017), while consistently maintaining high precision values (0.89 to 0.98) across all evaluated conditions ( Fig. 2b ). Subsequent performance evaluations were stratified by mutation subtypes and VAFs. Initially, mutations within the gold-standard dataset were categorized into six single-nucleotide substitution classes. scVar exhibited superior detection performance independent of subtype, achieving significantly higher recall rates (0.78–0.82) compared to the next-best performers, SComatic (0.44–0.48; P < 0.0004, two-sided t-test) and Monopogen (0.39–0.46; P < 0.0004), with no significant detection bias observed across single-nucleotide substitution subtypes (Supplementary Fig. 3a). Subsequently, we analyzed the VAF values of gold-standard mutations in WES tumor samples, comparing the recall rates of various software across different VAF intervals. In the critical ≤10% VAF range mutations, scVar achieved comparable performance to cellSNP. Furthermore, in the 10%–20% VAF range, scVar maintained superior recall value compared to SComatic, highlighting scVar’s efficacy in detecting low-frequency mutations (Supplementary Fig. 3b). Cosine similarity analysis of the mutation spectra further validated the correctness of mutation detection by scVar, demonstrating a high concordance with those identified by WES, with a cosine similarity score of 0.79 ( Fig. 2c ). To validate scVar’s capability not only to detect a greater number of somatic mutations in cancer but also to resolve these mutations at the level of individual cell types, we analyzed mutation profiles from paired normal and tumor samples. This comparison revealed elevated mutation rates in tumor samples relative to matched normal tissues, with the exception of samples P2 and P10 ( Fig. 2d ). Additionally, mutation rates estimated from scRNA-seq data in epithelial cells exhibited a strong positive correlation with mutation burdens inferred from the WES data (R = 0.75, P = 0.032, linear regression; Fig. 2e ), indicating that scVar facilitates precise quantification of mutation burdens at cell-type resolution. Overall, these results demonstrate that scVar is a robust and effective tool for detecting somatic mutations in scRNA-seq data, outperforming established methods including cellSNP, Monopogen, and SComatic in terms of recall and F1 score. Its proficiency in accurately identifying mutations, particularly those with low VAF, alongside its consistent mutation spectrum concordant with WES data, underscores its reliability. Furthermore, scVar’s ability to quantify mutation burdens at the resolution of individual cell clusters offers novel insights into tumor heterogeneity, thereby advancing the understanding of cancer biology. Explore the heterogeneity at the sampling sites level and single-cell level in the non-small cell lung cancer (NSCLC) samples Mutations Landscape of Different Sampling Sites We analyzed five treatment-naïve, non-metastatic NSCLC patients undergoing curative-intent lobectomy[ 34 ]. Tumor specimens were dissected into core, mddle, and edge regions, alongside distal normal tissue located more than 5 cm from the tumor margin. These samples were dissociated into single-cell suspensions and subjected to droplet-based scRNA-seq. Due to insufficient yield (<500 cells) from the tumor core regions of patients P1 and P2, subsequent analyses were restricted to patients P3–P5: patient P3 (68-year-old male, stage IIIB lung adenocarcinoma), patient P4 (64-year-old female, stage IIB lung adenocarcinoma), and patient P5 (60-year-old male, stage IA3 large cell carcinoma). The scVar computational tool demonstrated superior performance in resolving spatial heterogeneity across these three NSCLC cases. Clustering analysis of the mutation-cell binary matrix generated by scVar revealed distinct mutational landscapes among alveolar cells and T cells from patient P5 across different sampling sites ( Figs. 3a , 3b ). While transcriptional profiles exhibited limited spatial compartmentalization, mutational patterns displayed marked spatial segregation, indicating that genomic evolution occurs independently of transcriptional states. To thoroughly characterize spatial genomic heterogeneity, TMB was measured across distinct anatomical compartments within alveolar cells (tumor cells) and T cells. In both cell populations, the intermediate and peripheral regions demonstrated elevated TMB compared to the tumor core. Specifically, alveolar cells exhibited TMB values of 1.67 mutations per megabase (mutations/Mb) in the core, increasing to 2.33 and 2.37 mutations/Mb in the intermediate and edge regions, respectively. Similarly, T cells showed TMB values of 2.84 mutations/Mb in the core, rising to 3.99 and 3.80 mutations/Mb in the intermediate and edge compartments, respectively. These findings indicate the presence of distinct microenvironmental selection pressures across different tumor subregions ( Fig. 3c ). Download figure Open in new tab Figure 3. Analysis of spatial intra-tumor heterogeneity across tumor regions based on somatic mutations detected by scVar in scRNA-seq data. a. Hierarchical clustering of alveolar cells from all regions (represented in columns) based on somatic mutations (indicated in rows). Mutations identified in the cells are highlighted in red, while white signifies the absence of mutations. b. Hierarchical clustering of T cells across all regions (shown in columns) with somatic mutations (represented in rows). c. A line plot illustrating TMB in both alveolar cells and T cells across different sampling regions. d. Box plots displaying entropy and Simpson’s index for T cells at various tumor sites. e. GO enrichment results for mutated genes in alveolar cells (on the right) and T cells (on the left) across different sampling regions. Quantitative assessment using entropy and Simpson’s diversity index revealed a progressive increase in T cell clonal diversity from the tumor core (mean entropy: 0.24; mean Simpson’s diversity index: 0.93) to the tumor edge (mean entropy: 0.33; mean Simpson’s diversity index: 0.89), with intermediate zones exhibiting transitional values ( Fig. 3d ). Comparable trends were observed in the other two patients (Supplementary Fig. 4). This gradient pattern aligns with microenvironmental pressure gradients, wherein hypoxic tumor cores favor dominant clones, while immune-enriched tumor margins support greater subclonal diversity, potentially mediated by differential exposure to microenvironmental factors and external selective pressures. Functional partitioning was further elucidated through GO enrichment analysis of compartment-specific mutations. In alveolar cells, mutation-associated pathways exhibited distinct spatial heterogeneity across tumor subregions, with location-specific molecular signatures correlating with key oncogenic processes. Across the three tumor regions, pathways involved in the positive regulation of I κ B kinase (IKK)/ nuclear factor kappa B (NF- κ B) signaling were significantly enriched (adjusted p -values ranging from 0.005 to 0.037, hypergeometric test), consistent with prior findings that underscore NF- κ B’s established role in promoting tumor-associated inflammation and metastasis through sustained activation in lung carcinogenesis[ 35 , 36 ]. Both tumor core and intermediate regions exhibited enrichment of negative regulation of apoptotic signaling, suggesting systemic suppression of cell death mechanisms to sustain malignant survival. The unique enrichment of pathways associated with the positive regulation of mRNA processing in the tumor core may be attributable to its characteristic hypoxic and acidic microenvironment, as hypoxia has been shown to promote pre-mRNA splicing in the tumor microenvironment, leading to oncogene activation, tumor suppressor inactivation, metabolic adaptation, and angiogenic reprogramming[ 37 ]. In contrast, the tumor intermediate regions exhibited enrichment in multiple pathways linked to tumor microenvironment remodeling and metastasis, including: (1) epidermal growth factor response (EGF R ), indicative of EGFR-mediated proliferative/migratory pathway hyperactivation characteristic of adenocarcinomas; (2) chemokine production, associated with epithelial-mesenchymal transition (EMT) induction and cancer stem cell niche formation[ 38 ]; and (3) co-enrichment of c-Jun-N-terminal kinase (JNK) cascade and p38 mitogen-activated protein kinases (p38MAPK) cascade, mechanistically linked to hypoxia-driven invasion via NOTCH4-mediated pathway crosstalk[ 39 ]. These activated pathways collectively enhance the invasive potential of lung adenocarcinoma cells. T cells isolated from the middle region of the tumor exhibited a specific enrichment in pathways associated with cancer development, particularly those regulating apoptotic signaling and intrinsic apoptotic pathways. Additionally, this region also demonstrated enhanced autophagy of mitochondria, which is essential for driving effector T cell differentiation through metabolic reprogramming[ 40 ]. This mechanism plays a significant role in the immunological regulation of tumor progression and the effectiveness of therapeutic responses[ 41 ] ( Fig. 3e ). In summary, the application of scVar enabled detailed investigation of the genomic features in three NSCLC patients, thereby validating previously published observations at the single-cell level—an analytical resolution unattainable with conventional bulk genomic sequencing approaches. Mutations Landscape of Different Cell Types The mutational landscape across distinct compartments of the tumor (core, middle, and edge) exhibited cell type-specific functional enrichment patterns that closely align with the intrinsic biological roles of these cells and their contributions to tumorigenesis. In patient P3, mutated genes within five predominant cell types—alveolar cells, B cells, fibroblasts, monocytes, and T cells—each comprising over 500 cells, were enriched in pathways related to the regulation of DNA damage response, mRNA splicing, and mRNA processing. These perturbations collectively contribute to genomic instability and aberrant transcriptional programs, which are hallmark features of cancer progression. Notably, alveolar cells demonstrated distinct enrichment of the EGFR signaling pathway and ErbB signaling pathway. The EGFR pathway, a well-established oncogenic driver in lung cancer, operates synergistically with other ErbB family members to activate critical intracellular signaling cascades that promote cell survival and proliferation[ 42 ]. Furthermore, alveolar cells showed functional enrichment in pathways governing mitochondrial gene expression, a process intimately linked to tumor progression. Dysregulated expression of mitochondrial respiratory genes— commonly suppressed in lung cancer tissues relative to normal tissues—has been directly correlated with tumor aggressiveness and poor prognosis. Notably, alveolar cells were also enriched for processes involved in surfactant homeostasis. Homeostasis of lung surfactants is closely related to alveolar cell function. Type II alveolar epithelial cells synthesize and secrete lipids and proteins of lung surfactants that maintain lung function[ 43 ]. While immune response-regulating signaling pathways were universally enriched among immune cells. Monocytes uniquely activated toll-like receptor 2 (TLR2) signaling and positively regulated NF- κ B-inducing kinase (NIK)/NF- κ B signaling, which are associated with pro-tumorigenic inflammation and immune evasion[ 44 , 45 ]. TLR2 activation via NF-κB has been reported to paradoxically promote proliferation in lung adenocarcinoma cells. Furthermore, monocytes displayed enrichment in pathways regulating phagocytosis and monocyte differentiation, which are closely related to monocyte function[ 46 ]. T cells are also enriched with antigen receptor-mediated signaling pathways and alpha-beta T cell activation associated with that cell type ( Fig. 4a ). Download figure Open in new tab Figure 4. Analysis of intra-tumor heterogeneity across different cell types using somatic mutations detected by scVar in scRNA-seq data. a. GO enrichment analysis results of different cell types of mutant genes in three tumor samples from patient P3. b. The top inset presents the t-SNE plot of the 52,698 profiled cells annotated by associated cell type, while the bottom inset shows cells colored according to KRAS G12A mutation status as identified by scVar. c. UMAP visualization of alveolar cells. The left panel depicts the subclustering of alveolar cells, while the right panel shows the same cells colored according to KRAS G12A mutation status as identified by scVar. d. Pseudotime analysis of alveolar cells and proportion of cell sources at critical nodes shown in left. The right inset shows the position of KRAS G12A mutation cells in the trajectory and proportion of cell clusters at different nodes. e. Split violin plots illustrating the expression levels of genes associated with the RAF-MEK-ERK and PI3K-AKT-mTOR pathways in two alveolar cell subclusters. f. Stacked violin plots showing differentially expressed genes in KRAS G12A mutation cells compared to other alveolar cells. Through analysis of mutations across distinct cell types, we observed that certain cancer driver mutations may exhibit cell type–specificity, thereby shaping their oncogenic potential within distinct cellular contexts. The KRAS gene, a pivotal member of the RAS oncogene family, represents the most frequently mutated oncogene across various malignancies[ 47 ]. Oncogenic KRAS mutations, including G12A, G12C, and G12D, result in constitutive activation of KRAS , driving aberrant signaling through pathways like MAPK and PI3K-AKT[ 48 ]. These mutations are associated with tumorigenesis, therapeutic resistance, and poor clinical outcomes in multiple cancers, including lung, colorectal, and pancreatic carcinomas[ 49 , 47 , 50 ]. The KRAS G12A mutation, although less common, is clinically significant and has been linked to reduced progression-free and overall survival in metastatic colorectal cancer, underscoring its role as a potent oncogenic variant[ 51 ]. According to the Cancer Hotspots database, 166 cancer patients, including six with lung cancer, harbor the KRAS G12A mutation[ 29 , 30 ]. In tumor samples from three NSCLC patients, 64 alveolar cells bearing KRAS G12A mutations were identified exclusively in patient P3 ( Fig. 4b ). Importantly, all KRAS -mutated cells were confined to the alveolar cell population, highlighting the cell type-specific nature of this oncogenic driver. Spatial analysis revealed an uneven distribution of KRAS G12A-mutated cells: 43 (67.2%) localized to the tumor core, 18 (28.1%) to the intermediate region, and 3 (4.7%) to the tumor periphery, suggesting enrichment within hypoxic or high-stress microenvironments. Transcriptional profiling demonstrated that KRAS G12A-mutated alveolar cells possess distinct gene expression signatures relative to their wild-type counterparts, indicative of mutation-driven transcriptional reprogramming. All KRAS G12A-mutated cells clustered uniquely within Cluster 0, exhibiting clear separation from wild-type alveolar cells in UMAP visualization ( Fig. 4c ). Pseudotime trajectory analysis positioned KRAS G12A-mutated alveolar cells predominantly at the terminal branch of the developmental continuum ( Fig. 4d ), suggesting a later differentiation stage compared to wild-type cells. This observation is similar to previous findings in KRAS G12D-mutated alveolar cells[ 52 ]. Canonical KRAS signaling activates the RAF-MEK-ERK cascade and also modulates the PI3K-AKT-mTOR pathway, a critical axis regulating cell proliferation, differentiation, apoptosis, and glucose metabolism, with profound implications for tumorigenesis and therapeutic resistance. Comparative analysis revealed that clusters harboring the KRA S G12A mutation exhibited elevated expression of genes associated with both the RAF-MEK-ERK and PI3K-AKT-mTOR pathways relative to control clusters ( Fig. 4e ), suggesting hyperactivation of these oncogenic signaling networks[ 47 ]. Differential gene expression analysis identified multiple genes implicated in lung cancer and other malignancies as significantly upregulated in KRAS G12A-mutant cells compared to other alveolar cells ( Fig. 4f ). Notably, ATP1A1 , NMT1 , IFITM1 , and BPTF were markedly overexpressed. These genes have established associations with lung cancer pathogenesis. Specifically, BPTF promotes tumor growth and correlates with poor prognosis in lung adenocarcinoma[ 53 ], while IFITM1 drives tumor cell metastasis—effects demonstrably inhibited by its silencing in lung cancer models[ 54 ]. Additionally, KRTCAP2 displays the most pronounced differential expression, with established roles in hepatocellular carcinoma (HCC)[ 55 ]. Other cancer-relevant differentially expressed genes (DEGs) include the recognized oncogene RPN2 , which enhances the malignant properties of tumor cells and affects multidrug resistance[ 56 ]; REEP3 , overexpressed in pancreatic cancer and correlated with reduced overall survival and disease-free survival[ 57 ]; The ARF family is generally highly expressed across pan-cancer, and ARF5 is significantly associated with poor prognosis in HCC patients[ 58 ]. This comprehensive profile underscores the multifaceted oncogenic signatures associated with KRAS G12A mutations. In summary, relative to conventional single-cell transcriptome analyses, scVar enables the integration of somatic mutation data with transcriptomic profiles at single-cell resolution. This approach facilitates a more comprehensive investigation of the heterogeneity among diverse cell types and subpopulations within tumors. Consequently, it yields novel insights into tumor heterogeneity and presents promising opportunities for the development of targeted therapeutic interventions. Discussion This study presents scVar, a comprehensive framework capable of detecting SNVs from 10X scRNA-seq data and performing joint analyses of gene expression and genetic variation at single-cell resolution. Building on softwares established for bulk transcriptomics, scVar extends mutation detection to 10x Genomics scRNA-seq data. Importantly, scVar demonstrates enhanced sensitivity in identifying somatic mutations, particularly low-frequency variants, without the necessity for paired normal samples or predefined cell-type annotations. The software further incorporates comprehensive downstream analytical modules that facilitate the exploration of heterogeneity among distinct cellular subpopulations. These modules include the characterization of mutational signatures, quantification of tumor mutational burden and clonal diversity, and functional enrichment analyses of mutated genes. Through the joint analysis of mutational and transcriptomic features, scVar can reveal differential gene expression programs in mutated cells, delineate their trajectories within clonal and developmental lineages, and correlate transcriptional heterogeneity with genetic variation at single-cell resolution, thereby providing a new perspective for the study of tumor heterogeneity. SNVs, particularly low-frequency variants, often play an essential role in tumor evolution and heterogeneity. In benchmarking with both simulated datasets and eight cSCC tumors with matched WES and scRNA-seq, scVar consistently outperformed comparator methods, achieving higher recall and F1 scores with enhanced sensitivity to rare mutations. Integrating genomic and transcriptomic information at single-cell resolution facilitates direct mapping of genotype-to-phenotype relationships, thereby advancing the understanding of tumor heterogeneity. scVar can be seamlessly incorporated into standard 10X Genomics workflows, enabling the concurrent analysis of gene expression and mutation data, which enhances mutation-level insights and broadens its practical applicability. Application of scVar to real datasets has enabled detailed investigations of mutational characteristics across diverse cell populations, facilitated the examination of transcriptional programs in mutated cells such as those harboring oncogenic KRAS mutations, and supported the tracking of mutated cells along cellular trajectories, thereby providing deeper insights into intratumoral heterogeneity. Despite these advances, several limitations warrant consideration and motivate future developments. First, somatic mutations present in a limited subset of single cells may evade detection by WES or WGS, resulting in reduced cosine similarity between the mutation spectra derived from scVar and WES data. Second, current scRNA-seq technologies face challenges including data sparsity, elevated dropout rates, and inadequate gene expression levels. The lack of full-length transcript coverage in platforms such as 10X Genomics platform further complicates data completeness and mutation detection accuracy. To overcome these issues, future developments of scVar will focus on enhancing compatibility with full-length single-cell transcriptome protocols and long-read sequencing technologies, thereby improving variant identification and enabling the detection of alternative splicing events. Additionally, systematic sequencing artifacts and platform-specific errors continue to contribute substantially to false-positive variant calls. Consequently, the development of a Panel of Normals (PON) tailored for 10X scRNA-seq data, which catalogs recurrent artifact sites and models background noise, represents a promising approach to enhance variant calling accuracy. Moreover, copy number variations (CNVs) and structural variations (SVs) are critical contributors to the genomic instability characteristic of cancer cells, influencing tumorigenesis, tumor behavior, and therapeutic response. To address these challenges, future iterations of scVar aim to incorporate CNV and SV detection capabilities, thereby enabling a more comprehensive mutational analysis and enhancing the overall applicability and analytical robustness. Conclusions scVar is an integrated workflow designed to detect SNVs from 10x Genomics scRNA-seq data and to perform joint analyses of gene expression and genetic variation at single-cell resolution. Evaluation using simulated datasets as well as matched WES and scRNA-seq data demonstrates that scVar outperforms existing methods in identifying somatic mutations from scRNA-seq, exhibiting notably enhanced sensitivity for variants with low VAF. Application of scVar to cancer datasets enables the seamless integration of somatic mutation detection with transcriptomic profiles at the single-cell level, thereby providing more nuanced insights into tumor heterogeneity. Methods Detection of somatic mutations in single-cell data sets using scVar Initial read mapping and scRNA-seq expression analysis During the initial stage of the scVar protocol, sequenced reads are aligned to the reference genome using the count tool of Cell Ranger[ 59 ] Count (v.7.1.0) under default parameters. To reduce the impact of cell-free RNAs that constitute background noise, scVar employs soupX[ 60 ] (R v.4.1.2 and soupX v.1.6.2) to estimate and remove ambient RNA contaminations. A seurat object is created using CreateSeuratObject in Seurat[ 61 ] (v.4.3.0.1) based on the cleansed output from soupX with the parameters min.cells and min.features set to 3 and 10 respectively. Cells that contain greater than 10 and less than 15,000 expressed genes, fewer than 50% ribosomal transcripts, and fewer than 10% mitochondrial transcripts are saved. The data matrix undergoes a standardization process through NormalizeData, utilizing a scale.factor of 1e4. Top 2,000 variable features are identifird by FindVariableFeatures using variance stabilization transformation (“vst”). After further scaling data by ScaleData, Principal component analysis (PCA) is performed with the variable features by RunPCA. Unsupervised graph-based clustering is presented by FindNeighbors with 30 reduction dimensions and FindClusters with the resolution set to 0.7. SingleR[ 62 ] (v.1.8.1) is used to detect cell types based on Human Primary Cell Atlas data as the reference. Mutation detection and annotation To minimize the impact of low-quality cells on subsequent analyses, the reads are extracted from BAM files to create new BAM files based on the barcodes of high-quality cells identified after quality control. First, scVar separates the BAM files by chromosome using the split tool from Bamtools[ 63 ] (v.2.5.1), and employs the pysam package (Python v.3.11.4) to extract the reads corresponding to high-quality cells identified by the ’CB’ tag. The Multiprocessing package is utilized for parallel processing. The extracted BAM files are sorted and duplicates are marked using the SortSam, BuildBamIndex, and MarkDuplicates functions from Picard ( https://broadinstitute.github.io/picard/ ) (Java v.1.8.0_382 and Picard v.2.25.5). The SplitNCigarReads tool from GATK[ 24 – 26 ] (v.4.2.2.0) is used to split reads containing Ns. To mitigate systematic bias, scVar utilizes GATK BaseRecalibrator and ApplyBQSR to recalibrate base quality, referencing the high-confidence SNPs from the 1000 Genomes Project (1000G) Phase 1 (1000G_phase1.snps.high_confidence.hg38.vcf.gz), as well as Mills and 1000G gold standard indels (Mills_and_1000G_gold_standard.indels.hg38.vcf.gz), and dbSNP version 138 (downloaded as part of the GATK bundle). Candidate mutations are identified using GATK HaplotypeCaller. Due to the relatively low data coverage associated with 10X, scVar specifically sets the ploidy parameter to 40 to facilitate the detection of low-frequency mutations. Both the VEP[ 27 ] (v.107) and Annovar[ 28 ] (Perl v.5.34.0) are used to prioritize and annotate the genomic variants and assess their effects. Concurrently, scVar compiles several clinical databases using Vcfanno[ 64 ] (v.0.3.3), including DOCM[ 65 ], CancerHotspot[ 29 , 30 ], CIVIC ( https://civicdb.org/welcome ), COSMIC[ 31 ], CGI[ 66 , 67 ], and ClinVar[ 32 ], and incorporates these annotations into the INFO field of the VCF file. Mutation genotype assignment and filter To investigate mutations at the single-cell level, scVar quantifies reference and alternate reads for individual cells at each mutation site leveraging the pysam package (v.0.23.0). This process utilizes the unique molecular identifier (UMI) and ’CB’ tags embedded in sequencing reads. To ensure robust genotyping, users can configure critical quality thresholds: (1) a minimum read count per barcode to exclude cells with insufficient sequencing coverage, (2) a minimum mapping quality threshold to filter ambiguously aligned reads, and (3) a minimum base quality score to enhance confidence in base calling. To optimize computational efficiency, the mutation file is partitioned by chromosome, with each chromosomal dataset further subdivided into multiple .tbl files based on user-defined thread specifications. These segmented files are processed in parallel using the multiprocessing package, enabling scalable analysis across large mutation datasets. Ultimately, scVar presents the results detailing the counts of reference and alternate reads for all cells at different mutation sites. A mutation is considered present in a cell if at least one alternate read supports it. To reduce noise, scVar retains only mutations that are present in at least two cells. Mutations catalogued in gnomAD[ 68 ] (v.4.1) or ExAC[ 69 ] (v.0.3) with a population allele frequency greater than 1%, as well as those identified within a PON constructed from WGS data of the 1000 Genomes Project ( https://gatk.broadinstitute.org/hc/en-us/articles/360035890631-Panel-of-Normals-PON -), are systematically excluded from analysis. Additionally, users can set a maximum VAF threshold to refine mutations further, excluding variants exceeding a user-defined VAF limit to reduce potential germline contamination. Finally, scVar generates a sparse cell-by-mutation matrix where -1 indicates no coverage, 0 denotes wild-type, and 1 represents mutation. The chunk processing is employed to optimize memory efficiency. For multi-omics integration, results are structured as a MuData object containing transcriptomic data in the rna modality and mutation data with genomic annotations in the mutation modality, enabling comprehensive joint analysis of genomic variants and gene expression profiles. Downstream analysis in scVar Mutation signature, TMB, Entroy and Simpson’s index analysis scVar consolidates mutations within each cell type to comprehensively analyze the mutation processes across diverse cellular populations. Initially, scVar computes the proportion of mutations occurring in each of the 96 trinucleotide contexts and evaluates the similarity of mutational signature derived from COSMIC[ 31 ] (v.3.3.1) for both the overall results and cell type-specific subsets, employing deconstructSigs[ 70 ] (v.1.8.0). To assess tumor burden, scVar incorporates mosdepth[ 71 ] to assess the length of the genomic regions covered by the BAM file, focusing specifically on the autosomes and the X and Y chromosomes. Subsequently, based on the mutation results, scVar quantifies the number of non-synonymous mutations from the mutation results, including Frame_Shift, Splice_Site, Translation_Start_Site, Nonsense_Mutation, Nonstop_Mutation, and Missense_Mutation types. The tumor burden for each patient is subsequently calculated using the following formula. To quantify intertumoral genomic heterogeneity, scVar calculates Shannon entropy and Simpson’s diversity index for each mutation. Entropy( H ) is calculated according to the following formula: where P 0 denotes the probability of no mutation and P 1 denotes the probability of mutation. Simpson’s diversity index( D ) is calculated according to the following formula: where N 0 equals the number of non-mutated cells, N 1 equals the number of mutated cells and N corresponded to the total number of cells, including both mutated and non-mutated populations. Cell type-specific mutations detection To examine the relationship between genetic mutations and cellular functional heterogeneity, scVar detects mutations that are specific to distinct cell types. For each identified mutation, the counts of mutated and non-mutated cells are determined across predefined cell types or subclusters. A two-sided Fisher’s exact test is then performed using the fisher_exact function from the Python package scipy.stats (version 1.15.2) to assess the enrichment of mutations within specific cell populations. Subsequently, p -values are adjusted for multiple comparisons using the Benjamini-Hochberg (BH) method implemented through the multipletests function from the statsmodels package (version 0.14.4). A comparable statistical approach is employed at the gene level: for each mutated gene, counts of mutated and non-mutated cells are stratified by cell type or subcluster, followed by the application of Fisher’s exact test and subsequent multiple testing correction. This approach enables the systematic identification of mutations and genes exhibiting significant specificity to particular cell types, thereby establishing a connection between genetic variation and functional cellular heterogeneity. Cluster cells based on mutations The final output of scVar consists of a binary matrix delineating cells by mutations, which forms the foundation for mutation-driven cell clustering. Two primary feature extraction methodologies are utilized: (1) direct calculation of variance on the binary matrix to identify mutations with high variability across cells, and (2) variance-based feature selection following TF-IDF normalization, intended to mitigate noise arising from the matrix’s pronounced sparsity. The TF-IDF transformation is performed using the TfidfTransformer function from the scikit-learn library (version 1.5.2). Subsequent to feature selection, pairwise Jaccard distances between cells are computed to quantify dissimilarities in their mutational profiles. The hierarchical clustering framework implemented in scVar supports multiple linkage methods, including the unweighted pair-group method with arithmetic mean (UPGMA) and the weighted pair-group method with arithmetic mean (WPGMA), to construct dendrograms. Finally, clustering results are visualized via the pheatmap function from the ComplexHeatmap package (version 4.4.4), which incorporates dendrograms to highlight cell subpopulations distinguished by unique mutational signatures. Oncogenic pathway and GO function enrichment analysis scVar performs comprehensive enrichment analysis for mutated genes both across the entire cell population and within distinct cell types. The analysis focuses on functionally significant mutation types including: Frame_Shift, Splice_Site, Translation_Start_Site, Nonsense_Mutation, Nonstop_Mutation, Missense_Mutation, as well as variants located in intergenic regions (IGR), 5’ untranslated regions (5’UTR), and 3’ untranslated regions (3’UTR). Mutations classified as non-synonymous but lacking definitive functional implications are excluded from the analysis. For the investigation of oncogenic pathways, scVar employs the OncogenicPathways function from the maftools package (version 2.12.0) in R, with the resulting data visualized through heatmaps generated using the seaborn library (version 0.13.2) in Python. GO enrichment analysis is performed utilizing the enrichGO function from the clusterProfiler package (version 4.4.4), referencing the org.Hs.eg.db database. The analysis permits user specification of statistical significance thresholds via adjustable pvalueCutoff and qvalueCutoff parameters, with p -values corrected for multiple comparisons using the BH method. Visualization of enrichment outcomes is achieved through point plots created with matplotlib (version 3.10.1). Pseudotime analysis To investigate the distribution of mutated cells along differentiation trajectories, scVar implements time-series analysis using Monocle[ 72 ] (v2.24.0), with default parameters. Genes are included in the analysis only if they are expressed in a minimum of 10 cells. The setOrderingFilter function is applied to sort cells based on differentially expressed genes (qval < 0.01) among cell subsets. The reduceDimension function is used to generate traces, employing the DDRTree method. Finally, the trace plot is produced using ggplot. The core workflow of scVar is managed through the Snakemake pipeline management system (version 3.13.3), enabling automated and reproducible execution of analytical processes. Subsequent analysis components are developed as standalone scripts, providing flexibility for tailored customization of individual analytical tasks. scVar generates an extensive, interactive HTML report dynamically rendered using Vue.js, which incorporates visualizations, statistical summaries, and relevant metadata. The software is fully containerized via Docker technology, ensuring a platform-independent environment that supports computational reproducibility across various computing systems. Method of simulating data To rigorously assess the performance of variant calling software, we implemented a benchmark using simulated data. Initially, mutation site selection was performed as follows: for each input BAM file, genome-wide coverage profiles were generated using mosdepth[ 71 ]. Candidate insertion sites (n=5000 per sample) were then identified from genomic regions exhibiting sequencing depths of at least 10X, utilizing a multi-step filtering strategy implemented via pysam. This filtering procedure involved: (a) exclusion of population variants catalogued in gnomAD[ 68 ] (v.4.1) or ExAC[ 69 ] (v.0.3) with a population frequency exceeding 1%, as well as variants present in a PON generated from WGS data of the 1000 Genomes Project, (b) retention of sites with canonical reference nucleotides (A/T/C/G), and (c) enforcement of stringent quality criteria, specifically mapping quality ≥200 and base quality ≥35. The distribution of mutations across chromosomes was proportionally allocated based on chromosome length, and multiprocessing techniques were employed to optimize computational efficiency. Subsequently, BAM file modification was conducted by introducing controlled mutagenesis at the selected candidate sites. Initially, all reads were reverted to their reference-aligned states using pysam. Site-specific variants were then incorporated such that the allele fractions of the modified reads precisely matched predetermined VAFs of 10%, 20%, 30%, 40%, and 50%, resulting in five simulated samples per input file. The performance of four variant callers—cellSNP, Monopogen, SComatic, and scVar—was assessed using the F1-score metric. Precision and recall were systematically computed by comparing the ground-truth mutations with the variants identified by each caller across all VAF levels. Valuation with WES To validate the performance of the scVar software, we leveraged the paired WES and scRNA-seq datasets, employing mutations identified from the WES data as the reference gold standard. The WES data were aligned to the human reference genome using BWA MEM[ 73 ] (v.0.7.17). Subsequent processing of the resulting BAM files included duplicate removal and base quality score recalibration in accordance with the GATK Best Practices workflow. SNVs from WES data were identified using Strelka2[ 74 ] with default parameters, incorporating matched normal samples as germline controls. Genomic regions were then filtered to retain those with a minimum sequencing depth of 1X in the WES data and at least 10X in the scRNA-seq data, as determined by mosdepth[ 71 ]. To define a set of gold standard mutation candidates, overlapping regions between these datasets were determined via the Intersect function of Bedtools[ 75 ]. Regenotyping was performed on candidate mutations present in both the WES and scRNA-seq BAM files. Gold standard mutations were defined based on the following criteria: absence of alternative allele reads in the WES normal sample, presence of at least one alternative allele reads in the WES tumor sample, and more than two alternative allele reads in the scRNA-seq tumor sample. For performance evaluation, true positives (TP) were defined as mutations detected by scVar that were included in the established true mutation set. False negatives (FN) comprised mutations within the true set that scVar failed to identify. False positives (FP) were mutations uniquely reported by scVar but subsequently determined to be absent upon regenotyping in both WES and scRNA-seq BAM files. Methods of analyzing real data Single-cell gene expression quantification and determination of the major cell types Raw gene expression matrices were generated for each sample using Cell Ranger Count and SoupX, and merged using Scanpy[ 76 ] (v.1.10.2) concat method. Cells were filtered to retain those expressing more than 10 and fewer than 15,000 genes, with ribosomal transcript content below 50% and mitochondrial transcript content below 10%. Normalization was conducted by scaling the total counts per cell to 10,000 transcripts (normalize_total(target_sum=1e4)) followed by log(1+x) transformation. Highly variable genes (HVGs) were identified using Seurat’s variance-stabilizing method (highly_variable_genes(n_top_genes=2000, flavor=’seurat’)), selecting the top 2,000 genes for downstream analyses. Prior to PCA, data were scaled with a maximum value of 10 (scale(max_value=10)). PCA was conducted with default parameters, retaining the first 50 principal components. A k-nearest neighbor graph was constructed using 10 neighbors and 50 principal components (neighbors(n_neighbors=10, n_pcs=50)), followed by UMAP embedding with a minimum distance parameter of 1 (umap(min_dist=1)). Cell clustering was performed using the Leiden algorithm at a resolution of 1, implemented via Leidenalg (leiden(resolution=1, random_state=0, flavor="leidenalg", key_added="leiden_res1")). UMAP visualizations were generated using the utils.embedding function from the Omicverse[ 77 ] package (v.1.6.7). Cell type annotation was conducted by integrating marker genes from the CellMarker database[ 78 ] (CellMarker Augmented (2021)) and markers reported in the literature. For annotation based on database-derived marker genes, pathway enrichment analysis was carried out using the ov.single.pathway_aucell_enrichment function from the Omicverse package, which mapped gene sets from a predefined pathway dictionary (pathways_dict) to cellular activities utilizing eight parallel workers (num_workers=8). Differential gene expression analysis was subsequently executed using sc.tl.rank_genes_groups with the Wilcoxon rank-sum test (method="wilcoxon"), grouping cells according to Leiden clustering results. The top three differentially expressed genes per cluster were visualized as a dot plot through sc.pl.rank_genes_groups_dotplot, generated by variance-scaled expression values (standard_scale=’var’). Subclustering and trajectories analysis of alveolar cells To delineate subclusters within alveolar cells, we conducted a focused reanalysis of cells classified as alveolar. Initially, raw count data specific to alveolar cells were extracted. Subsequent data preprocessing involved normalization and logarithmic transformation. HVGs were identified employing the Seurat variance-stabilizing approach (sc.pp.highly_variable_genes) with the following parameters: min_mean = 0.0125, max_mean = 3, min_disp = 0.5, and flavor set to ’seurat’. PCA was then performed using default settings. A k-nearest neighbor graph was constructed utilizing 10 neighbors and 40 principal components (sc.pp.neighbors(n_neighbors=10, n_pcs=40)). UMAP embedding was subsequently applied with a minimum distance parameter of 0.5 (sc.tl.umap(min_dist=0.5)). To identify transcriptionally distinct subpopulations, Leiden clustering was implemented at a resolution of 0.1 (sc.tl.leiden). Finally, pseudotime analysis was conducted across all alveolar cells to infer dynamic cellular trajectories. To investigate pathway activity within alveolar subclusters, we compared the expression patterns of genes associated with the RAF– MEK – ERK cascade and the PI3K – AKT – mTOR signaling pathway. Visualization of gene expression level was performed using violin plots generated with the seaborn (sns.violinplot) package. Differential expression analysis was conducted between cells harboring the KRAS G12A mutation and those without detectable KRAS G12A mutation, using the sc.tl.rank_genes_groups function in Scanpy with the t-test method. From the top 20 genes ranked in this analysis, we concentrated on those implicated in tumorigenesis and cancer progression, particularly oncogenes and driver genes previously documented in lung cancer and other malignancies. Notable genes identified included KRTCAP2, ATP1A1, RPN2, NMT1, REEP3, IFITM1, ARF5 and BPTF. Availability of data and materials scVar is available on GitHub https://github.com/WeiZhaoLuck/scVar and Biocode https://ngdc.cncb.ac.cn/biocode/tool/BT008000 at National Genomics Data Center [ 79 ]. The code is released under the version 3 of the GNU General Public License. Matched WES and scRNA-seq datasets from eight cutaneous squamous cell carcinoma (cSCC) tumors were retrieved from the NCBI Gene Expression Omnibus database under accession numbers GSE144237 ( https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144237 ) and GSE144240 ( https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144240 ). Public single-cell lung cancer data were downloaded from ArrayExpress under accession E-MTAB-6149 ( https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-6149/files ). Funding Genomics Data Center Operation and Maintenance of Chinese Academy of Sciences [CAS-WX2022SDC-XK05]; the Alliance of International Science Organizations [Grant No. ANSO-PA-2023-07]; the Open Biodiversity and Health Big Data Programme of IUBS. Contributions Y.B., X.Z., and W.Z. designed the project. W.Z. and X.Z built the software. W.Z. conducted the software benchmark, application, and analysis. Y.B., X.Z., and W.Z. wrote the manuscript. All authors contributed to editing the manuscript. All authors read and approved the final manuscript. Ethics declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Competing interests Not applicable. Additional information Supplementary Information Supplementary Note To rigorously assess the efficacy of mutation detection tools in identifying low-frequency variants from 10X scRNA-seq data, we developed a simulation framework grounded in real scRNA-seq datasets. This framework involved the computational generation of hybrid tumor samples by mixing sequencing reads from two distinct individuals (referred to as sample1 and sample2) at predetermined ratios (10%:90%, 20%:80%, 30%:70%, 40%:60%, and 50%:50%). The simulation process comprised four principal stages. Initially, genomic regions exhibiting overlapping coverage between the two samples were identified using Bedtools (employing the bamtobed, merge, and intersect functions). Subsequently, reads corresponding to these shared regions were extracted from each BAM file utilizing Samtools. To simulate varying mutation frequencies, differential downsampling of the original BAM files was performed via Picard’s DownsampleSam function, retaining 10–50% of reads for sample1 and 90–50% for sample2, followed by merging with Samtools. For the establishment of a gold-standard mutation set, a stringent filtering pipeline was applied: (1) initial variant calling was conducted using Strelka2 (version 2.9.10), requiring variants to exhibit 100% VAF and coverage exceeding five reads in both samples; (2) homozygous germline mutations specific to sample1 were selected as candidate variants; (3) regenotyping validation across all samples was performed to retain only mutations demonstrating complete allele segregation (i.e., 100% alternative alleles in sample1 and 100% reference alleles in sample2); and (4) final VAF thresholding was applied, allowing a ±2% deviation from the expected downsampling ratio to define the benchmark mutation set. We conducted a systematic comparison of six mutation detection methodologies, including the GATK Best Practices workflow (incorporating Picard MarkDuplicates, SplitNCigarReads, BaseRecalibrator, ApplyBQSR, and HaplotypeCaller with ploidy set to 40), Strelka2 (v2.9.10), and VarScan (v2.3.9). Recall rates were computed as the ratio of TP to the sum of true positives and false negatives (TP + FN), where true positives corresponded to gold-standard variants detected by each tool, and false negatives represented benchmark mutations that were not identified. Notably, the optimized GATK pipeline exhibited superior performance across all mixing ratios, achieving recall rates ranging from 84.5% to 96.3%, in contrast to 62.1% to 89.4% observed for alternative tools (Supplementary Figure 1a). This performance advantage persisted when mutations were stratified by substitution type (C>A, C>G, C>T, T>A, T>C, T>G), with GATK and Strelka2 demonstrating comparable detection rates across all categories (Supplementary Figure 1b). The demonstrated robustness of GATK, particularly at low mutation frequencies (10–30% mixtures), informed its selection as the primary variant caller within the scVar analytical framework. Acknowledgements We would like to thank Yingke Ma, Song Wu, Tong Jin, Mochen Zhang, and anonymous reviewers for helpful feedback and discussions. Funder Information Declared Genomics Data Center Operation and Maintenance of Chinese Academy of Sciences , CAS-WX2022SDC-XK05 the Alliance of International Science Organizations , ANSO-PA-2023-07 the Open Biodiversity and Health Big Data Programme of IUBS Footnotes https://github.com/WeiZhaoLuck/scVar References 1. ↵ Dagogo-Jack I , Shaw AT . Tumour heterogeneity and resistance to cancer therapies . Nat Rev Clin Oncol . 2018 ; 15 ( 2 ): 81 – 94 . OpenUrl CrossRef PubMed 2. ↵ Grzywa TM , Paskal W , Włodarski PK . Intratumor and Intertumor Heterogeneity in Melanoma . Translational Oncology . 2017 ; 10 ( 6 ): 956 – 975 . OpenUrl PubMed 3. ↵ Liu J , Dang H , Wang XW . The significance of intertumor and intratumor heterogeneity in liver cancer . Exp Mol Med . 2018 ; 50 ( 1 ): e416 – e416 . OpenUrl CrossRef PubMed 4. ↵ Morris LGT , Riaz N , Desrichard A , Şenbabaoğlu Y , Hakimi AA , Makarov V , et al. Pan-cancer analysis of intratumor heterogeneity as a prognostic determinant of survival . Oncotarget . 2016 ; 7 ( 9 ): 10051 – 10063 . OpenUrl CrossRef PubMed 5. ↵ Lim ZF , Ma PC . Emerging insights of tumor heterogeneity and drug resistance mechanisms in lung cancer targeted therapy . J Hematol Oncol . 2019 ; 12 ( 1 ): 134 . OpenUrl CrossRef PubMed 6. Goyette MA , Lipsyc-Sharf M , Polyak K . Clinical and translational relevance of intratumor heterogeneity . Trends Cancer . 2023 ; 9 ( 9 ): 726 – 737 . OpenUrl PubMed 7. ↵ Safri F , Nguyen R , Zerehpooshnesfchi S , George J , Qiao L . Heterogeneity of hepatocellular carcinoma: from mechanisms to clinical implications . Cancer Gene Ther . 2024 ; 31 ( 8 ): 1105 – 1112 . OpenUrl CrossRef PubMed 8. ↵ Landau DA , Carter SL , Stojanov P , McKenna A , Stevenson K , Lawrence MS , et al. Evolution and impact of subclonal mutations in chronic lymphocytic leukemia . Cell . 2013 ; 152 ( 4 ): 714 – 726 . OpenUrl CrossRef PubMed Web of Science 9. Suzuki H , Aoki K , Chiba K , Sato Y , Shiozawa Y , Shiraishi Y , et al. Mutational landscape and clonal architecture in grade II and III gliomas . Nat Genet . 2015 ; 47 ( 5 ): 458 – 468 . OpenUrl CrossRef PubMed 10. ↵ Dentro SC , Leshchiner I , Haase K , Tarabichi M , Wintersinger J , Deshwar AG , et al. Characterizing genetic intra-tumor heterogeneity across 2,658 human cancer genomes . Cell . 2021 ; 184 ( 8 ): 2239 – 2254.e39 . OpenUrl CrossRef PubMed 11. ↵ Wang S , Sun ST , Zhang XY , Ding HR , Yuan Y , He JJ , et al. The Evolution of Single-Cell RNA Sequencing Technology and Application: Progress and Perspectives . IJMS . 2023 ; 24 ( 3 ): 2943 . OpenUrl PubMed 12. ↵ Yamawaki TM , Lu DR , Ellwanger DC , Bhatt D , Manzanillo P , Arias V , et al. Systematic comparison of high-throughput single-cell RNA-seq methods for immune cell profiling . BMC Genomics . 2021 ; 22 ( 1 ): 66 . OpenUrl CrossRef PubMed 13. Wang X , He Y , Zhang Q , Ren X , Zhang Z . Direct Comparative Analyses of 10X Genomics Chromium and Smart-Seq2 . Genomics, Proteomics & Bioinformatics . 2021 ; 19 ( 2 ): 253 – 266 . OpenUrl CrossRef PubMed 14. ↵ Jovic D , Liang X , Zeng H , Lin L , Xu F , Luo Y . Single-cell RNA sequencing technologies and applications: A brief overview . Clinical & Translational Med . 2022 ; 12 ( 3 ): e694 . OpenUrl CrossRef 15. ↵ Wang Y , Flowers CR , Wang M , Huang X , Li Z . CASi: A framework for cross-timepoint analysis of single-cell RNA sequencing data . Sci Rep . 2024 ; 14 ( 1 ): 10633 . OpenUrl PubMed 16. ↵ Jarczak J , Bujko K , Ratajczak MZ , Kucia M . scRNA-seq revealed transcriptional signatures of human umbilical cord primitive stem cells and their germ lineage origin regulated by imprinted genes . Sci Rep . 2024 ; 14 ( 1 ): 29264 . OpenUrl PubMed 17. ↵ Zhang Q , He Y , Luo N , Patel SJ , Han Y , Gao R , et al. Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma . Cell . 2019 ; 179 ( 4 ): 829 – 845.e20 . OpenUrl CrossRef PubMed 18. ↵ Bischoff P , Trinks A , Obermayer B , Pett JP , Wiederspahn J , Uhlitz F , et al. Single-cell RNA sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma . Oncogene . 2021 ; 40 ( 50 ): 6748 – 6758 . OpenUrl CrossRef PubMed 19. ↵ Zheng GXY , Terry JM , Belgrader P , Ryvkin P , Bent ZW , Wilson R , et al. Massively parallel digital transcriptional profiling of single cells . Nat Commun . 2017 ; 8 ( 1 ): 14049 . OpenUrl CrossRef PubMed 20. ↵ N. M. P , Liu H , Dillard C , Ibeawuchi H , Alsaeedy T , Chan H , et al. Improved SNV Discovery in Barcode-Stratified scRNA-seq Alignments . Genes . 2021 ; 12 ( 10 ): 1558 . OpenUrl 21. ↵ Huang X , Huang Y . Cellsnp-lite: an efficient tool for genotyping single cells . Bioinformatics . 2021 ; 37 ( 23 ): 4569 – 4571 . OpenUrl CrossRef PubMed 22. ↵ Dou J , Tan Y , Kock KH , Wang J , Cheng X , Tan LM , et al. Single-nucleotide variant calling in single-cell sequencing data with Monopogen . Nat Biotechnol . 2024 ; 42 ( 5 ): 803 – 812 . OpenUrl CrossRef PubMed 23. ↵ Muyas F , Sauer CM , Valle-Inclán JE , Li R , Rahbari R , Mitchell TJ , et al. De novo detection of somatic mutations in high-throughput single-cell profiling data sets . Nat Biotechnol . 2023 ; DOI: 10.1038/s41587-023-01863-z . OpenUrl CrossRef 24. ↵ McKenna A , Hanna M , Banks E , Sivachenko A , Cibulskis K , Kernytsky A , et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data . Genome Res . 2010 ; 20 ( 9 ): 1297 – 1303 . OpenUrl Abstract / FREE Full Text 25. DePristo MA , Banks E , Poplin R , Garimella KV , Maguire JR , Hartl C , et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data . Nat Genet . 2011 ; 43 ( 5 ): 491 – 498 . OpenUrl CrossRef PubMed Web of Science 26. ↵ Van Der Auwera GA , Carneiro MO , Hartl C , Poplin R , Del Angel G , Levy-Moonshine A , et al. From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline . CP in Bioinformatics . 2013 ; 43 ( 1 ). 27. ↵ McLaren W , Gil L , Hunt SE , Riat HS , Ritchie GRS , Thormann A , et al. The Ensembl Variant Effect Predictor . Genome Biol . 2016 ; 17 ( 1 ): 122 . OpenUrl CrossRef PubMed 28. ↵ Wang K , Li M , Hakonarson H . ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data . Nucleic Acids Research . 2010 ; 38 ( 16 ): e164 – e164 . OpenUrl CrossRef PubMed 29. ↵ Chang MT , Asthana S , Gao SP , Lee BH , Chapman JS , Kandoth C , et al. Identifying recurrent mutations in cancer reveals widespread lineage diversity and mutational specificity . Nat Biotechnol . 2016 ; 34 ( 2 ): 155 – 163 . OpenUrl CrossRef PubMed 30. ↵ Chang MT , Bhattarai TS , Schram AM , Bielski CM , Donoghue MTA , Jonsson P , et al. Accelerating Discovery of Functional Mutant Alleles in Cancer . Cancer Discovery . 2018 ; 8 ( 2 ): 174 – 183 . OpenUrl Abstract / FREE Full Text 31. ↵ Sondka Z , Dhir NB , Carvalho-Silva D , Jupe S , Madhumita , McLaren K , et al. COSMIC: a curated database of somatic variants and clinical data for cancer . Nucleic Acids Research . 2024 ; 52 ( D1 ): D1210 – D1217 . OpenUrl CrossRef PubMed 32. ↵ Landrum MJ , Lee JM , Benson M , Brown GR , Chao C , Chitipiralla S , et al. ClinVar: improving access to variant interpretations and supporting evidence . Nucleic Acids Research . 2018 ; 46 ( D1 ): D1062 – D1067 . OpenUrl CrossRef PubMed 33. ↵ Ji AL , Rubin AJ , Thrane K , Jiang S , Reynolds DL , Meyers RM , et al. Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma . Cell . 2020 ; 182 ( 2 ): 497 – 514.e22 . OpenUrl CrossRef PubMed 34. ↵ Lambrechts D , Wauters E , Boeckx B , Aibar S , Nittner D , Burton O , et al. Phenotype molding of stromal cells in the lung tumor microenvironment . Nat Med . 2018 ; 24 ( 8 ): 1277 – 1289 . OpenUrl CrossRef PubMed 35. ↵ Karin M , Greten FR . NF-κB: linking inflammation and immunity to cancer development and progression . Nat Rev Immunol . 2005 ; 5 ( 10 ): 749 – 759 . OpenUrl CrossRef PubMed Web of Science 36. ↵ Chen W . NF-kappaB in lung cancer, a carcinogenesis mediator and a prevention and therapy target . Front Biosci . 2011 ; 16 ( 1 ): 1172 . OpenUrl CrossRef PubMed 37. ↵ Song X , Li X , Ge Y , Song J , Wei Q , He M , et al. Alternative splicing events and function in the tumor microenvironment: New opportunities and challenges . International Immunopharmacology . 2023 ; 123 : 110718 . OpenUrl PubMed 38. ↵ Wu YC , Tang SJ , Sun GH , Sun KH . CXCR7 mediates TGFβ1-promoted EMT and tumor-initiating features in lung cancer . Oncogene . 2016 ; 35 ( 16 ): 2123 – 2132 . OpenUrl CrossRef PubMed 39. ↵ Li X , Cao X , Zhao H , Guo M , Fang X , Li K , et al. Hypoxia Activates Notch4 via ERK/JNK/P38 MAPK Signaling Pathways to Promote Lung Adenocarcinoma Progression and Metastasis . Front Cell Dev Biol . 2021 ; 9 : 780121 . OpenUrl PubMed 40. ↵ Xia H , Green DR , Zou W . Autophagy in tumour immunity and therapy . Nat Rev Cancer . 2021 ; 21 ( 5 ): 281 – 297 . OpenUrl CrossRef PubMed 41. ↵ Giansanti M , Theinert T , Boeing SK , Haas D , Schlegel PG , Vacca P , et al. Exploiting autophagy balance in T and NK cells as a new strategy to implement adoptive cell therapies . Mol Cancer . 2023 ; 22 ( 1 ): 201 . OpenUrl PubMed 42. ↵ Engelman JA , Cantley LC . The Role of the ErbB Family Members in Non–Small Cell Lung Cancers Sensitive to Epidermal Growth Factor Receptor Kinase Inhibitors . Clinical Cancer Research . 2006 ; 12 ( 14 ): 4372s – 4376s . OpenUrl Abstract / FREE Full Text 43. ↵ Whitsett JA , Wert SE , Weaver TE . Alveolar Surfactant Homeostasis and the Pathogenesis of Pulmonary Disease . Annu Rev Med . 2010 ; 61 ( 1 ): 105 – 119 . OpenUrl CrossRef PubMed Web of Science 44. ↵ Di Lorenzo A , Bolli E , Tarone L , Cavallo F , Conti L . Toll-Like Receptor 2 at the Crossroad between Cancer Cells, the Immune System, and the Microbiota . Int J Mol Sci . 2020 ; 21 ( 24 ): 9418 . OpenUrl PubMed 45. ↵ Zhang T , Ma C , Zhang Z , Zhang H , Hu H . NF-κB signaling in inflammation and cancer . MedComm . 2021 ; 2 ( 4 ): 618 – 653 . OpenUrl PubMed 46. ↵ Gergen AK , Kohtz PD , Halpern AL , Li A , Meng X , Reece TB , et al. Activation of Toll-Like Receptor 2 Promotes Proliferation of Human Lung Adenocarcinoma Cells . Anticancer Res . 2020 ; 40 ( 10 ): 5361 – 5369 . OpenUrl Abstract / FREE Full Text 47. ↵ Huang L , Guo Z , Wang F , Fu L . KRAS mutation: from undruggable to druggable in cancer . Sig Transduct Target Ther . 2021 ; 6 ( 1 ): 386 . OpenUrl 48. ↵ O’Sullivan É , Keogh A , Henderson B , Finn SP , Gray SG , Gately K . Treatment Strategies for KRAS-Mutated Non-Small-Cell Lung Cancer . Cancers (Basel ). 2023 ; 15 ( 6 ): 1635 . OpenUrl PubMed 49. ↵ Luo J . KRAS mutation in pancreatic cancer . Semin Oncol . 2021 ; 48 ( 1 ): 10 – 18 . OpenUrl CrossRef PubMed 50. ↵ Lee JK , Sivakumar S , Schrock AB , Madison R , Fabrizio D , Gjoerup O , et al. Comprehensive pan-cancer genomic landscape of KRAS altered cancers and real-world outcomes in solid tumors . npj Precis Onc . 2022 ; 6 ( 1 ): 91 . OpenUrl 51. ↵ Fiala O , Buchler T , Mohelnikova-Duchonova B , Melichar B , Matejka VM , Holubec L , et al. G12V and G12A KRAS mutations are associated with poor outcome in patients with metastatic colorectal cancer treated with bevacizumab . Tumour Biol . 2016 ; 37 ( 5 ): 6823 – 6830 . OpenUrl PubMed 52. ↵ Han G , Sinjab A , Rahal Z , Lynch AM , Treekitkarnmongkol W , Liu Y , et al. An atlas of epithelial cell states and plasticity in lung adenocarcinoma . Nature . 2024 ; 627 ( 8004 ): 656 – 663 . OpenUrl CrossRef PubMed 53. ↵ Dai M , Lu JJ , Guo W , Yu W , Wang Q , Tang R , et al. BPTF promotes tumor growth and predicts poor prognosis in lung adenocarcinomas . Oncotarget . 2015 ; 6 ( 32 ): 33878 – 33892 . OpenUrl CrossRef PubMed 54. ↵ Sakamoto S , Inoue H , Kohda Y , Ohba S ichi , Mizutani T , Kawada M. Interferon-Induced Transmembrane Protein 1 (IFITM1) Promotes Distant Metastasis of Small Cell Lung Cancer . IJMS . 2020 ; 21 ( 14 ): 4934 . OpenUrl PubMed 55. ↵ Sun P , Zhang H , Shi J , Xu M , Cheng T , Lu B , et al. KRTCAP2 as an immunological and prognostic biomarker of hepatocellular carcinoma . Colloids and Surfaces B: Biointerfaces . 2023 ; 222 : 113124 . OpenUrl 56. ↵ Han Z , Wang Y , Han L , Yang C . RPN2 in cancer: An overview . Gene . 2023 ; 857 : 147168 . OpenUrl PubMed 57. ↵ Liu GH , Tan XY , Xu ZY , Li JX , Zhong GH , Zhai JW , et al. REEP3 is a potential diagnostic and prognostic biomarker correlated with immune infiltration in pancreatic cancer . Sci Rep . 2024 ; 14 ( 1 ). 58. ↵ Li Q , Li F , Song X , Lu N , Jing X , Wen H , et al. Pan-cancer analysis of ARFs family and ARF5 promoted the progression of hepatocellular carcinoma . Heliyon . 2024 ; 10 ( 7 ): e29099 . OpenUrl 59. ↵ Zheng GXY , Terry JM , Belgrader P , Ryvkin P , Bent ZW , Wilson R , et al. Massively parallel digital transcriptional profiling of single cells . Nat Commun . 2017 ; 8 ( 1 ): 14049 . OpenUrl CrossRef PubMed 60. ↵ Young MD , Behjati S . SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data . GigaScience . 2020 ; 9 ( 12 ): giaa151 . OpenUrl CrossRef PubMed 61. ↵ Hao Y , Hao S , Andersen-Nissen E , Mauck WM , Zheng S , Butler A , et al. Integrated analysis of multimodal single-cell data . Cell . 2021 ; 184 ( 13 ): 3573 – 3587.e29 . OpenUrl CrossRef PubMed 62. ↵ Aran D , Looney AP , Liu L , Wu E , Fong V , Hsu A , et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage . Nat Immunol . 2019 ; 20 ( 2 ): 163 – 172 . OpenUrl CrossRef PubMed 63. ↵ Barnett DW , Garrison EK , Quinlan AR , Strömberg MP , Marth GT . BamTools: a C++ API and toolkit for analyzing and managing BAM files . Bioinformatics . 2011 ; 27 ( 12 ): 1691 – 1692 . OpenUrl CrossRef PubMed Web of Science 64. ↵ Pedersen BS , Layer RM , Quinlan AR . Vcfanno: fast, flexible annotation of genetic variants . Genome Biol . 2016 ; 17 ( 1 ): 118 . OpenUrl CrossRef PubMed 65. ↵ Ainscough BJ , Griffith M , Coffman AC , Wagner AH , Kunisaki J , Choudhary MN , et al. DoCM: a database of curated mutations in cancer . Nat Methods . 2016 ; 13 ( 10 ): 806 – 807 . OpenUrl CrossRef PubMed 66. ↵ Tamborero D , Rubio-Perez C , Deu-Pons J , Schroeder MP , Vivancos A , Rovira A , et al. Cancer Genome Interpreter annotates the biological and clinical relevance of tumor alterations . Cancer Biology , 2017 ; 67. ↵ Muiños F , Martinez-Jimenez F , Pich O , Gonzalez-Perez A , Lopez-Bigas N . In silico saturation mutagenesis of cancer genes . Cancer Biology , 2020 ; 68. ↵ Karczewski KJ , Francioli LC , Tiao G , Cummings BB , Alföldi J , Wang Q , et al. The mutational constraint spectrum quantified from variation in 141,456 humans . Nature . 2020 ; 581 (7809): 434 – 443 . OpenUrl CrossRef PubMed 69. ↵ Karczewski KJ , Weisburd B , Thomas B , Solomonson M , Ruderfer DM , Kavanagh D , et al. The ExAC browser: displaying reference data information from over 60 000 exomes . Nucleic Acids Res . 2017 ; 45 ( D1 ): D840 – D845 . OpenUrl CrossRef PubMed 70. ↵ Rosenthal R , McGranahan N , Herrero J , Taylor BS , Swanton C . deconstructSigs: delineating mutational processes in single tumors distinguishes DNA repair deficiencies and patterns of carcinoma evolution . Genome Biol . 2016 ; 17 ( 1 ): 31 . OpenUrl CrossRef PubMed 71. ↵ Pedersen BS , Quinlan AR . Mosdepth: quick coverage calculation for genomes and exomes . Bioinformatics . 2018 ; 34 ( 5 ): 867 – 868 . OpenUrl CrossRef PubMed 72. ↵ Qiu X , Mao Q , Tang Y , Wang L , Chawla R , Pliner HA , et al. Reversed graph embedding resolves complex single-cell trajectories . Nat Methods . 2017 ; 14 ( 10 ): 979 – 982 . OpenUrl CrossRef PubMed 73. ↵ Jung Y , Han D . BWA-MEME: BWA-MEM emulated with a machine learning approach . Bioinformatics . 2022 ; 38 ( 9 ): 2404 – 2413 . OpenUrl CrossRef PubMed 74. ↵ Kim S , Scheffler K , Halpern AL , Bekritsky MA , Noh E , Källberg M , et al. Strelka2: fast and accurate calling of germline and somatic variants . Nat Methods . 2018 ; 15 ( 8 ): 591 – 594 . OpenUrl CrossRef PubMed 75. ↵ Quinlan AR , Hall IM . BEDTools: a flexible suite of utilities for comparing genomic features . Bioinformatics . 2010 ; 26 ( 6 ): 841 – 842 . OpenUrl CrossRef PubMed Web of Science 76. ↵ Wolf FA , Angerer P , Theis FJ . SCANPY: large-scale single-cell gene expression data analysis . Genome Biol . 2018 ; 19 ( 1 ): 15 . OpenUrl CrossRef PubMed 77. ↵ Zeng Z , Ma Y , Hu L , Tan B , Liu P , Wang Y , et al. OmicVerse: a framework for bridging and deepening insights across bulk and single-cell sequencing . Nat Commun . 2024 ; 15 ( 1 ): 5983 . OpenUrl CrossRef PubMed 78. ↵ Hu C , Li T , Xu Y , Zhang X , Li F , Bai J , et al. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data . Nucleic Acids Research . 2023 ; 51 ( D1 ): D870 – D876 . OpenUrl CrossRef PubMed 79. ↵ CNCB-NGDC Members and Partners , Bao Y , Bai X , Bu C , Chen H , Chen H , et al. Database Resources of the National Genomics Data Center, China National Center for Bioinformation in 2025 . Nucleic Acids Research . 2025 ; 53 ( D1 ): D30 – D44 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted October 14, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Joint analysis of expression and variation at single cell resolution by scVar Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Joint analysis of expression and variation at single cell resolution by scVar Wei Zhao , Fei Yang , Wenting Zong , Xinchang Zheng , Yiming Bao bioRxiv 2025.10.12.681753; doi: https://doi.org/10.1101/2025.10.12.681753 Share This Article: Copy Citation Tools Joint analysis of expression and variation at single cell resolution by scVar Wei Zhao , Fei Yang , Wenting Zong , Xinchang Zheng , Yiming Bao bioRxiv 2025.10.12.681753; doi: https://doi.org/10.1101/2025.10.12.681753 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7633) Biochemistry (17680) Bioengineering (13888) Bioinformatics (41927) Biophysics (21445) Cancer Biology (18585) Cell Biology (25491) Clinical Trials (138) Developmental Biology (13373) Ecology (19897) Epidemiology (2067) Evolutionary Biology (24308) Genetics (15606) Genomics (22494) Immunology (17736) Microbiology (40385) Molecular Biology (17175) Neuroscience (88583) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4822) Physiology (7641) Plant Biology (15149) Scientific Communication and Education (2045) Synthetic Biology (4293) Systems Biology (9822) Zoology (2271)