The GSE43458, GSE83456, and GSE42834 datasets were obtained from NCBI Gene Expression Omnibus (GEO) ( https://www.ncbi.nlm.nih.gov/geo/ ). GSE43458 consists of 80 lung cancer samples and 30 control samples run on an Affymetrix Human Gene 1.0 ST Array [HuGene-1_0-st]. GSE83456 contains 45 pulmonary tuberculosis samples and 61 control samples run on an Illumina HumanHT-12 V4.0 Expression BeadChip. GSE42834 contains 20 controls, eight patients with lung cancer, and 19 pulmonary tuberculosis patients, also run on an Illumina HumanHT-12 V4.0 Expression BeadChip. The R packages affy and annotate were used to process the raw data, make an expression matrix, and match probes to their gene symbol. The R package Limma was used to screen the DEGs based on the GSE42834 data. All DEGs were screened with q-value < 0.001 and |log2FC| > 0.5 as thresholds. The common differential genes in these results were selected for functional analysis.

This study was performed according to the Declaration of Helsinki (2013) of the World Medical Association. The study protocol was approved by the Ethics Committees of The Xinjiang Uygur Autonomous Region Chest Hospital (approval number 2021BAT011). Fresh primary tissues were collected from 15 pathologically confirmed patients undergoing surgery for lung cancer, pulmonary tuberculosis, or lung cancer combined with pulmonary tuberculosis at The Eighth Affiliated Hospital of XinJiang Medical University (Urumqi, China). Histological diagnosis of tumors was performed and confirmed by two pathologists. Samples were immediately frozen in liquid nitrogen and stored at −80°C until further analysis. The clinicopathological features of the 15 patients are presented in Supplementary Table S5.

Based on the expression and clinical pathological data of the GSE43458 and GSE83456 datasets, the genes showing the top 60% variance were selected for weighted gene co-expression network analysis (WGCNA). This was used to calculate the correlation coefficients between genes for clustering and to construct a co-expression weighted network. The hierarchical clustering function was used to classify genes with similar expression profiles into modules based on the topological overlap matrix (TOM) dissimilarity with a minimum size of 30 for the gene dendrogram. The blue modules were significantly associated with TB, and the turquoise and yellow modules were significantly associated with LC. The dissimilarity of module eigengenes (MEs) was calculated to choose a cut-off to merge some modules. Lastly, the network of eigengenes was also visualized.

To determine the appropriate modules, we computed the association between MEs and clinical traits. The log10 transformation of the P-value (GS = lgP) in the linear regression between gene expression and clinical data was then used to establish gene significance (GS). The average GS for every gene in a module was also defined as the module membership (MM). The module associated with a clinical feature was the one with the highest absolute MM ranking out of all the modules that were chosen.

The protein–protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org), which covered almost all functional interactions between the expressed proteins, and interactions with a combined score greater than 0.4 were considered statistically significant. The results of this investigation were shown using Cytoscape software (version 3.7.0).

The R packages clusterProfiler and org.Hs.eg.db were used for Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the hub genes for all data.

Genomic DNA (gDNA) was extracted from the recently frozen tissues using the QIAmp DNA Tissue Kit (TIANGEN, Beijing, China) following the manufacturer’s instructions. The quantity and purity of the gDNA were assessed using a Qubit® 2.0 fluorometer (Thermon Fisher Scientific. Waltham, MA, USA) and a NanoDrop 2000 (Thermo Fisher Scientific, Inc.). The fragmentation status was evaluated using the Agilent 2200 TapeStation system using the Genomic DNA ScreenTape assay (Agilent Technologies, Santa Clara, CA, USA) to produce a DNA integrity number. An additional quality control step to determine the DNA integrity was performed using a multiplex PCR approach.

A total of 300 ng of each gDNA sample, based on the Qubit quantification, was mechanically fragmented (duty factor, 10%; peak incident power, 175 W; cycles per burst, 200; treatment time, 240 s; bath temperature, 2–8°C) on an M220 focused ultrasonicator (Covaris, Inc.). The target DNA fragment size was 350 bp. Subsequently, 200 ng of sheared gDNA was used to perform end repair, A‑tailing, and adapter ligation with a library preparation kit (Agilent Technologies, Inc.) according to the manufacturer’s protocol. Subsequently, the libraries were captured using Agilent SureSelect XT custom 0.5–2.9M probes (Agilent Technologies, Inc.) and amplified. The captured libraries were sequenced on an Illumina NovaSeq 6000 PE150 (Illumina Inc.) for 2 × 150 paired end reads, using Cycle Sequencing v3 Reagents (Illumina).

Clean data were obtained after filtering out the adapters and reads with a proportion of N > 10%, with N being the unidentified bases in the sequencing process, using fastp (fastp 0.20.0). Low-quality bases (Phred score < 15) were excised from the 3′ ends of reads. Only sequenced samples with >80% of data with a quality score ≥ Q30 (95% of base call accuracy) were used in the analysis. Reads with length < 50 bp after excision were removed. The clean data were mapped to the human reference genome (University of California Santa Cruz ID: hg38) using the Burrows–Wheeler Alignment algorithm (BWA 0.7.17). The alignment in SAM format was converted to BAM files using SAMtools (SAMtools 1.9.0). Next, the genome analysis toolkit (GATK; v4.0.2.1) was used for sorting, duplicate marking, and base recalibration. The final BAM files were analyzed using QualiMap v.2.2.1 to provide a global overview of the data, including mapped reads, mean mapping quality, and mean coverage. The variants (single nucleotide variants (SNV) and insertion–deletion mutations (InDels)) were called for unpaired tumor sequences with 40 pooled blood samples (from healthy individuals) using the GATK mutect2 tumor-only mode (parameter: af-alleles-not-in-resource, 0.00025%), and germline mutations and contaminations were filtered out using GATK FilterMutectCalls (parameter: max-germline-posterior, 0.995). Somatic variants were annotated using the ANNOVAR software tool.

The following filter conditions were used to identify candidate somatic alterations: i) all variations with COSMIC evidence (http://cancer.sanger.ac.uk/cosmic) were retained; and ii) SNV acquisition conditions: (1) tumor samples require at least 10× coverage; (2) at least 5× coverage supports mutant alleles in tumor DNA; (3) mutant allele frequency (AF) ≥ 0.05; and (4) genomic locations with mutant allele frequencies greater than 0.1% in the Thousand Genomes Project and Exome Aggregation Consortium (ExAC) were filtered out (AF ≥ 0.001).

The mutational landscape across the cohort was created using the maftools package in R software (R 4.3.0, R Core Team; https://www.R-Project.org). A cut-off value of an adjusted P-value (p.adjust) < 0.05 was used to identify significantly enriched terms.

WGCNA analysis was performed for both GSE43458 and GSE83456. To construct a scale-free topology network, soft threshold powers (β) of 12 in GSE43458 (scale-free R2 = 0.80) ( Figure 1A ) and 12 in GSE83456 (scale-free R2 = 0.80) ( Figure 1B ) were estimated. For GSE43458, the hierarchical clustering tree revealed that 17 co-expression modules were clustered ( Figure 1C ), and the orange module was negatively correlated with the LC proportion (Cor = 0.9, P = 1×10−22) ( Figure 1E ). For GSE83456, dynamic hybrid cutting clustered 20 co-expression models ( Figure 1D ), with the black module having the strongest positive correlation with the TB proportion (Cor = 0.83, P = 1×10−14), and the salmon module showing the strongest negative correlation (Cor = 0.85, P = 6×10−16) ( Figure 1F ). In the orange module, scatter plots showed strong negative correlations between MM and GS for LC (Cor = 0.92, P = 1×10−200) ( Figure 1G ); strong positive correlations were also observed between MM and GS for TB in the black module (Cor = 0.89, P = 1×10−200), as well as large negative correlations between MM and GS in the salmon module (Cor = 0.89, P = 1×10−200) ( Figure 1H ). Hence, these three modules were selected for in-depth investigation. A total of 5157 and 1042 genes were incorporated in the black and salmon modules, respectively, while the orange module contained 2743 genes.

To identify differential genes in LC, TB, and TB patients with LC, we analyzed the DEGs in three modules using another dataset, GSE42834 from the GEO database. We set the cut-off as |log2FC| > 0.5 and q-value < 0.001 to screen DEGs from GSE42834. Figure 2A shows a volcano plot of the DEGs. We overlapped the DEGs and the genes in the three modules (LC vs control, LC vs TB, and TB vs control) by Venn diagram and found that 3606 common genes were present in all three modules ( Figure 2B ). Figures 2C, D demonstrate the GO and KEGG analyses of these 3606 genes. Extracellular matrix was the most enriched cellular component (CC) term, G protein-coupled receptor activity was the most enriched molecular function (MF) term, and neuroactive ligand–receptor interaction was the most enriched KEGG pathway.

To assess the function of the hub genes, we extracted the common genes derived from WGCNA and DEGs. As shown in Figure 3A , 278 common genes at the intersection of the three hub gene sets were detected and visualized via a Venn diagram. Subsequently, we performed GO and KEGG analyses on these 278 genes. Cell–substrate adhesion was the major enriched biological process (BP) term, and collagen-containing extracellular matrix and growth factor activity were the major enriched CC and MF terms, respectively ( Figure 3B ). The RAS and PI3K–AKT signaling pathways were the main enriched KEGG pathways ( Figure 3C ). The PPI network was constructed with 278 genes. The EGFR pathway is an oncogenic pathway in human non-small cell lung cancer (NSCLC), which also affects the levels of some pathway-related binding proteins or downstream activities (14, 15). Pathway analyses further revealed that the levels of some EGFR-related genes were altered, suggesting that EGFR might serve as a regulatory signal node in the development of TB-associated lung cancer ( Figure 3D ).

To further investigate the role of EGFR in lung cancer and pulmonary TB, we performed exon sequencing on samples from lung cancer, lung cancer associated with pulmonary tuberculosis, and pulmonary tuberculosis patients. Supplementary Figure S1 presents a summary of the MAF files ( Supplementary Table S6 ) for the 15 patients. Totally, 2094 meaningful variations in 497 genes were identified. Totally, 2094 meaningful variations in 497 genes were identified. A waterfall plot depicts 25 of the genes containing indel mutations ( Figure 4 ). For the SNVs, T > C was the most frequent SNV class. The median number of variants identified in the 15 samples was 5 (range, 1–24).

Figure 4 presents the mutational profile of the 15 patients with LC and pulmonary TB, including 25 mutated genes, organized by the TB, LC, and TB+LC groups. The mutated gene with the highest frequency was ADCK5 (67%). EGFR, one of the most frequently mutated genes in lung cancer, had a mutation rate of 20%, similar to the previously reported EGFR mutation rate of 5%–30% in LC (16). The other recurrently altered genes were mucin-3A (MUC3A; 60%), keratin-associated protein 5–7 (KRTAP5–7; 53%), killer cell immunoglobulin-like receptor 2DL4 (KIR2DL4; 33%), zonadhesin (ZAN; 27%), fatty acid desaturase 6 (FADS6; 27%), cytochrome P450 family 1 subfamily B member 1 (CYP1B1; 27%), and ataxin-3 (ATXN3; 20%). Supplementary Figure S2 illustrates the frequency of mutations and the resulting protein structure.

To further investigate the biological functions of the mutated genes, KEGG signaling pathway enrichment analyses were performed. Supplementary Table S7 shows enriched signaling pathways associated with the mutated genes. The results of the KEGG signaling pathway analysis are presented in Figure 5 , showing that the genes are involved in the PI3K–AKT signaling pathway and non-small cell lung cancer, which is consistent with our previous analysis.

EGFR had the highest degree in the aforementioned PPI, with a mutation rate of 20% in lung cancer. This implies that EGFR is involved in the progression of lung cancer. Thus, we investigated the relationship between EGFR and other genes. The Venn diagram in Figure 6A depicts the genes common to both the 278 hub genes and all mutated genes, including EGFR, HSPA2, CECR2, and LAMA3 ( Figure 6A ). To investigate the effects of EGFR expression on HSPA2, CECR2, and LAMA3, we performed the CIBERSORT algorithm on 15 tumor samples to calculate the expression of EGFR and the three key genes in each sample. As shown in Figures 6B–D , CECR2, LAMA3, and HSPA2 were positively correlated with EGFR expression.