Using a stem cell-like feature classifier to predict outcomes and immunotherapy response in liver cancer with computer analysis

The gene expression data of patients with HCC were obtained from “the Cancer Genome Atlas Liver Hepatocellular Carcinoma” (TCGA-LIHC; https://portal.gdc.cancer.gov/projects/TCGA-LIHC↗) and the GEO (https://www.ncbi.nlm.nih.gov/geo/↗) databases. Patients without complete information were excluded. Finally, we selected 365 and 221 patients from TCGA (9) and GSE14520↗, respectively, and GSE14520↗ was used for external validation. Since data were obtained from two databases, to ensure the data can be compared and reduce batch effects, transcript per million values were log2 transformed and normalized using the R package. We determined the correlation between gene signatures and the response of patients from the ERP117672 and GSE202069↗ cohorts to immunotherapy.

We constructed the mRNAsi model using the “one-class logistic regression (OCLR)” algorithm (5). Next, we used “Spearman correlation analysis” to determine the similarity between cancer cells of patients with HCC and stem cells. Finally, we determined the mRNAsi by calculating the correlation coefficient value between 0-1. A stemness index closer to 1 indicates a higher stemness capacity.

The optimal cut-off value for the calculated mRNAsi score was determined using the “survminer” package. Next, patients with HCC from the TCGA-LIHC cohort were categorized into high- and low-mRNAsi subgroups based on this cut-off value. We used the “linear models for microarray data (limma)” package to screen for mRNAsi-related DEGs in tumor and normal tissues in both databases (10) based on the following criteria: “|R| > 0.3” and “P < 0.05.” Next, the “false discovery rate” (FDR) method was used for correcting the results. Finally, we identified 19 mRNAsi-related DEGs for predicting prognosis using the “univariate Cox regression analysis.”

We performed “unsupervised Consensus clustering” using the “ConsensusClusterPlus” R package to determine a stemness-based classification (11). First, we employed the “k-means algorithm” for subsampling the items at 80% and categorizing them into various groups. The process was repeated 500 times. Next, we used the “consensus matrix” and “cumulative distribution function” plots for determining the optimal cluster numbers. The overall survival (OS) of patients in the stemness subtypes was determined using Kaplan-Meier (KM) survival curve analysis. Additionally, we employed the “limma” package to identify DEGs with “|log2 fold change (FC) |>1” and “P < 0.05” in the three subgroups. Finally, we performed “functional enrichment analysis” using the “WebGestaltR” package on these DEGs.

To determine the relative abundance of immune cell types in a mixed cell population in patients with HCC, we imported and analyzed the unnormalized RNA-Seq data using the “Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT)” tool (12). Next, we predicted the infiltration of stromal/immune cells in patients with HCC using the “Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE)” algorithm (13). The “ESTIMATE” algorithm calculates the stromal, immune, and ESTIMATE scores, which show stromal abundance, immune cell infiltration, and tumor purity, respectively. We performed “single sample Gene Set Enrichment Analysis (ssGSEA)” using the “Gene set variation analysis (GSVA)” R package and assessed 29 immune gene signatures to quantify the enrichment scores of immune-related terms. Finally, we analyzed the expression of several immune checkpoint genes to determine the response of patients in the two stemness subtypes to immunotherapy.

We performed an unsupervised “GSVA” R package to determine the differential activity of pathways in patients in the two stemness subtypes. The “h.all.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” gene sets were retrieved from the “Molecular Signatures database” (14) and were selected as the background gene set. In addition, the “limma” R package was used to perform differential analysis of the “Kyoto Encyclopedia of Genes and Genomes (KEGG)” and “HALLMARK” pathways in patients in the two stemness subtypes. The significantly enriched pathways were identified based on the following criteria: “FDR < 0.05” and “log2 FC > 1.”

We used “univariate Cox regression analysis” to analyze 816 DEGs in the three subgroups. Genes with “P < 0.001” and “hazard ratio (HR)> 1” were classified as risk genes, and genes with “P < 0.001” and “HR < 1” were classified as protective genes. To accurately predict the status of the stemness subtype, we ranked risk and protective genes based on importance using “the least absolute shrinkage and selection operator (LASSO) regression analysis” via the “glmnet” R package. “Multivariate logistic regression analysis” was used for constructing a stemness-based classifier based on the following formula: riskscore = Coefficient (stemness gene i) * Expression (Stemness gene i). Next, the classifier was normalized to the 0-1 range. We calculated the normalized riskscore using the “survminer” package, and the patients were categorized based on the threshold value of the risk score into the high-risk group (HRG) and the low-risk group (LRG). Finally, we used the “time-dependent receiver operating characteristic (ROC) curve” to determine the predictive ability of the stemness-based classifier model, and the area under the ROC curve (AUC) value was determined using the “pROC” package. We also analyzed the differences in tumor immune microenvironment (TIME) and the pathways enriched in patients in the HRG and LRG.

We obtained data from the “Cancer Cell Line Encyclopedia” database and used the “pRRophetic” software to predict the patient’s response to erlotinib by determining the half-maximal inhibitory concentration (IC₅₀) value of the drug for all patients (15, 16). The IC₅₀ value is the measure of the drug efficacy in suppressing a specific biological function in patients with HCC. Therefore, a lower IC₅₀ value indicates higher sensitivity of patients to a specific drug. We used the “ridge regression” model on the expression data of patients with HCC to predict the response of patients, and the accuracy of the prediction was evaluated using 5-fold cross-validation.

Human hepatoma cell lines including MHCC-97H, PLC/PRF/5, Huh7, HepG2, and Hep3B as well as Hepa1–6 cell line (C57BL/6‐derived hepatoma) were all commercially obtained from Cell Bank of Type Culture Collection of the Chinese Academy of Sciences. DMEM or RPMI‐1640 media (Gibco, Thermo Fisher Scientific, Waltham, MA, USA) added with 1% penicillin and streptomycin (Solarbio, Beijing, China) and 10% FBS (Biological Industries, CT, USA) was used to grow the above mentioned cell lines, which were incubated in a 5% CO₂ incubator at 37 °C. With the use of Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA), the sequences of human UCK2 siRNAs were transfected into hepatoma cell lines.

A total of 5x10⁶ Hepa1-6 cells in Matrigel were subcutaneously implanted into 5-week-old male C57BL/6 mice (Beijing HFK Bioscience) via the right flank to develop immunocompetent murine models of HCC. When tumors grew to a volume of 200-300 mm³, the mice were given PBS or Anti-PD1 (Hamster anti-murine PD1 mAb J43 (BioXCell) i.p. at 10 mg/kg every 3 days for a total of 5 doses) in a randomized manner. Finally, the mice tumors were collected for further studies.

We performed RNA isolation and qRT-PCR as described previously (17). Supplementary Table S1 shows the primer sequences. The study protocol was approved by the Ethics Committee of Peking University Shenzhen Hospital. Informed consent was obtained from all patients included in the study.

We used the “Pearson correlation” test to determine the correlation between two variables that are not linearly correlated. The categorical and pairwise characteristics of different subgroups were compared using the “Chi-squared test”. The data of skewed distribution and ordinal data in the subgroups were compared using the “Wilcoxon test”. Student t test was used to compared between two groups. We performed a “KM survival curve” analysis to determine the OS of patients with HCC for a specific period. The significant difference in survival was determined using the “log-rank” test. We used the “R (version 4.2.0)” software for statistically analyzing the data. P < 0.05 (two-tailed) was considered statistically significant.

First, we used the “OCLR” algorithm on the transcriptomic data of patients in TCGA-LIHC cohort to obtain the mRNAsi scores. Next, we categorized patients into the mRNAsi-high and mRNAsi-low groups based on the optimal cutoff value calculated using the “survminer,” Next, we performed the “KM survival curve analysis” to determine the impact of mRNAsi on the OS of patients with HCC. The KM plot showed that the OS of patients in the mRNAsi-high group was poor compared to the mRNAsi-low group (). Next, we ranked all patients from low to high mRNAsi scores and determined the correlation between the mRNAsi scores and the clinicopathological features of patients. A significant difference in gender, survival status, and tumor grade was observed in patients in both mRNAsi groups (). Additionally, the mRNAsi scores of patients with stage II HCC were significantly higher compared to stage I HCC; however, no such trend was observed in patients with more advanced clinical stages (). Moreover, the OS of patients with stage I and II HCC in the mRNAsi-high group was poor compared to the mRNAsi-low group (). We observed a similar trend in patients with stage III and IV HCC; however, the difference was not significant (). There was a trend between mRNAsi score and tumor grade (). The OS of patients with grade I and II HCC in the mRNAsi-high group was better compared to the mRNAsi-low group (). Similarly, the patients with stages III and IV exhibited a significant trend (). 1 1 1 1 1 1 1 1

We constructed a novel stemness-based classifier using multiple bioinformatic tools to investigate the mechanism of the two mRNAsi groups and the functions of the mRNAsi model. First, we identified DEGs in tumor and adjacent tissues of patients in TCGA-LIHC and GSE14520↗ cohorts based on the following criteria: “∣FC∣> 1.2” and “P < 0.05.” Next, we performed the correlation analysis on DEGs and the mRNAsi scores using “Pearson’s correlation coefficient” based on the following criteria “∣R∣> 0.3” and “P < 0.05.” We identified 19 mRNAsi-related DEGs, of which 13 genes were upregulated (risk genes, HR > 1) and six genes were downregulated (protective genes, HR < 1) in patients in the mRNAsi-high group (Figure 1A).

We determined the heterogeneity of stemness characteristics by constructing a novel stemness-based classifier based on 19 mRNAsi-related DEG expressions using “unsupervised consensus clustering.” The patients were categorized into the stemness subtype C1 (109 cases, 29.9%), the stemness subtype C2 (138 cases, 37.8%), and the stemness subtype C3 (118 cases, 32.3%). Figures 1B–D show a significant difference in the expression pattern of 19 mRNAsi-related DEGs. Furthermore, the “KM survival curve” analysis showed that the patients in C3 had a significantly poor prognosis compared to patients in C1 and C2 (Figure 1E) in TCGA-LIHC cohort. Similar results were observed in the GSE14520↗ cohort (Figure 1F). The mRNAsi scores of patients in C3 were higher compared to C1 and C2 (Figure 1G). Figure 1H shows an increase in JPT1, CEP131, SNRPA1, CKS1B, SNRPE, CCT3, MSTO1, UCK2, MRTO4, SRM, PUS1, RRP1, and RPL8 expression (risk genes) in patients in C3. Furthermore, an increase in DNASE1L3, CCL14, GPR182, GHR, ANO1, and KLRB1 (protective genes) expression was observed in patients in C3 (Figure 1H). These results showed that the stemness-based classifier could predict the patient’s prognosis and should be analyzed further.

We determined the clinical features of patients in the three stemness subtypes. The percentage of younger patients in C3 was high (). The results revealed no significant difference in the gender of patients from TCGA-LIHC cohort in the three stemness subtypes (). Additionally, a significant difference in the clinical features, such as tumor stage, grade, and stemness, was observed in patients in the three stemness subtypes (). The patients with high grades were classified in C3, which could be the cause of the poor prognosis of patients in C3 (). 2 2 2 2

DNA damage assessment, including aneuploidy, homologous recombination deficiency (HRD), the fractions of genome altered, and the number of segments, is associated with malignancy and immune cell infiltration. Aneuploidy occurs in various cancers and indicates an increased evasion of immune surveillance by cancer cells and reduced patient response to immunotherapy (18). The aneuploidy score of patients in C3 was significantly higher compared to C1 and C2 (Supplementary Figure 3A). Furthermore, the incidences of HRD were higher in patients in C3 compared to C1 and C2 (Supplementary Figure 3B). Furthermore, the fractions of genome altered were lower in patients in C1 and C2 compared to C3 (Supplementary Figure 3C). The number of segments in patients in C1 was low compared to C2 (Supplementary Figure 3D). A significant difference was observed in the rate of TMB in patients in C3 and C1 (Supplementary Figure 3E). We also compared our signature to a previously published classifier. We obtained the HCC subtypes from a previous study (19). The study showed that the prognosis of patients in C3 was better (19), and these patients were primarily enriched in C1 in our study (Supplementary Figure 3F). TMB score of patients in C3 was significantly higher compared to C1. The rate of mutations in TP53, CSMD3, LRP1B, and DNAH7 was lower in patients C3 compared to C1 (Supplementary Figure 3G). The rate of mutations in TP53 (50% vs. 24%, 11%), CSMD3 (14% vs. 6%, 4%), LRP1B (9% vs. 12%, 2%), and DNAH7 (10% vs. 7%, 2%) was higher in patients in C2. Therefore, the degree of malignancy in patients in C3 was higher.

We analyzed the abundance of immune cell infiltration in patients from TCGA-LIHC cohort using CIBERSORT (Figure 2A). A difference in the abundance of macrophages was observed in patients in these three stemness subtypes. A significant increase in the abundance of M0 and a significant decrease in M1 and M2 macrophage abundance was observed in patients from TCGA-LIHC cohort in C3. Additionally, the results revealed a significant decrease in the abundance of resting mast cells in patients in C3 compared to C1. A significant decrease in the abundance of resting memory CD4 cells, which are activated T cells, was observed in patients from TCGA-LIHC cohort in C3. Boxplots were used to explore the relationship between the three immune cell clusters based on ssGSEA and mRNAsi scores. The results showed low stromal scores and a high immune score in patients in C3 (Figure 2B). We used ssGSEA to assess 29 immune gene signatures concluding the primary function and immune, stromal, and other cell (20). The infiltration of pro and anti-tumor immune cells was high in patients in C3. Moreover, the proliferation of cancer cells in patients in C3 was high (Figure 2C). The pro and anti-tumor immune score of patients in C3 was high. Finally, we performed PROGENy analysis (21) to identify enriched signaling pathways. The results revealed that the Trail, TGF-β, NFKB, TNFA, MAPK, HYPOXIA, P53, and EGFR signaling pathways were enriched in patients in C3, and the VEGF and PI3K signaling pathways were enriched in patients in C1 (Figures 2D, E).

The results showed an increase in therapeutic signatures, such as base excision repair, cell cycle, the Fanconi anemia pathway, DNA replication, miRNA in cancer, nucleotide excision and mismatch repair, homologous recombination, the p53 signaling pathway, oocyte meiosis, oocyte maturation mediated by progesterone, pyrimidine metabolism, proteasome, spliceosomes, antiviral carcinogenesis in patients in C3 compared to C1 and C2 (Figure 3A). An increase in the WNT-β-catenin signaling pathway, IDH1, and KDM6B expression was observed in some therapeutic signature scores. We also found that the radiotherapy-related signatures were significantly differentially among hypoxia, cell cycle, and DNA replication. Based on a previous study (22), patients in C3 could benefit from anti-PPARG therapy instead of anti-EGFR and anti-FGFR3 therapies (Figure 3B). Immunotherapy is widely used for treating patients with HCC. Hence, we used the HisgAtlas dataset for analyzing immune checkpoint genes. The results revealed a significant increase in immune checkpoint gene expression in C2 compared to C1. An increase in the expression of several immune checkpoint genes, such as VSIR, PDCD1LG2, CD274, HAVCR2 (TIM-3), BTLA, LAG-3, TIGIT, CTLA-4, and PDCD1, was observed in patients in the high-risk subgroups of C3 (Figure 3C). Furthermore, we used the “Tumor Immune Dysfunction and Exclusion (TIDE)” software to determine the response of patients in the two clusters to immunotherapy. Patients with higher TIDE scores could evade immune surveillance, which indicates that these patients may not benefit from immunotherapy. The TIDE scores of patients from TCGA-LIHC cohort in C3 were higher compared to C1 or C2, which suggests that patients in C3 could evade immune surveillance and may not benefit significantly from treatment with immunotherapy (Figure 3D). There was an increase in IFN-β, T cell exclusion, and dysfunction, and MDSC scores, and low T cell dysfunction scores in patients in C3 (Figure 3D). This could be the underlying cause of the high degree of malignancy in patients in C3.

First, we used the “limma” package to identify DEGs in these three stemness clusters. We identified 166 upregulated and 303 downregulated genes in patients in C1 compared to C2 and C3 (). No DEGs were observed in patients in C2 compared to C1 and C3. However, 524 upregulated and 242 downregulated genes were observed in patients in C3 compared to C1 and C2 (). Next, we performed the “Gene ontology” and “KEGG pathway enrichment” analyses to annotate gene functions. The results revealed that these 166 upregulated genes in C1 were enriched in the protein activation pathway and high-density lipoprotein particle (). Further, 242 downregulated genes in C3 were enriched in the small molecule catabolic processes and high-density lipoprotein particles. The “KEGG pathway enrichment” analysis showed the enrichment of the retinol metabolism, the biosynthesis of the steroid hormones, and the metabolism of xenobiotics by cytochrome P450 in C1 or C3 (). 4 4 4 4

We identified 816 DEGs in three stemness subtypes and performed “COX regression analysis” to identify the prognostic factors. A total of 272 risk genes and 26 protective genes associated with the OS were identified (Figure 4A, P < 0.001). We categorized patients based on the expression of mRNAsi-related DEGs and heterogeneity of stemness characteristics using unsupervised consensus clustering. The “glmnet” package was used to perform the “LASSO-COX regression” analysis. Figure 4B shows a gradual decrease in the regression coefficients of predictors towards zero with an increase in lambda. We determined the optimal lambda value using 10-fold cross-validation, and the results showed that the partial likelihood deviance of the model was minimized when lambda = 0.0433 (Figure 4C). We selected 20 genes associated with this lambda value for further analysis. The model reached the best when the lambda = 0.0433. We categorized patients into stemness subtype I (193 patients, 52.9%) and stemness subtype II (172 patients, 47.2%, Figures 4A-C). We calculated the riskscore as follows = 0.23 * UCK2 + 0.138 * HMMR + 0.094 * CABYR + (-0.072 * CFHR3) + (-0.148* PPARGC1A) + 0.129 * EPO + (-0.094 * FTCD) + 0.071 * CXCL8 + 0.069 * MAGEA6 (Figure 4D). A “time-dependent ROC” analysis was used to determine the accuracy of the model in predicting the prognosis. The AUC values of 1-, 2-, 3-, 4-, and 5-year OS were 0.8, 0.75, 0.76, 0.78, and 0.76, respectively, indicating the good predictive significance of the model (Figure 4E). The KM survival curve showed that the prognosis of patients in stemness subtype II was poor compared to patients in stemness subtype I (Figure 4F). Finally, we validated the model in GSE14520↗, and the results were similar (Figures 4G, H).

A significant difference in the gender, T stage, tumor stage, grade, stemness, and previous miRNAsi model was observed in patients in HRG and LRG (). The results showed that mRNAsi was positively correlated with the degree of risk score (). The ssGSEA method of GSVA package was used to calculate the scores of the stemness related gene sets downloaded from the molecular signature database (MSigDB) in the TCGA-LIHC cohort samples, including “WONG EMBRYONIC STEM CELL CORE”, “YAMASHITA LIVER CANCER STEM CELL UP” and “YAMASHITA LIVER CANCER STEM CELL DN”, and then correlation analysis was performed with nine-gene signature. We found that the wone embryonic stem cells core score and yamashita liver cancer stem cell up score of HCC samples were both significantly positively correlated with the risk score, respectively. While yamashita liver cancer stem cell down score was significantly negatively correlated with risk score (). We also found that the high-risk group had a significantly higher wone embryonic stem cells core score and yamashita liver cancer stem cell up score, as well as a significantly lower yamashita liver cancer stem cell down score than the low-risk group (). In addition, the KM survival curve showed that the prognosis of patients with high riskscores was poor compared to patients with riskscores in both mRNAsi groups (). Moreover, the OS of patients with stage 1 and II HCC in HRG was better compared to LRG (). A similar survival trend was observed in patients with stage III and IV HCC; however, the difference was not significant (P = 0.055). The OS of patients with grade I and II HCC in HRG was better compared to LRG. A similar survival trend was observed in patients with stage III and IV HCC (). These results show that our stemness-based classifier could predict the patient’s prognosis and should be analyzed further. 5 5 6 6 5 5 5

To investigate the difference between the riskscore and immune characteristics, we compared the 22 immune cell types in the TIME. The infiltration of several immune cells differed in patients in HRG and LRG (Figure 5A). The differences in the immune cell infiltration in patients in the two stemness subtypes were analogous to the previous mRNAsi model. A difference in the abundance of macrophages was observed, wherein a significant increase in the abundance of M0 macrophages and a significant decrease in M1 and M2 macrophages in patients from TCGA-LIHC cohort in HRG (Figure 5A). The stromal score of patients in HRG was remarkably increased in LRG. Furthermore, the results revealed a significant increase in the immune scores of patients in HRG compared to LRG. However, no significant difference was observed in the ESTIMATE scores of patients in both riskscore groups (Figure 5B). Next, we analyzed the riskscore and pathway activity. The riskscore was significantly associated with the VEGF and PI3K signaling pathways, etc. (Figure 5C). Moreover, several immune cells were significantly associated with the riskscore (Figure 5D). Additionally, the riskscore was negatively associated with angiogenesis, endothelium, cancer-associated fibroblasts (CAFs), and EMT signature (Figure 5E). All genes from the HALLMARK database were subjected to GSEA to analyze the enrichment of the riskscores. The results revealed enrichment of the mitotic spindle, G2M checkpoint, protein secretion, unfolded protein response, and mTORC1 signaling by genes in HRG (Figure 5F).

First, we determine immune checkpoint gene expression in patients in HRG and LRG. The results showed an increase in PDCD1, CD274, CTLA4, LAG3, PDCD1LG2, BTLA, HAVCR2, TIGHT, and VISIR expression in patients in HRG (Figure 6A). Furthermore, there was a strong association between riskscore and immune checkpoints expression (Figure 6B). Next, we used the TIDE database to determine the response of patients in different riskscore groups to immunotherapy (Figure 6C). A significant difference was observed in the TIDE scores and MDSC abundance in patients in LRG and HRG. A significant correlation was observed between the degree of riskscore and the infiltration of MDSC, CAFs, IFNG, and the TIDE scores. Additionally, the degree of riskscore was negatively correlated with T cell dysregulation (Figure 6D). Cyclopamine, AZ628, sorafenib, VX-680, and imatinib are used to treat cancer patients. The results revealed that patients with high riskscore could benefit from these drugs (Figure 6E). Additionally, to explore the correlation between riskscore and immunotherapy, we determined the ability of riskscores to predict the patient’s response to ICB. GSE202069↗ dataset included a total of 41 samples, of which 17 patients received PD-1 therapy, including 9 non-responders and 8 responders. Riskscore was calculated in GSE202069↗ dataset using the same method as in TCGA, and the difference in riskscores between responders and non-responders was compared. The results showed that the risk score of the non-responder group was higher and the prognosis was worse (Figures 6F, G). The survival curve between the risk groups showed that the prognosis was worse in the high-risk group (Figure 6H). ERP117672 included a total of 40 RNA-seq samples from HCC patients treated with pabolizumab, including 29 non-responder samples, 6 responder samples, and 5 non-evaluable samples. The results showed that Riskscore was higher in PD and SD group than PR group (Figure 6I).

To validate the robustness of the nine gene signature, we collected tumor tissues and matched non-tumor tissues from 19 patients with HCC and determined the expression of nine genes using qRT-PCR. A significant increase in UCK2, HMMR, and MAGEA6 expression and a significant decrease in CXCL8, EPO, PPARGC1A, FTCD, and CFHR3 expression was observed in tumor tissues. In addition, an increase in CABYR expression was observed in tumor tissues; however, the increase was not statistically significant (Figure 7A). Because previous analyses have shown that gene signature could predict immunotherapy outcomes. Herein, we performed an anti-PD1 treatment in an HCC mice model. Macroscopically, smaller size and fewer numbers of HCC nodules were found in the mice with anti-PD1 treatment compared to the untreated mice (Figure 7B). Significant decreases of tumor UCK2 expression was found in the anti-PD1 treated mice compared to the untreated mice (Figure 7C). Furthermore, we knockdown UCK2 expression using siRNA in 5 HCC cell lines, but only effective reduction in MHCC-97H and Hep3B cells (Figure 7D). To understand whether UCK2 regulates the expression of immune genes, we determined immune genes expressions in MHCC-97H and Hep3B cells treated siUCK2. Those immune genes were reported to regulated HCC immune infiltration or immunotherapy (23–27). RT-qPCR results showed that siUCK2 inhibited HMGB1 expression in Hep3B cells (Figure 7E), and significant reduction of PGAM1, PHF19 and upregulation of HMGB1, PRRC2A in MHCC-97H cells (Figure 7F).