Identification of molecular biomarkers associated with non-small-cell lung carcinoma (NSCLC) using whole-exome sequencing

2.1Sample collection and diagnosis

A total of 19 NSCLC cases (EGFR, ALK and ROS1 negative) with available clinical follow-up were retrieved from the Department of Pathology, AIIMS, New Delhi. The haematoxylin and eosin-stained slides were analysed and histological type of the tumour was determined according to World Health Organization (WHO) 2021 classification of thoracic tumours. Blocks showing more than 80% tumour component in their respective sections were used for WES. Treatment and follow up details were retrieved from case record files from the Department of Medical Oncology, AIIMS, New Delhi.

2.2Formalin fixed paraffin embedded (FFPE) DNA isolation and repair

DNA extraction was performed using FFPE DNA tissue extraction kit (A2352, Promega, USA) according to the manufacturer’s instructions. FFPE DNA was repaired and purified using Gene JET FFPE DNA Purification Kit (K0881, Thermo Scientific, USA) according to the manufacturer’s instructions. Quantity and purity of gDNA were assessed by Qubit® 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and NanoDrop ND-1000 (Thermo Scientific, USA).

2.3Sample sequencing

Sequencing libraries were prepared using the SureSelect All Human Exon V5 Kit (USA) according to the manufacturer’s instructions. The final enriched pooled were sequenced on Illumina HiSeq 2000 platform (Illumina Inc., USA) generating 2 × 150 bp paired-end reads. Image analysis and base calling were carried out by Illumina software (CASAVA) with default parameters. Demultiplexing and FASTQ file generation from Illumina basecall (BCL) files was performed using Bcl2fastq conversion software.

2.4Variant calling and quality control

Quality of raw reads (FASTQ files) was examined using FastQC [5]. Adaptor and low-quality sequences were trimmed using Trimmomatic software [6]. Paired clean reads with longer than 50 bases were aligned against the Genome Reference Consortium Human Build 38 patch release 7 (GRCh38.p7) using BWA-MEM algorithm of Burrows-Wheeler Aligner (BWA 0.6.1) [7]. SAM/BAM post-processing steps including SAM to BAM conversion, sorting, adding read group information, mark duplicates, and base quality score recalibration were performed using the Genome Analysis Toolkit (GATK 4.0.6.0) [8, 9]. The quality of the recalibrated BAM files was checked with QualiMap v2.0.2 [10]. Finally, a genomic variation, including single-nucleotide polymorphisms (SNPs) and small INDELs (insertion and deletion) were detected for each sample individually using GATK HaplotypeCaller in GVCF mode (-ERC GVCF), and the results were combined using GenotypeGVCFs. Raw variant calls were soft filtered using GATK VariantFiltration based on the following parameters: LowCoverage (DP < 5), LowQual (Q < 50), StrandBias (FS P-value> 60), SNV cluster (three or more SNVs within 10 bp), Poor Mapping Quality (> 10% of reads have nonunique alignments).

2.5Variant annotations

The impacts of variants on protein structure or function were predicted using SnpEff [11] and ANNOVAR [12]. It compiles prediction scores from multiple algorithms including PhyloP, SIFT, LRT, SiPhy, Polyphen-2, GERP++, MutationAssessor, Fathmm, MutationTaster, CADD, and MetaSVM. In addition to these tools, variants were re-annotated using the germline/population databases (dbSNP, 1000 Genomes and ExAC) [13, 14], and cancer/somatic databases (COSMIC, TCGA and ICGC) [15, 16]. ClinVar [17] and My Cancer Genome (http://www.mycancergenome.org) were used to determine the clinical significance of each variant, while drug databases (PharmGKB, OncoKB) [18, 19] were utilized to gain information about the treatment implications of specific cancer gene mutations and how these mutations affect treatment response.

2.6Additional filters to reduce false positives somatic variants

A high-confidence somatic variant for tumor samples without a matched normal control was selected based on the following criteria: 1) mutations were considered true positives if they have a) QUAL ⩾ 20, b) genotype quality (GQ) ⩾ 20, c) mapping quality (MQ) ⩾ 20, d) coverage depth at candidate site (DP) ⩾ 20, e) QualByDepth (QD) ⩾ 2.0, and g) frequency ⩾ 25% in tumor samples [20], 2) all common variants with minor allele frequency (MAF) of > 1% in the germline/population databases (ExAC, and 1000 Genomes) were filtered out since those variants are deemed polymorphic/benign rather than pathogenic somatic driver mutations, 3) known germline variants reported at dbSNP (version 151) were excluded, and alterations listed as known somatic variations in COSMIC, The International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA) were retain [21], 4) MAF threshold of 0.0001 was used in the gnomAD or TopMed database to filter variants for the somatic mutation [22], 5) variants were considered if the variant allele frequency (VAF; also known as variant allele fraction) is deviating from germline polymorphisms (∼0.5/1 for heterozygous/homozygous) [23, 24], 6) variants were considered if they are truncating variants (nonsense mutations, frameshift deletions/insertions, mutations located at exon-flanking regions, and highly conserved intronic splice sites), or apparent missense mutations predicted to be pathogenic by in-silico prediction tools, 7) synonymous variants that were not previously reported as pathogenic and not predicted to alter splicing were filtered out.

2.7Enrichment analysis and candidate gene prioritization

Known or predicted variants to be involved in the lung or in related cancers were predicted using DisGenNET [25]. The Human Gene Damage Index server (http://lab.rockefeller.edu/casanova/GDI) was used to predict LoF-intolerant genes. The gene product physically interacts with a protein encoded by a known disease gene was explored using NetworkAnalyst [26]. The DAVID tool was used to perform KEGG pathway and GO functional-enrichment analyses of DEGs. Gene product in a pathway associated with the disease and gene expressed in the tissue or organ of interest was retrieved from the literature.

2.8Survival and expression analysis

A web-based tool, GEPIA (Gene Expression Profiling Interactive Analysis; http://gepia.cancer-pku.cn), was used in the study of lung cancer patients to assess gene expression between tumor and normal data from the Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx). During the expression analysis, the threshold included an expression fold change of 1.5 between cancer and normal tissues, as well as a p value of 0.05. This relationship was then visualized with a boxplot. A survival analysis as well as a correlation analysis between two genes were also conducted using GEPIA.

2.9Sanger sequencing

PCR was carried out on ABI Palm thermal cycler (Applied Biosystem, California), using both forward and reverse primers for greater accuracy and the results were analysed using SeqMan II software (DNASTAR 5.07). The mutations with both base counts more than 10% and QV (quality value) more than 20 were considered to be trusted mutations.

3.1Patient characteristics and sequencing statistics

A case study of 19 patients was included in this study, where 10 were ADCA (EGFR, ALK and ROS1 wild type), and nine were SQCC histological subtypes of NSCLC. Patients were early onset with the average diagnostic age of 56 years, where male: female ratio was 5.6: 1. Among, three patients were the non-smoker, whereas one case was with unknown smoking history. Nearly 75% of patients were present with co-morbidities, where all patients were either in stage III or IV. A total of 13.48 GB raw data and 11.98 GB processed data were generated per exome for the tumor sample of ADCA, whereas on an average 13.76 GB raw data and 12.24 GB processed data were obtained for the tumor sample of SQCC using WES. A higher percentage of reads were aligned to the human reference genome (GRCh38) in the tumors of ADCA patients (98.87%; range 96.48–99.95%) compared to tumors of SQCC patients (97.40%; range 85.03–99.51%), indicating that the generated dataset was highly relevant with the reference genome. The average GC content (49.74%) in the tumors of ADCA patients ranged from 42.10 to 64.63%, whereas average GC content (48.50%) in the tumors of SQCC patients ranged from 41.81 to 52.28%. Clinical characteristics and sequencing summary of lung cancer patient’s participants in this study are listed in Table S1 in Supplemental File 1.

3.2Detection and characterization of SNVs

After initial variant filtering (as described in methods), a total of 1,157,921 (single-allelic 1,157,792 and multi-allelic 129) variants were retained in the tumors of ADCA patients (n= 10) which was slightly higher than total variants (1,076,209; single-allelic 1,076,069 and multi-allelic 140) detected in SQCC patients (n= 9). The number of variants per chromosome ranged from 6,227 (chrY) to 98,788 (chr4) in the tumors of ADCA subtype, whereas ranged from 8,385 (chrY) to 104,645 (chr4) in the tumors of SQCC subtype. The variants rate per chromosomes varied from 1,925 (chr4) to 9,190 (chrY) and revealed on average 1 variants after every 2,667 bases in ADCA subtype, whereas it was after every 2,869 bases in SQCC subtype. We observed higher known variants i.e., 616,999 (53.285%) in ADCA subtype as compared to 537,980 (49.988%) in SQCC subtype. The distribution of variants by their type disclosed 1,066,489 SNPs, 38,751 insertions and 52,681 deletions in ADCA subtype. However, the different distributions of insertions/deletions (35,799/49,010) and SNP (991,400) was observed in SQCC subtype. The identified SNPs from NSCLC were categorized into two clusters based on nucleotide substitutions i.e., transitions (A/G and C/T) and transversions (A/C, A/T, C/G, and G/T). The transition-to-transversion (Ts/Tv) ratio was slightly higher (1.698 for all SNPs and 2.214 for known SNPs) in ADCA subtype compared to Ts/Tv ratio of 1.4859 (all SNPs) and 2.137 (known SNPs) in SQCC subtype. Among detected SNPs, 23.05% were heterozygous, and 76.93% were homozygous in ADCA subtype, whereas it was 21.35% and 78.64% in SQCC subtype, respectively. The ratio of heterozygous SNVs to homozygous SNVs (Het/Hom ratio) was 0.29 and 0.27 in ADCA and ADCA, respectively where the lower value was associated with true positive variants. The ratio of nonsense to missense mutations (0.007), and missense to silent mutations (0.829) in ADCA subtype was nearly similar to 0.007 and 0.863, respectively in SQCC subtype. The ratio of nonsense to missense and missense to silent mutations in the human genome may reflect a role for natural selection, especially purifying selection. The ‘GAT’ codons have been replaced maximum times by ‘GAC’ codons in both ADCA (588) and SQCC (665) subtype. The characterization of SNVs in ADCA and SQCC subtype, revealed that ADCA was more genetically unstable compared to SQCC. Variant’s summary identified by whole exome sequencing is listed in Table S2 in Supplemental File 1.

Figure 1.

Nonsynonymous somatic SNVs and INDELs identified in lung cancer patients by whole-exome sequencing. The outer-coloured ring and number indicate chromosome number and partition; the middle green ring and letters represent genes with non-synonymous SNVs and their corresponding chromosomes; the inner violet ring shows non-synonymous INDELs and their corresponding chromosomes. The star symbol signifies that gene associated with more than one mutation.

Figure 2.

Oncoplot for the somatic variants in non-small cell lung cancer (NSCLC). The graph depicting top 18 mutated genes ordered by decreasing frequency. The right barplot shows overall frequency in ADCA and SQCC-subtype. The colour box represents the type of mutations including SNV/nonsynonymous (violate), SNV/ stop gain (light blue), deletion/nonframeshift (dark blue), insertion/frameshift (light green) and insertion nonframeshift (red). The top stacked barplot shows a number of somatic mutations per sample.

3.3Detection of somatic variants in tumor only samples

To detect somatic SNV, the present study focused on missense variants in the exonic regions or splice sites. The downstream filtering by genomic location revealed a total of 6,712 and 8,000 exonic variants in ADCA and SQCC subtype, respectively. Among exonic SNVs (ADCA subtype), 2,113 were synonymous, 1,985 were nonsynonymous, 28 were frame-shift indels, 51 were nonframe-shift indels, 15 were stop-gain, 2 were stop-loss, and 2,518 were non-coding SNVs. In SQCC-subtype, 3,249 were synonymous, 2,656 were nonsynonymous, 52 were frame-shift indels, 59 were nonframe-shift indels, 36 splicing-variant, 1 was gene fusion, and 1,947 were non-coding SNVs. Exonic missense, nonsense, stop-loss, frameshift and splice site variants all have potential to affect protein function. Therefore, we excluded the synonymous variants that have no functional impact and retained 4,599 and 4,751 variants from ADCA and SQCC subtype, respectively. As rarity [27] is the key criterion to have a functional effect on the encoded protein, the filtered variants were used to eliminate the common germline mutations (minor allele frequency below 5% in population/germline databases) and as an outcome 1,642 and 2,141 variants were retained in ADCA and SQCC subtype, respectively. After excluding false positive mutations (based on additional filter criteria’s 1a-given in methodology), 500 and 734 variants in ADCA and SQCC subtype, respectively was observed. To exclude deemed polymorphic/benign variants, high-quality rare variants (MAF ⩽ 1% and QUAL ⩾ 500) were excavated which revealed 94 variants in ADCA subtype, whereas 87 variants in SQCC subtype. To identify candidates likely to have deleterious effects, combination of multiple variant annotation tools was applied that revealed the impact of amino acid changes on protein function based on the combine scores. The variants were classified as damaging (predict pathogenic by maximum number of tools), probably damaging (predict pathogenic or benign by an equal number of tools), benign (predict benign by maximum number of tools) and uncertain significance (unknown) as per the variant assessment guidelines by the American College of Medical Genetics. The alterations listed in COSMIC, ICGC and TCGA were considered as known somatic variations in this study. The final outcome revealed a total of 24 somatic variants (ADCA = 14, SQCC = 10) associated with 18 genes and were classified as known somatic variant (n= 10), deleterious variant (n= 8), and variant of uncertain significance (VUS) (n= 6) (Table 1, Fig. 1). The gene (n= 11) namely CTBP2, ESRRA, FDFT1, FOLR3, GPRIN2, HRNR, KCNJ18, KRT18, LILRA2, MTRNR2L8, and TEKT4 were observed to be mutated in both ADCA and SQCC histologic subtype of lung cancer. In addition, mutated gene (n= 4) namely LRP2, MPRIP, NYX, and TBP were observed in histologic subtype of ADCA only, whereas FBXO6, MIR3689F, and UMPS mutated genes were found in SQCC subtype only (Fig. 2). Out of 24 somatic variants, 19 (79.17%) variants had a previously known dbSNP ID while the remaining 5 (20.83%) was unassigned, new variants. The 17 variants had higher mutations frequency (⩾ 25%) in tumor samples, whereas 21 variants were observed to be true somatic mutation based on variant allele frequency which is used to infer whether a variant comes from somatic cells or inherited from parents when a matched normal sample is not provided.

Figure 3.

Hub genes in significant network modules. The hub genes with high degree and high betweenness were denoted with red colour.

Figure 4.

Scatter plot illustrating enriched KEGG pathways and gene ontology. Top 40 KEGG pathways are depicted in the Fig. 3A. The rich factor was determined by dividing the number of genes enriched in a pathway by the total number of genes annotated in that pathway. Figure 3B, C, and D show the top 10 Biological Processes (BP), Molecular Functions (MF), and Cellular Components (CC), respectively. The colour and size of the dots denote the range of the -log P-value and the number of genes in the shown pathways, respectively. The scatter plot was made using R software v4.0.3.

3.4Knowledge-driven variant prioritization

Figure 5.

Prognostic roles of potential key genes in the NSCLC patients. Survival curves are plotted for NSCLC cancer patients. (A) Overall survival: CTBP2, GPRIN2, KRT18, LILRA2, LRP2, UPMS; (B) Disease free survival: CTBP2, KRT18, LRP2.

Figure 6.

Analysis of potential key genes expression level in NSCLC patients. The red and gray boxes represent cancer and normal tissues, respectively. (A) LRP2; (B) GPRIN2; (C) FOLR3; (D) MPRIP; (E) KRT18 and (F) LILRA2; LUAD: Lung adenocarcinoma; LUSC: Lung squamous cell carcinomas.

Though, we followed the guidelines suggested for experimental design and variant filtering, yet we obtained more candidates with likely functional effects than can be verified experimentally. In this study, high priority candidates were selected based on the biological hypothesis. The analysis of known or predicted variants to be involved in the lung or in a related cancer revealed nine genes (CTBP2, ESRRA, FBXO6, FDFT1, KRT18, MPRIP, TBP, LILRA2 and UMPS) associated with the lung carcinoma. In addition, four genes i.e., LRP2 (bone, breast, colorectal, pancreatic, prostate, and renal cell carcinoma), HRNR (breast, liver and pancreatic carcinoma), FOLR3 (breast and ovarian carcinoma) and TEKT4 (breast and thyroid carcinoma) were involved in other’s cancer types. The genes namely MIR3689F, MTRNR2L8, NYX, GPRIN2 and KCNJ18 were observed to be not associated with any type of cancer and considered as low priority genes for the validation. The encoded proteins of 16 genes, have loss-of-function (LoF) variants that damage or eliminate them. The LoF-intolerant genes (CTBP2, GPRIN2, HRNR and TEKT4) were classified as extremely loss of function intolerant (pLI ⩾ 0.9), whereas gene (MTRNR2L8) with low pLI scores (⩽ 0.1) was considered as LoF-tolerant (common loss-of-function variants) and was not selected for the validation. We also prioritized genes based on the interactome of known disease-associated proteins. The analysis disclosed 13 seeds (out of 18 genes) associated with 584 nodes and 622 edges in the network. The genes, FBXO6 (degree 153; betweenness centrality 77253.04), TBP (degree 152; betweenness centrality 78171.64), KRT18 (degree 89; betweenness centrality 44603.52), CTBP2 (degree 64; betweenness centrality 33732.54), ESRRA (degree 44; betweenness centrality 21895.15), LRP2 (degree 38; betweenness centrality 20031.73), MPRIP (degree 24; betweenness centrality 10184.34), FDFT1 (degree 17; betweenness centrality 7452.33), UMPS (degree 15; betweenness centrality 7266.93), and HRNR (degree 14; betweenness centrality 6258.26) were observed to be the most highly ranked hub genes in this study (Fig. 3). Moreover, relevant information for the genes of interest was retrieved from the literature. Ten gene/proteins including HRNR, KCNJ18, ESRRA, MPRIP, FBXO6, FOLR3, FDFT1, UMPS, KRT18, and LILRA were observed to be overexpressed in lung cancer which might have a potential role in cancer development, proliferation, and metastasis. Remarkably, most of the genes were involved in important cancer-related pathways including pathways in cancer, small cell lung cancer, and non-small cell lung cancer. The PI3K-Akt signaling pathway, ECM-receptor interactions, cell adhesion, focal adhesion and the cell cycle may also play important roles in lung cancer pathogenesis (Fig. 4). Outcomes of GO enrichment analysis showed that 1) for biological processes (BP), genes were significantly enriched in DNA dependent transcription, transcription from RNA polymerase II promoter, initiation and transcription initiation from RNA polymerase II promoter; 2) for cell components (CC), genes were significantly enriched in nucleoplasm, organelle lumen and membrane enclosed lumen; 3) for molecular function (MF), genes were enriched in transcription from RNA polymerase II promoter, transcription factor binding and positive regulation of transcription, DNA dependent (P< 0.05, Fig. 4). Close observation showed that the variant in the gene MIR3689F (miRNA) and MTRNR2L8 (it is unclear if this is a transcribed protein-coding gene, or if it is a nuclear pseudogene of the mitochondrial MT-RNR2 gene) was incorrect for this study and hence not selected for the validation. The characteristics of candidate variants from public resources and published literature are given in Table 2.

Figure 7.

Correlation analysis of potential key genes in NSLCL. TCGA and/or GTEx expression data show significant correlations between LRP2, FOLR3, LILRA2, and GPRIN2.

3.5Expression analysis of potential biomarkers for NSCLC

CTBP2, GPRIN2, KRT18, LILRA2, LRP2, and UMPS mutations are associated with lower overall survival (Fig. 5a), and CTBP2, KRT18, and LRP2 mutations are associated with lower disease-free survival (Fig. 5b). It is therefore possible to identify these genes as potential biomarkers for NSCLC. Next, we applied GEPIA to check gene expression levels in NSCLC tissues and in normal tissues. The expression levels of the 5 genes (LRP2, GPRIN2, FOLR3, MPRIP, LILRA2) decreased significantly, but the expression level of one gene (KRT18) increased in comparison with normal tissues. It was observed that LRP2 and MPRIP gene expression was significantly decreased in the LUSD, whereas KRT18 gene expression increased in the LUAD (Fig. 6). These genes appear to be promising therapeutic targets. Additionally, we examined the correlation between mRNA expression of LRP2 and prognosis in patients with MPRIP (P-value = 1.1e-06) and UMPS (P-value = 9.1e-10). Furthermore, there is a significant positive correlation between GPRIN2-LILRA2 (P-value = 3.4e-10), GPRIN2-LRP2 (P-value = 0), GPRIN2-MPRIP (P-value = 0), FOLR3-MPRIP (P-value = 1.6e-09) and LILRA2-LRP2 (P-value = 0) and this may play an influential role in lung cancer prognosis (Fig. 7).

3.6Somatic variants validation by Sanger sequencing

To eliminate false-positive rates of the identified somatic mutations from WES data, we selected 10 genes for Sanger sequencing validations. Interestingly we observed that seven genes were mutated in more than 60% samples and three genes were mutated in either one or two samples. Further, Sanger sequencing results showed 100% concordance in seven genes and the remaining three genes concordances were found only in 80% cases. The mutations observed along with WES and Sanger sequencing data have been depicted in the Table 3 and Fig. S1 in Supplemental File 1.

The gene-wise results of Sanger sequencing are given below:

4.1Mutated genes present in all samples of ADCA and SQCC subtype

Of interest, the three genes, KCNJ18, TEKT4, and GRIPN2 are mutated in all NSCLC samples and can serve as common diagnostic markers for both subtypes.

4.2Mutated genes present in 80% samples of ADCA and SQCC subtype

In addition to the three aforementioned genes, HRNR, FOLR3, CTBP2 and ESSRA are significantly mutated in more than 80% of the NSCLC samples. All these mutated genes are directly or indirectly play a role in tumorigenesis and can additionally serve as common pathogenetic link for subtypes of NSCLC.

4.3Histological subtypes specific mutated genes

In the current study 3 histology specific genes were emerged from WES analysis. TBP and MPRIP genes were solely associated with the ADCA subtype whereas FBOX6 was found in one case of SQCC.