Protocol for Chinese lung cancer evolution and microenvironment tracking under therapy study

Ethics approval and consent to participate

The CLEVER study has received approval from the Ethics Committee of Peking University People’s Hospital (approval number 2022PHB140-001). It is registered with ClinicalTrials. gov under the identifier NCT05352035. Patient recruitment for the study began in April 2022.

Patients enrollment

Sample collection

Data collection and management

According to the medical record numbers of enrolled patients, clinical, imaging, and pathological data were completely, correctly, and clearly entered into the statistical information table through medical record inquiry. Clinical data are archived and stored as a database in numbered order, with multiple backups to ensure proper storage. The original files are retained for the period specified by relevant guidelines. Original data are managed by designated personnel and stored on the hospital’s internal computer hard disk, which is not accessible from external networks. The final archived data of this study are stored only in serial numbers and do not contain any personal privacy information of patients. No patient privacy information was involved when the data were used for publication.

Nucleic acid extraction

Genomic DNA and total RNA are purified from the frozen tumor and normal tissues using the AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany). DNA in paired white cell samples is extracted using the TIANamp Genomic DNA Kit (Tiangen, Beijing, China). Cell-free DNA in the plasma samples is extracted using the QIAamp Circulating Nucleic Acid Kit (Qiagen, Hilden, Germany). Nucleic acid concentrations are quantified using either the Qubit HS dsDNA Kit or the Qubit RNA HS Assay Kit (Invitrogen, Carlsbad, USA).

Whole exome sequencing and analysis

The DNA library is prepared with an input of 200 ng genomic DNA using the VAHTS Universal Plus DNA Library Prep Kit for Illumina V2. The subsequent hybridization capture is performed using the SureSelect XT target enrichment system. The final DNA library is quantified using the Qubit 3.0 and fragment detection is conducted with the Agilent 2100. Finally, the DNA library is sequenced using paired-end 150 bp reads on the Illumina NovaSeq 6000 platform.

The sequencing data is analyzed using the WES pipeline. The sequencing data from white cells is used to exclude germline mutations. Clean reads in FASTQ format are aligned to the human UCSC reference genome (hg19/GRCh37) using the Burrows-Wheeler Aligner (BWA, v0.7.15). Sambamba (v0.6.8) is used to process PCR duplicates in the mapped BAM files. The Genome Analysis Toolkit (GATK v4.0.11.0) is employed for local realignment and base quality recalibration to compute sequencing coverage and depth. Single nucleotide variations (SNVs) and small insertions and deletions (Indels: < 50 bp) are identified using GATK MuTect2 (v1.1.4). Multiple testing correction for gene-level mutation frequency comparisons will be performed using the Benjamini-Hochberg procedure, with an adjusted p-value (FDR) cutoff of < 0.05.

Based on multi-region sequencing, region-based tumor mutation calling is combined to recall tumor-informed mutations, allowing for the acquisition of new recalled mutations for each region. Using these recalled region-level mutations, the phylogenetic tree is reconstructed using the CNIPHER method,^[24] enabling the identification of clonal and subclonal mutations. Somatic copy number alterations (SCNAs) are processed using ASCAT (V2.3), which allows for the determination of intratumoral heterogeneity in both SNVs and SCNAs, specifically assessing subclonal mutation fractions and subclonal SCNAs fractions. This comprehensive analysis facilitates a deep understanding of evolutionary events within the tumor, characterized by intratumoral heterogeneity, clonal illusion, subclonal expansion, whole genome doubling, weighted genome integrity index (wGII), and loss of heterozygosity (LOH). These insights provide a detailed depiction of the tumor’s evolutionary landscape and its underlying genetic dynamics.

RNA sequencing and analysis

To construct a transcriptome sequencing library, 1–3 μg of total RNA is prepared using the VAHTS Universal V6 RNA-Seq Library Prep Kit for Illumina (NR604-01/02). The final RNA library is quantified using the Qubit 3.0, and fragment analysis is conducted with the Agilent 2100. Additionally, the Bio-RAD CFX 96 fluorescence quantitative PCR is used to quantify the effective concentration of the RNA library (ensuring it is greater than 10 nmol/L) using the Bio-RAD KIT iQ SYBR Green. The RNA library is then sequenced using paired-end 150 bp reads on the Illumina NovaSeq 6000 platform.

The sequencing data is aligned to the human UCSC reference genome (hg19/GRCh37) using the STAR aligner (version 2.7.4a) and annotated with ANNOVAR software. Subsequent analyses include immune cell infiltration assessment, identification of differentially expressed genes (DEGs), enrichment analysis, and subtypes clustering. DEGs will be identified with DESeq2, applying the Benjamini–Hochberg FDR correction. Genes with FDR < 0.05 and |log₂FoldChange| > 1 will be considered significant. These analyses aim to explore the heterogeneous tumor microenvironment, providing insights into the complex interactions and variations within the tumor.

Single-cell RNA-seq sequencing

Fresh or frozen tissue samples are enzymatically digested to obtain single-cell suspensions, achieving a concentration of 700–1200 living cells per ml. Libraries are prepared using the Chromium Next GEM Single Cell 3’ Reagent Kits v3.1, following the manufacturer’s protocol. The purified library is sequenced on the Illumina NovaSeq 6000 platform with a sequencing length of 150 bases in a paired-end format.

The resulting FASTQ files are processed and quantified using the GRCh38 human reference genome with the Cell Ranger single-cell toolkit (v.3.1.0). Quality control steps include filtering out: (1) low-quality cells expressing < 200 genes; (2) > 25,000 or < 500 UMIs; (3) Low-quality cells with > 20% mitochondrial gene count; (4) Cells judged as twin by DoubletFinder. All samples were integrated using the Seurat or Seurat-wrappers package (v0.1.0). Raw counts were normalized with a scale factor of 10,000 per cell and subsequently log-transformed (log [count+1]). The top 2000 highly variable genes (HVGs) were identified for integration purposes. Data scaling and principal component analysis (PCA) were conducted using default settings. Significant principal components (PCs) were selected based on the elbow plot of the standard deviations of the PCs and were utilized for neighbor detection, followed by Louvain clustering ^[25] and uniform manifold approximation and projection (UMAP).^[26] The same analytical pipeline was applied for the detection of sub-clusters.

Spatial enhanced resolution omics-sequencing (stereo-seq)

The tissue was embedded in an OCT embedding agent (Sakura, 4583) and stored at -80°C. It was then sectioned into 10 μm thick slices using a cryotome (Leica, CM1950). Tissue structure was identified through H & E staining. Appropriately sized sections were placed on the Stereo-seq chip, and tissue integrity was assessed using ssDNA staining (Invitrogen, Q10212). Following the manufacturer’s protocol, chip permeation, RNA capture, reverse transcription, and cDNA library construction were completed.

Spatial transcriptome data were processed using Spaceranger (v.1.2.2). To achieve spatial localization of cell subsets on the spatial transcriptomics (ST) slides, we integrated the Stereo-seq data with single-cell RNA sequencing (scRNA-seq) data using the Seurat package. Seurat estimated the expression abundance of each cell subset at each spatial location by analyzing the similarity of expression patterns, thereby mapping the complete spatial distribution of cell types. These results are suitable for downstream analysis.

GeoMx digital spatial profiling

A whole transcriptome atlas comprising 18,000 genes and a proteome atlas consisting of 570 proteins were generated using the GeoMx Digital Spatial Profiler (NanoString) on formalin-fixed, paraffin-embedded (FFPE) tissue sections. The slides were stained with morphological markers, including PanCK (epithelial cell marker), CD3 (lymphocyte marker), CD68 (macrophage marker), and DAPI (DNA marker). Based on these markers, 178 circular and polygonal regions of interest (ROIs), along with two negative controls, were selected for analysis.

Digital counts for each gene were standardized following quality control procedures. Protein digital counts were normalized using internal spike-in controls (ERCC) and the signal-to-noise ratio (SNR). DESeq-based differential expression analyses for both transcript and protein counts will include FDR control via the Benjamini–Hochberg method, reporting only targets with adjusted p-values < 0.05 as significant.

Imaging mass cytometry

The paraffin-embedded tissues were sectioned into 4 μm slices and heated at 68°C for 1 h. The slices underwent dewaxing twice in xylene, followed by rehydration through a graded series of ethanol solutions at concentrations of 95%, 85%, and 75%. Antigen retrieval was conducted by heating the slices at 100°C for 30 min, followed by natural cooling. Subsequently, the slices were washed twice with PBS-TB for 5 min each and then blocked with SuperBlock for 30 min. After blocking, the slices were washed three additional times with PBS-TB and incubated overnight at 4°C with a metal-labeled antibody cocktail (Maxpar® X8 Kit). Post-incubation, the slices were washed and treated with 1.25 μmol/L Intercalator-Ir in PBS-TB at room temperature for 30 min, followed by washes with PBS-TB and distilled water. The prepared sections were then scanned using a Hyperion Imaging System to generate multiple images.

The IMC data processing involved four steps: overflow compensation, image denoising, contrast enhancement, and cell segmentation. The range of label expression was normalized to the 99th percentile for each channel, and batch effects were corrected using Harmony (v0.1.0). Clustering of cells was performed using Rphenograph (v0.99.1) with 100 nearest neighbors, and the cluster means were visualized as an annotated heat map.^[27] The K-means clustering method with (k = 15) was employed to identify cell neighborhoods (CNs) using a clustering window of the 20 nearest cells, based on Euclidean distance. Cell-cell interactions were analyzed using an arrangement assay facilitated by the imcRtools package (v1.0.2).^[28]

Minimal residual disease detection

Upon obtaining the region-specific exome sequencing data, we integrated the data from each region to design a personalized MRD panel. Our surveillance strategy includes up to 150 somatic mutations identified in each patient, along with 21 well-established cancer-related driver mutations. In cases where a large number of somatic mutations are detected, we prioritize screening clonal mutations over subclonal mutations. To ensure the sensitive detection of low-frequency mutations, we employ ultra-deep sequencing at a depth of 100,000×, using a minimum input of 10 ng of circulating tumor DNA (ctDNA). ctDNA-based MRD is dynamically monitored to evaluate the clinical significance of ctDNA in plasma for detecting low levels of residual disease.

Multi-omics data integration

For the CLEVER cohort, we implemented a multidimensional profiling strategy encompassing three complementary analytical axes: molecular stratification (genomic whole-exome sequencing, transcriptomic RNA-seq, and proteomic GeoMx DSP protein), resolution scaling (bulk tissue analysis, single-cell resolution, and spatial mapping), and cross-platform integration. This framework was operationalized through molecularly barcoded aliquots derived from identical tumor regions, ensuring cellular and spatial concordance across omics layers.

Machine learning model development and validation

The patients recruited in the first phase of this study will be randomly divided into training set, test set and validation set according to the ratio of 6: 2: 2. Considering the “events per variable” (EPV) rule, we apply the following steps to develop a robust recurrence prediction model. First, we will not treat all molecular features as separate predictors. Instead, we will compute summary metrics to reduce dimensionality (see section 5.1). Then Least Absolute Shrinkage and Selection Operator (LASSO) regression will be applied within the training subset to identify the most informative predictors, where the penalty parameter will be selected via grid search in nested 5-fold cross-validation. This procedure typically yields less than 9 variables, satisfying EPV≥10 given an expected 80 events (26.7% of the 300 patients).

All parameter tuning will be confined to the training set to avoid information leakage. Machine learning-based classifiers will be implemented on the training set with 1000 rounds of ten-fold cross validation. In each round, the training set will be split into 10 subsets (folds); in each fold, a model will be trained on nine of the subsets and then validated on the remaining subset to tune the hyperparameters. Then we will apply all obtained models to the test set and compare the performance measurements. The optimal method will be used to construct the final model. The tested algorithms included: adaptive boosting (AdaBoost), decision tree, logistic regression, linear support vector machine (Linear SVM), radial basis function kernel support vector machine (RBF SVM), naive Bayes, nearest neighbors, neural net, quadratic discriminant analysis (QDA), random forest, and extreme gradient boosting (XGBoost). All machine learning algorithms have a solid research foundation.^[35]

We will assess C statistics (area under the receiver-operator characteristic curve [ROC-AUC]), calibration, and clinical utility (decision-curve analysis) in the validation and test set. Finally, we determined an algorithm for developing predictive models, which has been widely popularized in the machine learning community and has excellent scalability, high running speed and state-of-the-art accuracy. It has better average performance than other algorithms in tests. Therefore, it is chosen for further optimization. The model willed be externally and respectively validated on two additional cohorts of NSCLC patients from published database^[10,11] and subsequent phases of the CLEVER study to confirm reproducibility and generalizability. If necessary, we will recruit a multi-center prospective multi-region cohort from collaborating institutions of medical alliance for further verification.

Statistical analysis

The Wilcoxon test is used to compare distribution differences in genomic and transcriptomic biomarkers between groups, while Fisher’s exact test is applied to compare clinical variables and the number of gene mutations between groups. All p-values are two-sided, with P < 0.05 considered statistically significant.