cells were then centrifuged for 5 min at 300 × g, 4 °c and the pellet was resuspended into 1× rbc lysis buffer ( ebioscience ) for 2 min at room temperature, after which 20 ml of cold pbs were added to stop the lysis reaction. cells were cryopreserved in 5 % dmso in knockout serum replacement ( kosr; gibco, catalogue number: 10828010 ) until further use. on the day of facs sorting, cells were rapidly thawed at 37 °c and transferred to complete rpmi media. live - cell enrichment was performed using macs dead cell removal kit ( miltenyi biotec ) following the manufacturer ’ s instructions. red blood cells were further depleted by negative selection using cd235a microbeads ( miltenyi biotec ) and macs ls columns ( miltenyi biotec ), following the manufacturer ’ s instructions. for facs sorting, cells were stained with zombie aqua to exclude dead cells and the cocktail of antibod ies for 30 min at 4 °c. cells were centrifuged for 5 min at 300 × g, 4 °c, resuspended in 500 μl of 5 % fbs in pbs and subsequently filtered into polypropylene facs tubes. immune cells were sorted as live, cd45 +; mdsc were sorted as live, cd45 +, lineage - ( lin: cd3, cd56, cd19 ), cd33 +, hla - dr - / low ( supplementary data 22 and supplementary fig. 1a ). cells were sorted into a 1. 5 - ml tube, counted and submitted for 10x scrna - seq library preparation. each cell suspension was submitted for 3 ’ single - cell rna sequencing using single cell g chip kit, chemistry v3. 1 ( 10x genomics pleasanton, ca, usa ), following the manufacturer ’ s instructions. libraries were sequenced on an illumina novaseq 6000, and mapped to the grch38 human reference genome using the cellranger toolkit ( version 3. 1. 0 ). integrating numerous samples, notably from diverse cancer subtypes and adjacent normal tissues, is challenging due to variations in gene programmes between samples. consequently, these differences often hinder a coherent biological alignment when attempting simultaneous embedding. most current integration techniques, primarily focused on batch correction, operate under the assumption of shared cell states across samples. however, while they aim to mitigate technical disparities, they might inadvertently erase genuine biological distinc tions. therefore, we applied the qc filtration and doublet removal on the merged dataset ( tumour + b / h ) but we split the datasets between tumour and b / h for hvg selection, pca, batch correction ( using harmony ), clustering and annotations. starting from the unnormalised, uncorrected gene expression matrices produced ( per sample ) by the cellranger protocol, we performed careful downstream analysis of the scrna - seq data. for each cellranger output ( corresponding to a specific technical and biological replicate of the separate tumour, background and healthy data ) we identified low - quality cells or empty droplets by applying the barcoderanks and emptydrops functions using the r package dropletutils. following per - sample droplets removal, the complete set of cell expression matrices was merged ( we merged tumour, background, and healthy samples ), and quality control ( qc ) was applied to the resultant merged matrix. the remaining analysis is implemented using standard approaches in the scanpy framework. the qc is based on three parameters: the total umi count ( lower - upper threshold ), the number of detected genes ( lower - upper threshold ), and the proportion of mitochondrial gene count per cell ( 20 % fraction upper bound ). we applied scrublet to remove potential doublets with 0. 06 as the expected doublet rate and then filtered the results using the parameter values ( 2 for minimum read count of cell, 3 for minimum detected cell of gene, 85 for minimum gene variability percentage, and 30 for the number of principal components used to embed the transcriptomes prior to k - nearest - neighbour graph construction ). the resulting merged and filtered expression matrix is then normalised using the scaling factor 10, 000, followed by log1p transformation. for dimensionality reduction, we first selected sets of highly - variable genes ( hvgs ) from the initial gene set of 25, 718. starting from the hvg selection, the merged matrix was split into two separate matrices: tumour, and combined background / healthy which we refer to as b / h. after hvg selection, 1604 genes were selected from the tumour matrix and 1486 from b / h. from these separate hvg sets, we applied dimensionality reduction using principal component analysis ( pca ). next, we performed pca separately for tumour and for b / h and retained the top 15 components, according to the scree plot elbow rule. the resulting matrix is then batch corrected to account for additional technical variations arising between samples which are non - biological in origin. we apply batch correction by using harmonypy ( a python version of the original harmonyr package ), based on recommended benchmarking against other procedures. following between - sample batch correction, we computed a neighbourhood graph and applied leiden clustering ( with leiden resolution being 1 ) to the 15 - dimensional harmonised pca space. for visualisation purposes, we used uniform manifold approximation and projection ( umap ) manifold embedding to capture the global features of the 15 - dimensional clustered manifold and represent the global structure in two and three dimensions. we identified top 100 representative genes for each cluster by performing the wilcoxon signed - rank test with the bonferroni correction, followed by a filtering to obtain genes overexpressed in the target group ( minimum log fold change as 0 ) and expressed in at least 30 % of cells within the group. we did not control the fraction of gene expression of other clusters, by setting the maximum threshold as 100 %. we then annotate d each cell cluster according to the the expression profile of these marker genes and the expression of other canonical genes significant for different lung cell types based on the literature ( see extended results ). the annotation procedure was done iteratively. with this approach we generated two separate annotate d umaps, together with associated marker genes, for the tumour and b / h datasets. to compare cell - type abundances, we calculated the proportion of each cell type within each patient and broad cell annotation in the unenriched ( cd235 - ) samples. we contrasted cell - type proportions between groups ( tumour vs. background or luad vs. lusc ) using a wilcoxon rank - sum test. finally, we corrected for multiple testing using a two - sided bonferroni correction independently for each group analysed. the association between the relative cell - type abundance for each immune cell type was evaluated on the pearson ’ s product - moment correlation coefficients. to test consistency in cell - type annotation performed separately in tumour and b / h, we performed reference - query mapping from tumour to b / h using scarches. for the 828, 191 immune cells ( 464, 952 in tumour and 363, 239 in b / h ) identified through our separate annotations, we selected a common set of 10, 000 hvgs. we first built an scvi model and trained it on the tumour dataset using broad cell types for reference, and applied schpl method ( provided in the scarches package, parameters set to use knn classifier, 100 neighbours and with pca dimensionality reduction ) to obtain the hierarchy for the tumour cell types. we then applied the b / h dataset to the pretrained reference model for a query, and predicted b / h broad cell types based on tumour hierarchy ( probability threshold set as 0. 2 ). finally, we compared the predicted cell types with our separate annotations in b / h using a heatmap to visualise the confusion matrix.