Pan-Cancer Analysis of Genomic and Prognostic Characteristics Associated With Coronavirus Disease 2019 Regulators

Background: Cancer patients are alleged to have poor coronavirus disease 2019 (COVID-19) outcomes. However, no systematic or comprehensive analyses of the role and mechanisms of COVID-19 receptor-related regulators in cancer are available. Methods: We comprehensively evaluated the genomic alterations and their clinical relevance of six COVID-19 receptor-related regulators [transmembrane serine protease 2 (TMPRSS2), angiotensinogen (AGT), angiotensin-converting enzyme 1 (ACE1), solute carrier family 6 member 19 (SLC6A19), angiotensin-converting enzyme 2 (ACE2), and angiotensin II receptor type 2 (AGTR2)] across a broad spectrum of solid tumors. RNA-seq data, single nucleotide variation data, copy number variation data, methylation data, and miRNA–mRNA interaction network data from The Cancer Genome Atlas (TCGA) of 33 solid tumors were analyzed. We assessed the sensitivities of drugs targeting COVID-19 receptor-related regulators, using information from the Cancer Therapeutics Response Portal database. Results: We found that there are widespread genetic alterations of COVID-19 regulators and that their expression levels were significantly correlated with the activity of cancer hallmark-related pathways. Moreover, COVID-19 receptor-related regulators may be used as prognostic biomarkers. By mining the genomics of drug sensitivities in cancer databases, we discovered a number of potential drugs that may target COVID-19 receptor-related regulators. Conclusion: This study revealed the genomic alterations and clinical characteristics of COVID-19 receptor-related regulators across 33 cancers, which may clarify the potential mechanism between COVID-19 receptor-related regulators and tumorigenesis and provide a novel approach for cancer treatments.


INTRODUCTION
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus has resulted in the ongoing coronavirus disease 2019 (COVID- 19) pandemic. As of July 7, 2021, there are 184,820,132 confirmed cases and 4,002,209 deaths, with the numbers still surging worldwide (1). With the continued increase in cases and affected regions, patients with chronic conditions, such as cancer, have been disproportionately affected (2)(3)(4)(5). The COVID-19 pandemic has been identified as a global health emergency by the World Health Organization (WHO).
Respiratory inflammation is activated by the reninangiotensin-aldosterone system (RAS), which maintains the blood pressure by angiotensin II (Ang II) and is catalyzed by the angiotensin-converting enzyme (ACE). Angiotensinogen (AGT) is the protein precursor of Ang II (6,7). Ang II receptor type 2 (AGTR2), a member of the G-protein coupled receptor 1 family, functions as a receptor for Ang II. ACE2 degrades Ang II, counteracting its chronic effects, and serves as the SARS-CoV-2 receptor. ACE2 is also a molecule present on the surface of various cell types, including type II alveolar cells, bronchial transient epithelial secretory cells, endothelial cells, intestinal epithelium cells, and uterine epithelial cells (8). The spike protein (S protein) of SARS-CoV binds to cell surface ACE2 receptors (9). ACE1, homologous to the ACE2 gene, may be involved in the progression of diseases caused by several human coronaviruses (10,11). Transmembrane serine protease 2 (TMPRSS2), a member of the serine protease family, facilitates human coronavirus infections (SARS-CoV and SARS-CoV-2) via proteolytic cleavage of the ACE2 receptor, which promotes viral uptake and cleavage of coronavirus spike glycoproteins, activating glycoproteins for host cell entry (12)(13)(14). Solute carrier family 6 member 19 (SLC6A19), a SARS-CoV-2 co-receptor, is a neutral amino acid transporter and forms a heterodimer with ACE2 (15). However, the genomic alterations and prognostic characteristics of COVID-19 receptor-related regulators in cancer are still unclear.
The clinical symptoms of COVID-19 range from asymptomatic to severe cardiopulmonary disease (16)(17)(18). Enhanced expression of ACE2 and immunosuppressive states caused by malignancies and anticancer treatments, such as chemotherapy or surgery, contribute to more severe disease in older patients with COVID-19 (19,20). Recent studies also identified that aberrant expression of ACE2 receptor-related regulators is associated with the activation of several cancerassociated pathways (21)(22)(23). Therefore, it is of great clinical significance to clarify the genomic and clinical characteristics of the six ACE2 receptor-related regulators among 33 solid tumors for the management and treatment of tumor patients with COVID-19.

Dataset Acquisition and Preprocessing
The Genotype-Tissue Expression (GTEx) dataset (V7.0) (https:// commonfund.nih.gov/GTEx/) was used for gene expression analysis in normal tissues from healthy individuals. The tumorassociated data are composed of mRNA Seq data, clinical data, single nucleotide variation (SNV) data, copy number variation (CNV) data, and methylation data, which were collected from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer. gov/). Reverse phase protein array (RPPA) data were obtained from The Cancer Proteome Atlas (TCPA) (https://tcpaportal. org/tcpa/index.html). The Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerrxgene.org) was used to investigate the correlation between gene expression and drug sensitivity.

mRNA Expression Analysis
For mRNA differential expression analysis between paired tumor and normal samples, TCGA mRNA expression was normalized using RNA-Seq by Expectation-Maximization (RSEM). The number of samples for each cancer type ranged from 48 to 1,098, where only 14 cancer types that harbored over 10 pairs of tumor and normal samples were incorporated into analyses, namely, BLCA, BRCA, COAD, ESCA, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, PRAD, STAD, and THCA. Gene expression values were represented as RNA-Seq by Expectation-Maximization (RSEM) normalized data (24). The genes with a fold change (FC) <2 and significance false discovery rate (FDR) <0.05 underwent further analysis.

Subtype Analysis
Expression subtype analysis was used to find clinically relevant genes that may affect cancer subtype. To make the analysis feasible, the number of subgroups in a given subtype was at least 10, leaving 11 cancer types for gene analysis. We analyzed 11 cancer types for ACE2 receptor-relevant genes using a Student's t-test (n_subtype = 2) and ANOVA test (n_subtype > 2). The method used for the clinically relevant analysis depends on the number of subgroups in each cancer subtype.

Survival Analysis
For expression survival analysis, mRNA expression and clinical survival data were merged by the sample barcode. Tumor samples were divided into "high" and "low" gene expression groups using the median RSEM value. The R package "survival" was used to fit the survival time and survival status for the two groups. A Cox Proportional-Hazards model was calculated for each gene using the R package. Genes that had a Kaplan-Meier log-rank test p-value <0.05 were retained.

SNV Analysis
SNV data of 33 cancer types (N = 8663) were investigated. SNV oncoplot (or waterfall plot) was generated by maftools (25). The TCGA SNV data includes the following variant type values: Missense_Mutation, Silent, 5' Flank, 3' UTR, RNA, In_Frame_Del, Nonsense_Mutation, Splice_Site, Intron, 5' UTR, In_Frame_Ins, Frame_Shift_Del, Nonstop_Mutation, 3' Flank, Frame_Shift_Ins, and Translation_Start_Site. The Silent, Intron, IGR, 3' UTR, 5' UTR, 3' Flank, and 5' Flank were filtered out for SNV percentage calculation. The percentage of SNVs in each gene's coding region was calculated by the number of mutated samples divided by the number of cancer samples. SNV data and clinical overall survival data were combined, and the R package was used to estimate the survival difference between mutated and non-mutated genes.

CNV Analysis
CNV raw data from 33 cancer types (N = 11,495) were investigated and processed with GISTICS2.0 (26). The CNV was divided into heterozygous and homozygous CNV subtypes, which represented the occurrence of CNV on one chromosome or two chromosomes, respectively. The homozygous or heterozygous CNV profile showed the percentage of homozygous or heterozygous CNV, including CNV amplification and deletion percentages for each gene in each cancer. The percentage of CNV subtypes was calculated using GISTIC-processed CNV data. Only genes with >5% CNV were considered significant. As the method has been employed by Schlattl et al. (27), the mRNA expression and CNV data were merged by a sample's TCGA barcodes. The association between paired mRNA expression and CNV percentage were detected based on a Pearson product-moment correlation coefficient and t-distribution.

Methylation Analysis
Methylation data of paired tumor and normal samples across 14 cancer types (N = 10,129) were investigated. The mRNA expression and methylation data were merged by a sample's TCGA barcode. The association between paired mRNA expression and methylation was tested based on a Pearson product-moment correlation coefficient and t-distribution. The mRNA expression and methylation data of the regulators were merged via the TCGA barcode. The association between paired mRNA expression and methylation data was calculated using the Pearson's product-moment correlation coefficient, followed by a t-distribution test. p-values were adjusted by the FDR, and genes with an FDR ≤ 0.05 were retained. Further analysis was carried out on genes that were significantly influenced by genome methylation. Methylation data and clinical overall survival data were combined, and the methylation level of a gene was divided into two groups by median methylation. Cox regression was performed to estimate the hazard (risk of death). If the Cox coefficient was <0, the high methylation group showed a poorer survival, the Hyper_worse defined as High risk, otherwise defined as Low risk.

Pathway Activity Analysis
Following the method used by Ye et al. (28), RPPA data from TCPA were used to calculate a score for 7,876 samples. RPPA data of replicates-based normalization (RBN) were median-centered and normalized by the standard deviation across all samples for each component to obtain the relative protein level. The pathway score is the sum of the relative protein levels of all positive regulatory components minus the sum of the relative protein levels of all negative regulatory components in a given pathway. Gene expression was divided into two groups (upregulation group or downregulation group) by the median expression. The difference in the pathway activity score (PAS) between the two groups was determined. When PAS (gene A, upregulation group) was greater than the PAS (gene A, downregulation group), we considered gene A as having an activating effect on a pathway; otherwise, gene A had an inhibitory effect on a pathway.
miRNA Regulation Network Analysis miRNA regulation data (N = 9,105) was collected from TCGA across 33 cancer types. miRNA expression and gene expression were merged by TCGA barcode. The association between paired mRNA and miRNA expression was tested based on a Pearson product-moment correlation coefficient and t-distribution. The p-value was adjusted by the FDR, and genes with an FDR of ≤0.05 (R < 0) were retained. The correlation was calculated for all paired samples. In addition, with consideration to the presence of positive regulators (including transcription factors), an miRNAgene pair with negative correlation was considered as a potential negative regulation pair. Network was generated by visNetwork R packages.

Drug Sensitivity Analysis
Following the method used by Rees et al. (29), 481 small molecules from the Cancer Therapeutics Response Portal (CTRP) were collected. To analyze the correlation between gene expression and drug sensitivity, the values from the area under the dose-response curve (AUC) for drug and gene expression profiles for all cancer cell lines were downloaded. The Pearson correlation coefficients of transcription levels and AUCs were normalized using Fisher's z transformation. The Pearson correlation coefficients of the transcript levels and AUCs were normalized using Fisher's z transformation. A Bonferronicorrected, two-tailed distribution test, with a family-wise error rate of <0.025 in each tail, was used for the z-score calculation. Pearson correlation coefficients of annotated drug-target pairs were compared with the same number of correlation pairs generated by random sampling of the correlations. The gene set drug resistance analysis was performed on IC 50 drug data.

Statistical Analysis
Correlations between gene expression were evaluated using the Spearman's correlation test. The prognostic significance of the indexes was estimated using Kaplan-Meier survival curves and compared by a log-rank test. The Cox proportional hazards model was used to calculate the adjusted hazard ratio (AHR). All statistical analyses were performed with SPSS version 23.0 (SPSS Inc, Chicago, IL, USA) and R version 3.4.4 (http://www.r-project.org). p < 0.05 was considered as statistical significant.
We further explored the effect of regulator expression on cancer survival and found that high expression of TMPRSS2 in ACC; ACE1 in UVM and LUSC; AGTR2 in KIRP, KICH, and LUSC; ACE2 in LGG; and AGT in UVM were associated with poor survival of cancer patients, while expression of TMPRSS2 in KIRP, KICH, PAAD, and LIHC; ACE1 in KIRC, LIHC, and OV; SLC6A19 in KIRC, KIRP, and ESCA; ACE2 in KIRC and UVM; and AGT in PAAD were associated with good survival in cancer patients ( Figure 1C and Supplementary Figure 1; p < 0.05). As shown in Figure 1D, the low expression of ACE2 was significantly associated with poor survival in KIRC (HR = 0.54; p = 6.4e−05). These results indicated that the expression of COVID-19 receptor-related regulators may play an important role in the progression and deterioration of cancer with COVID-19.

CNV of ACE2 Receptor Regulators
To identify the CNV change of ACE2 receptor regulators at the chromosome arm level, we analyzed the CNV data of ACE2 receptors from the TCGA database. We found that TMPRSS2, SLC6A19, ATGR2, AGT, ACE2, and ACE1 had >5% CNV amplification or deletion in 33 cancers. As shown in the CNV pie distribution in Figure 3A, TMPRSS2 had 80% heterozygous amplification in TGCT but 63% heterozygous deletion in ESCA; ACE2 had 51% heterozygous amplification in ACC and >50% heterozygous deletion in OV and KICH; and AGT in LUAD, UCS, BRCA, LIHC, CESC, LUSC, SKCM, ESCA, and CHOL; ACE1 in KIRP; and SLC6A19 in ACC and LUSC had almost 50% heterozygous amplification, whereas ACE1 in KICH; SLC6A19 in TGCT; and AGT in KICH had almost 50% heterozygous deletion. To identify the heterozygous/homozygous CNV profile in each cancer, we further analyzed heterozygous/homozygous amplification and heterozygous/homozygous deletion. As shown in Figure 3B, all regulators had heterozygous amplification and deletion. However, homozygous CNV analysis showed that SLC6A19 had homozygous amplification in 12 solid cancers, with TMPRSS2 homozygous deletion only found in PRAD ( Figure 3C).
Comparing the relationship between CNV and mRNA expression, the correlation analysis indicated that mRNA expression of each regulator was positively correlated with its CNV in most cancers (p < 0.05). However, the expression of TMPRSS2 in KIRC; SLC6A19 in ESCA and PAAD; and ACE1 in TGCT, SKCM, and LIHC were negatively correlated with the CNV (p < 0.05) (Figure 3D). These results indicated that the CNV of ACE2 receptor-related regulators mediated their abnormal expression, which may play an important role in cancer patients with COVID-19.

Methylation Analysis of ACE2 Receptor Regulators
We explored the methylation analysis of ACE2 receptor regulators to identify the corresponding epigenetic methylation levels. As shown in Figure 4A    good survival (p < 0.05; Figure 4C). As shown in Figure 4D, the hypermethylation of SLC6A19 was significantly associated with poor survival in KIRP (p = 8.6e−05; Figure 4D).

Pathway Activity Analysis
The pathway relation network indicated that ACE2 receptorrelated regulators were involved in TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, hormone ER, hormone AR, EMT, DNA damage response, cell cycle, and apoptosis pathways ( Figure 5A). The global percentage of cancers in which regulators have an effect on a pathway showed that ACE1 was involved in the activation of apoptosis, DNA damage, epithelialmesenchymal transition (EMT), hormone ER, hormone AR, RAS/MAPK, and RTK pathways and the inactivation of the cell cycle and TSC/mTOR pathways. ACE2 was associated with the activation of PI3K/AKT, RAS/MAPK, and TSC/mTOR pathways, and with the inactivation of the cell cycle, DNA damage, EMT, and hormone AR pathways. AGT was associated with the activation of the EMT pathway and the inactivation of apoptosis. AGTR2 was associated with the activation of RAS/MAPK, RTK, and TSC/mTOR pathways and the inactivation of apoptotic, cell cycle, DNA damage, hormone ER, and hormone AR pathways. SLC6A19 was involved in the activation of RTK and hormone AR pathways and the inactivation of hormone ER and TSC/mTOR pathways. TMPRSS2 was associated with the activation of the RTK pathway and inactivation of EMT ( Figure 5B). As ACE2 receptor-related regulators are often mutated in UCEC, we further analyzed the global percentage of pathway activity in UCEC. We found that ACE1 was mostly involved in the inhibition of the cell cycle (21% inhibition vs. 7% activation) and activation of RAS/MAPK (9% inhibition vs. 13% activation). ACE2 was mainly involved in the inhibition of hormone AR (12% inhibition vs. 7% activation) and activation of the RTK pathway (0% inhibition vs. 19% activation). AGT was associated with inhibition of apoptosis (18% inhibition vs. 0% activation) and activation of EMT (6% inhibition vs. 16% activation). TMPRSS2 was mainly involved in the inhibition of the DNA damage response (12% inhibition vs. 7% activation) and EMT (34% inhibition vs. 4% activation), while it was associated with the activation of hormone AR (9% inhibition vs. 13% activation), hormone ER (9% inhibition vs. 13% activation), and RTK (6% inhibition vs. 22% activation) pathways (Supplementary Figure 5). These results indicated that the abnormal expression of ACE2 receptor-related regulators
the expression of ACE2 (p < 0.05); hsa-miR-183-5p and hsa-miR-377-3p were negatively regulated with the expression of SLC6A19 (p < 0.05); and hsa-miR-24-3p were negatively regulated with the expression of ACE1 (p < 0.05). These results indicated that the miRNA regulation network mediated ACE2 receptor-related regulators, which may be involved in the progression of cancer in patients with COVID-19.

Drug Sensitivity Analysis
Genomic aberrations influenced clinical response to treatment and are potential biomarkers for drug screening in cancer.
To know the role of ACE2 receptor-related regulators on chemotherapy or targeted therapy, drug sensitivity and gene expression profiling data of cancer cell lines from the CTRP were integrated. Spearman's correlation analysis showed that drug sensitivity toward vincristine, teniposide, ouabain, docetaxel, doxorubicin, erlotinib, afatinib, AZD7762, and AT13387 correlated with the expression of AGTR2, SLC6A19, ACE2, and TMPRSS2 (negative correlation with IC 50 ). Drug resistance toward staurosporine correlated with the expression of TMPRSS2, JW55, FGIN-1-27, BRD-K96431673, BRD-K86535717, BRD-K75293299, and BRD-K49290616 (positive correlation with IC 50 ) (Figure 7). These results indicated that the abnormal expression of ACE2 receptor-related regulators may mediate sensitivity to chemotherapy and targeted drug therapy.

DISCUSSION
COVID-19 is a global health emergency problem with a large number of confirmed cases and deaths that are much greater than any infection in recent decades. Condition severity and mortality have been identified as being significantly higher in patients with other comorbidities, such as cancer (30)(31)(32). As care for chronic conditions, such as cancer, still needs to continue during the pandemic, it is necessary for healthcare providers to determine which type of cancer will put patients at a higher risk of exhibiting severe forms of the COVID-19 infection. Patients have also had to balance the risks and benefits of cancer-directed interventions within the context of the added risk of contracting COVID-19. In this study, we comprehensively analyzed the genomic and prognostic characteristics of six COVID-19 receptorrelated regulators, where we found that genetic and epigenetic alterations, and an miRNA network of COVID-19 receptorrelated regulators, led to their abnormal expression, which correlated significantly with the activation of cancer hallmarkrelated pathways and clinical survival. Targeting these COVID-19 receptor-related regulators may be an important method to treat cancer patients with COVID-19.
We firstly explored the genetic alterations and prognoses of these regulators and found that the abnormal expression was associated with clinical prognosis. Our results indicated that there were 14 tumor types that differentially expressed one or more of these regulators. The regulators were highly expressed in normal mucosal epithelial tissue, such as kidney, urinary bladder, and mucocutaneous and gastrointestinal tract, which was consistent with ACE2 protein distribution. This co-expression pattern further validated that the SARS-CoV-2 entry process requires the interaction of these regulators. By tracking the genetic differences in these six regulators, we found that missense mutations were the main mutation type in SNV, with ACE1 having the highest mutation frequency in cancer. In addition, ACE1 mutations in malignancies were distributed across all exons of ACE1 with several hot spot mutation sites. ACE mutations have been reported to be involved in a number of lymph node metastases of gastric cancer (33) and associated with a worse prognosis in prostate cancer (34). However, there was also non-conformity between genomics alternation and clinical prognosis. Thus, we speculated that genetic and epigenetic alteration of the regulators may cause gene dysfunction and promote tumorigenesis in certain contexts.
Further investigation into the biological function of the regulators identified several pathways, including TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, hormone ER, hormone AR, EMT, DNA damage response, cell cycle, and apoptosis pathways, that were significantly enriched in cancers. In UCEC, different ACE2 receptor-related regulators were associated with different cancer-related signaling pathways. For example, TMPRSS2 was involved in the activation of the RTK pathway and AGTR2 was associated with the inhibition of the cell cycle and the apoptosis pathway. Recent studies identified that Ang II can also promote cell growth and proliferation via the transforming growth factorbeta (35), RTK (36), and mTOR pathways (37). Activation of Ang II receptor in cancer cell lines resulted in increased MAPK activation, JAK-STAT signaling, and cell proliferation (38,39). Thus, activation and inhibition of cancer-related signaling pathways mediated by ACE2 receptor-related regulator molecular networks played different roles in tumorigenesis and prognosis.
In clinical applications, dexamethasone, which can reduce inflammation, and remdesivir, which can inhibit viral replication, have been widely used to decrease the mortality in cancer patients with COVID-19 (40,41). There are currently no effective drugs for COVID-19. There is an urgent need for therapeutic interventions, especially for cancer patients with weakened immune systems. Our drug sensitivity analysis identified that ACE2 receptor-related regulator expression levels were also involved in drug sensitivity. Vorinostat is an anticancer histone deacetylase (HDAC) inhibitor and has previously been shown to have anti-fibrotic effects and can reduce the risk of acute respiratory deterioration by upregulating ACE2 expression (42,43). Erlotinib, an epidermal growth factor receptor (EGFR) inhibitor, has been reported to inhibit endocytosis and intracellular trafficking of multiple viruses, including hepatitis C, dengue, and Ebola, exerting broadspectrum antiviral effects by increasing ACE2 expression (44). Thus, we speculate that targeting ACE2 receptorrelated regulators will become an ideal approach in cancer treatment. However, variations of ACE2 receptor-related regulators exist at all regulation levels, including genetics and epigenetic alterations, mRNA expression, miRNA networks, and pathway correlations. These variations may alter drug effects, treatment responses, and patient survival. Thus, the potential mechanisms of each drug's effect on ACE2 receptorrelated regulator expression and cancer progression require further investigation.

CONCLUSION
In conclusion, our findings indicate the need for precautions for and protection of cancer patients during the COVID-19 pandemic. However, the balance between the risks and benefits of cancer-directed interventions should be reassessed. Thus, targeting ACE2 receptor-related regulators could be a promising strategy against cancer patients with COVID-19.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JZ, HJ, and KD performed all experiments, prepared the figures, and drafted the manuscript. JZ, TX, BW, and BC participated in the data analysis and results interpretation. JZ, KD, CC, YY, and JY designed the study and participated in the data analysis. All authors contributed to the article and approved the submitted version.