2.1 Data source

The transcriptome sequencing and genetic mutation data for LUAD patients were acquired from The Cancer Genome Atlas (TCGA) database (Workflow Type: TPM; https://portal.gdc.cancer.gov/repository), and corresponding clinicopathologic info, including gender, age, American Joint Committee on Cancer (AJCC) pathologic stage, was originated from UCSC xena online program (https://xena.ucsc.edu/). As the validation set, GSE31210 dataset including the gene expression profiling and overall survival (OS) info of 227 LUAD patients and ICGC-LUAD dataset including the gene expression profiling were retrieved from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) and International Cancer Genome Consortium (ICGC; http://dcc.icgc.org) databases respectively. Besides, the 16 pairs LUAD and corresponding adjacent normal lung tissue samples for quantitative polymerase chain reaction (PCR) and the other 12 pairs for immunohistochemistry (IHC) from the First Affiliated Hospital of Guangxi Medical University served as the GX-LUAD cohort. The previously published procedures for tissue storage, RNA extraction, and quantitative PCR for reverse transcription were followed [10]. The designed sequence info of PCR primer is listed as below: GAPDH, forward: GTCAGCCGCATCTTCTTT, reverse: CGCCCAATACGACCAAAT. HOXC6, forward: AGACGAGAAAGAGAGTGGGAGA, reverse: GCTGGAACTGAACACGACATTCT. Moreover, the expression of HOXC6 was determined employing 2 ∆∆CT methodology [11].

2.2 Exploring the differential expression and prognostics value of HOXC6 in LUAD

The HOXC6 expression levels were compared between LUAD and normal or neighboring lung tissues in the TCGA-LUAD, ICGC-LUAD and Guangxi cohorts respectively. In addition, the Human Protein Atlas (HPA) database (https://www.proteinatlas.org) was applied to reveal the HOXC6 expression at a protein level in LUAD and normal lung tissues. In the survival analysis, we applied the “survival” and “survminer” R packages to obtain the optimized cut-off value of HOXC6-expression level for prognostic value evaluation. Then, high- and low-HOXC6 expression subgroups were categorized by cut-off value of HOXC6 expression. Kaplan–Meier survival plots for OS were depicted between two subgroups in the TCGA-LUAD and GSE31210 cohorts. Employing the R packages "ggpubr" and "limma", the differences in HOXC6 expression levels across various clinicopathological characteristics were examined in the TCGA-LUAD cohort and the results were illustrated as boxplot and heatmap.

To validate the protein expression level of HOXC6 in tumor tissue, we further performed the immunohistochemistry (IHC) of the 12 pairs LUAD and corresponding adjacent normal lung tissue samples in the Guangxi cohort. The paraffin section, staining and quantitative evaluation procedures of IHC was previously described [12, 13]. As directed by the manufacturer, all primary antibodies were diluted and incubated for the entire night at 4 °C. The following primary antibodies were employed in this investigation: HOXC6 (1:200, Bioss, China).

2.3 Nomogram construction

Adopting a multivariate Cox proportional risk regression model in “survival” and “rms” R packages [14, 15], the nomogram relying on the TCGA-LUAD data was developed to project the probabilities of 1-, 3-, and 5-year OS. The LUAD patients with insufficient clinicopathological info were eliminated. This prediction model combined HOXC6 expression with clinicopathological indicators, including gender, age, AJCC pathologic Tumor Node Metastasis (TNM) stage, which are routinely available in our clinical practice. The consistency levels of 1-, 3-, and 5-year calibration curves were utilized to evaluate the discrimination efficiency of nomogram.

2.4 Gene expression and TMB correlation analysis of HOXC6

Based on the transcriptome profile of TCGA-LUAD at mRNA level, the Pearson correlation analyses were performed to discover HOXC6-related genes. Top 10 HOXC6-related genes were identified under specific situations that the cut-off value of correlation coefficient was set to 0.6 and P-value was set to 0.001. The above comprehensive results were summarized in a circus plot. Moreover, Pearson correlation analyses were further performed to disclose the relationships between HOXC6 and the genes associated with immune checkpoints, including IDO1, LAG3, CTLA4, TNFRSF9, ICOS, CD80, PDCD1LG2, TIGIT, CD70, TNFSF9, ICOSLG, KIR3DL1, CD86, PDCD1, LAIR1, TNFRSF8, TNFSF15, TNFRSF14, IDO2, CD276, CD40, TNFRSF4, TNFSF14, HHLA2, CD244, CD274, HAVCR2, CD27, BTLA, LGALS9, TMIGD2, CD28, CD48, TNFRSF25, CD40LG, ADORA2A, VTCN1, CD160, CD44, TNFSF18, TNFRSF18, BTNL2, C10orf54, CD200R1, TNFSF4, CD200, and NRP1. Tumor mutation burden (TMB) makes a difference in predicting prognosis and immunotherapy response for LUAD patient [16]. We summarized the genetic mutation data of TCGA-LUAD and the TMB of each tumor sample was worked out. Then, the spearman correlation was implemented between the TMB value and HOXC6 expression.

2.5 GSEA function enrichment analysis

Based on TCGA-LUAD dataset, the Gene Set Enrichment Analysis (GSEA) was carried out between the low- and high-HOXC6 expression subgroups to identify the possible regulating pathways of HOXC6, in which R package, including “org.Hs.eg.db” and “clusterProfiler” and the gene set data, including “c2.all.v7.0.symbols.gmt” and “c5.all.v7.0.symbols.gmt” were employed [17]. Relevant enrichment pathways were determined by functional terms that satisfied the requirements of a nominal P < 0.05 and a false discovery rate (FDR) < 0.25. Then, the top 7 pathways with the smaller P-value were visualized.

2.6 Immune cell infiltration analysis

Cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) is a novel analytical methodology to estimate the relative abundance of different immune cell types in tumor utilizing gene expression data [18]. Using CIBERSORT (R script v1.03), the proportions of various immune cells were assessed for each tumor sample based on the TCGA-LUAD gene expression profiling. Then, to uncover the correlation between HOXC6 and tumor immune infiltration, Spearman correlation analyses and relative abundance difference analysis were conducted. Additionally, employing “estimate” R package and reference data, including “commonGenes.gct” and “estimateScore.gct”, the stromal score, immune score and ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) score of each TCGA-LUAD samples were calculated [19]. And we further explored whether there was a difference in these scores between the low- and high-HOXC6 expression subgroups.

2.7 Drug sensitivity prediction analysis

In the drug sensitivity assessment study, the “pRRophetic” R package was adopted to estimate the half-maximal inhibitory concentrations (IC50s) of TCGA-LUAD cohort to chemotherapeutic and targeted medicines linked to LUAD managements by ridge regression method [20]. In addition, The Cancer Immunome Atlas (TCIA, https://tcia.at/) online program provides the immunophenotypic score (IPS) of a tumor sample using the transcriptome info of solid cancers from TCGA dataset. To some extent, the therapeutic responsiveness to immune-checkpoint inhibitors (ICIs), including anti-programmed cell death protein 1 (PD-1) and anti- cytotoxic T lymphocyte antigen-4 (CTLA-4) antibodies, could be reflected in IPS value [21]. A higher IPS suggests a possibly finer reaction to immunotherapy. Next, we further probed into the differences of IC50s and IPS in the low- and high-HOXC6 expression subgroups in the TCGA-LUAD cohort.

2.8 Statistical analysis

The data analysis and result visualization were performed using R (v4.3.1). For continuous data, the Wilcoxon rank-sum test was utilized to compare data between various subgroups, whereas for categorical data, the chi-square test was applied. The difference in HOXC6 expression levels between the LUAD and the accordingly neighboring normal tissues of the TCGA and GX-LUAD cohorts was investigated using a paired t-test. P < 0.05 was given as the threshold if not specifically mentioned.