Dissecting LncRNA Roles in Renal Cell Carcinoma Metastasis and Characterizing Genomic Heterogeneity by Single-Cell RNA-seq

Long noncoding RNAs (lncRNA) have recently emerged as important regulators in cancer cell proliferation and metastasis. However, the role of lncRNAs in metastatic clear cell renal cell carcinoma (ccRCC) remains unclear. Here, single-cell RNA sequencing data were analyzed from primary renal cell carcinoma and paired metastatic renal cell carcinoma specimens, and characterized the expression profiles of over 10,000 genes, including 1,874 lncRNAs. Further analysis revealed that lncRNAs exhibit cancer type– and tissue–specific expression across ccRCC cells. Interestingly, a number of lncRNAs (n = 173) associated with ccRCC metastasis, termed ccRCC metastasis–associated lncRNAs (CMAL). Moreover, functional analysis based on a CMAL-PCG coexpression network revealed that CMALs contribute to cell adhesion, immune response, and cell proliferation. In combination with survival analysis, 12 CMALs were identified that participate in TNF and hypoxia-inducible factor 1 signaling to promote ccRCC metastasis. Further investigation on intratumoral heterogeneity showed that some CMALs are selectively expressed in different subpopulations. Implications: To explore ccRCC metastasis, the current study performed a global dissection of lncRNAs and a complex genomic analysis of ccRCC tumor heterogeneity. The data shed light on the discovery of potential lncRNA biomarkers and lncRNA therapeutic targets.


Introduction
Kidney cancer, or renal cell carcinoma (RCC), is characterized by comprehensive biological process. Contributing factors include aberrant expression of protein-coding genes (PCG) and noncoding genes (1). The most common type of RCC is clear cell renal cell carcinoma (ccRCC). Patients with ccRCC tend to have worse prognosis than patients with other RCC subtypes, and ccRCC is highly prone to metastasize to the lungs, bones, and lymph nodes (2). As one of the major drivers of mortality in patients with cancer, metastasis is a multistep process that involves dissemination of cancer cells to distant organ sites and subsequent adaptation to foreign tissue microenvironments (3). Previous studies have shown some different molecular and signaling pathway signatures between metastatic renal cell carcinoma (mRCC) and primary renal cell carcinoma (pRCC). For example, the expression of cell-surface protein CD44 was significantly higher in mRCCs compared with pRCCs (4); von Hippel-Lindau tumor suppressor (VHL) gene presents variant genotypes between the primary and metastatic tumor (5); a set of genes involved in the process of epithelial-mesenchymal transition (EMT), their expression levels also vary dramatically in aggressive tumors (6). Currently, some significant progress has been made in the development of targeted therapeutic approaches for mRCC. A vascular endothelial growth factor receptor inhibitor, Sunitinib, is now widely prescribed for the treatment of mRCC (7). However, other kidney cancer metastasis markers are not effectively applied to clinical treatment. Because of concerns with drug resistance and side effects, it would be imperative to find new markers for mRCC prognosis and therapy.
LncRNAs are defined as transcripts longer than 200 nucleotides with no protein-coding potential (8). During the past decade, accumulating evidence demonstrates that lncRNAs play important roles in tumorigenesis and cancer invasiveness. For example, overexpression of the lncRNA HOTAIR promotes breast cancer metastasis (9); oncogenic lncRNA PCAT-1 contributes to prostate cancer cell proliferation via regulating MYC stabilization (10). Thus, lncRNAs as functional transcripts could contribute to improve prognosis in patients with cancer (11). Although thousands of lncRNAs have been discovered, their functions in cancer remain poorly studied, and research on lncRNA signatures for cancer metastasis is still fragmentary. LncRNAs were always lower expressed than PCGs, which leads to a difficulty in identification (8). In recent years, RNA sequencing (RNA-seq) has successfully been applied to assess large-scale lncRNAs. In general, the expression of lncRNA is more cancer type-specific than PCGs (12).
The existence of heterogeneity between individual tumors has been recognized by experimental cancer research (13). Intratumor heterogeneity represents a major challenge to effective treatment and personalized medicine (14). The recent development of single-cell RNA-seq makes it feasible to detect genes at single-cell resolution (15). More recently, it was demonstrated that lncRNAs also show expression heterogeneity within glioblastoma cells (16). Therefore, the identification of distinct expression patterns and the dissection of functional roles of lncRNAs may be of vital importance for ccRCC metastasis.
In the present study, we evaluated the differences of lncRNA's expression profiles between primary and metastatic ccRCC. Moreover, a coexpression network was constructed to evaluate the function of these ccRCC metastasis-associated lncRNAs (CMAL). In addition, we found four cell subtypes in metastasis ccRCC single-cell samples and identified specific genes that were able to characterize ccRCC single-cell subtypes. Through further research, we found that that some CMALs were specifically expressed in different subtypes. Moreover, some genes involved in important pathways are also specifically expressed. These findings suggest that research in landscaping metastasis-associated lncRNA may yield to diagnostic accuracy even in therapeutic development.

Renal cell carcinoma RNA-seq datasets
Sequences data and expression data were obtained from three previously RNA-seq datasets: (i) 116 single-cell RNA-seq and 3 popular RNA-seq from NCBI-GEO under accession number GSE73122 (ref. 17 Processing of single-cell RNA-seq data Paired-end 100 bp reads from scRNA-seq were aligned to the human transcriptome (GENCODE GRCh38) using Bowtie2 (version 2.1.0; ref. 19). log 2 (TPM þ 1) was calculated as the expression level, where TPM represents transcripts per million, which was defined by RNA-Seq by Expectation-Maximization (RSEM) (version 1.2.31; ref. 20). According to the results from RSEM, we removed two single-cell samples (PDX-pRCC_SC_66 and Pt-mRCC_SC_05) with total map read count less than 10,000. Finally, 8,634 PCGs and 1,874 lncRNAs with most abundant expression were selected according to two quality filtering criterions: (i) genes expressed in more than three single cells for each tumor or detected by the bulk RNA-seq and (ii) genes with average expression level more than three [log 2 (TPMþ1) > 3].

Detection of differential PCGs and lncRNAs
We applied the single-cell differential expression (SCDE) software package ("single-cell differential expression") to identify altered lncRNAs and PCGs between mRCC and pRCC single cells. To avoid the noise-induced fluctuation in gene expression due to the xenograft and cell culture, we excluded the significant differentially expressed genes (DEG) between PDX-mRCC and Pt-mRCC (21). Note that 173 lncRNAs and 1,445 PCGs were selected as the significant altered genes by the following strategies: Choose the genes with corrected Z-score of expression difference (cZ) between mRCC and pRCC more than 2 and maximum likelihood estimate (mle) more than 1, meanwhile exclude the genes whose cZ between PDX-mRCC and Pt-mRCC is more than 2 and the mle is more than 2. These differentially expressed lncRNAs were identified as CMALs.

TCGA clear cell renal cell samples classify by CMAL
The lncRNA expression and the clinical data were downloaded from TANRIC (18) and TCGA (22), respectively. Forty CMALs that do not exist in the TANRIC dataset were removed. Totally, 448 TCGA RCC samples and 134 lncRNAs were considered in this part. Unsupervised clustering was performed based on the expression level of CMALs with "Euclidean" as the distance measurement. The "cutree" function in R was used to cut the clustering tree. A small amount of patients have shown their own peculiar expression patterns. These patients will no longer be considered in cut tree, and they are referred to as "non-cluster." Then, a Kaplan-Meier survival curve and a Cox proportional hazard regression model were performed to compare survivals between the different subclass of TCGA samples.

Functional prediction of differentially expressed lncRNAs
The "guilt-by-association" analysis was used to identify coexpressed lncRNA-PCG pairs, which Pearson correlation coefficients were calculated based on the expression level between CMALs and PCGs. The threshold of Pearson correlation coefficients was set to more than 0.5. Finally, 474 and 375 PCGs were identified to correspond to up-and downregulated CMALs, respectively. All functional enrichment analyses were performed in Database for Annotation, Visualization and Integrated Discovery (23). Redundant gene ontology (GO) terms were reduced by the Revigo web tool (http://revigo.irb.hr/; ref. 24). Cytoscape v3.4.0 (25) was used to visualize the function networks.

Quantitative real-time PCR
The human renal proximal tubule epithelial cell line HK-2 and renal cancer lines 786-O and Caki-1 were kind gifts from Dr. Xinqiang Zhu laboratory (Zhejiang University School of Medicine). All experiments performed on cells that were passaged less than 25 times. Mycoplasma testing was performed during the studies, and cell cultures were free of mycoplasma. HK-2 and 786-O were cultured in RPMI-1640 with 10% FBS, and Caki-1 in MEM with 10% FBS. Reverse transcription was carried out with a PrimeScript RT reagent kit with gDNA Eraser kit (RR047A; Takara) after total cellular RNA was extracted by using Trizol reagent (15596-026; Life Technologies). Quantitative PCR was performed by using SYBR Premix Ex Taq II (RR420A; Takara) according to the manufacturer's instructions with a LightCycler system (Roche Diagnostics). Each of the samples was all prepared in triplicate. The expression of GAPDH was detected for normalization. Sequences of primers used for qRT-PCR in this study are shown in Supplementary Table S5.

Clustering single cells and identifying specific genes
A total of 10,508 abundant genes were used in the principal component analysis (PCA). We ranked all genes based on their highest absolute values of gene loading scores across the first four principal components in PCA. Using the top 1,000 genes ranked by PCA loading scores, hierarchical clustering was performed for mRCC single cells. The dendrogram was divided into four clusters. Two scores were used to measure the specificity of genes according to the principle preciously described (26). Gene specificity score of each comprising cells was computed as the OR for each gene, within a cluster compared with outside a cluster, as shown in formula: OR ij ¼ log p ij ð1Àq ij Þ q ij ð1Àp ij Þ , where, for a given threshold, p ij is the fraction of cells expressing gene i above this threshold in cluster j, and q ij is the fraction of cells expressing gene i above this threshold that are not in cluster j. Here, we used each gene's own 75 percentiles for expression thresholds. Then gene specificity score was used to rank elements within each cluster, and we chose the top 15 PCGs and lncRNAs in each cluster as specific genes. Gene expression enrichment score for each gene in each cluster is calculated by: e ij ¼ u ij =v ij , where u ij is the mean expression level of gene i within cluster j, and v ij is the mean expression level of gene i outside of cluster j.

Composition of expressed genes and lncRNAs in single cells
We reanalyzed a previously reported scRNA-seq dataset that profiled 119 single cells and 3 population samples isolated from human mRCC, PDX-mRCC, and paired PDX-pRCC tumor tissues (Supplementary Table S1; ref. 17). In total, 8,634 PCGs and 1,874 lncRNAs were detected after filtering cells and genes with low coverage, which is comparable with previous results of lncRNA detection in KIRC tissue samples (18). We applied GENCODE lncRNAs categorization to our expressed lncRNAs, resulting in the following distribution: antisense (958), lincRNA (714), sense intronic (162), and sense overlapping (40). Among them, 450 (42%) lncRNAs and 4,947 (74%) PCGs were commonly abundant expressed in all three tissues ( Supplementary Fig. S1A).
Meanwhile, we calculated the SD of lncRNA's and PCG's expression across all single cells. Although the variation was not distinguishable between four lncRNA subtypes, total lncRNAs displayed significantly higher variation than PCGs (Fig. 1A). This result indicated that the heterogeneity of lncRNA was consistently higher than PCGs whether between different tissues or cells in same tissue. Moreover, to investigate whether these lncRNAs could distinguish metastasis tumor tissue from primary tumor tissue, PCA was performed on lncRNAs ( Fig. 1B) and PCGs ( Supplementary Fig. S1B). PCA results showed that mRCC and pRCC cells represented distinct populations. In addition, the cluster of PDX-mRCC cells overlapped with which of Pt-mRCC cells, suggesting that PDX-mRCC cells can approximately represent parental mRCC.

Characterization of lncRNA expression profiles in single cells
We compared the expressed lncRNAs and PCGs in each single cell. Genes with expression levels greater than zero were regarded as expressed genes. On average, 26.6% lncRNAs and 75.1% PCGs were expressed in single cells. We also found that the percentages of lncRNAs are positively correlated with PCGs across single cells and bulk tissues, and there are more expressed genes in three bulk tissue than single cells (Supplementary Fig. S2A). Consistent with the observation, we found that lncRNAs and PCGs that were lower abundant in bulk tissues were generally detected in fewer single cells (Fig. 1C). This result suggests that gene expression level in bulk tissue may reflect the expressed proportion in single cells. Pool tissues Single cells Furthermore, lncRNAs were detected at levels 1.14-fold lower than PCGs on average (P < 2.2e-16, Mann-Whitney U) in our single-cell datasets, and it showed no difference between the four categories of lncRNAs (P > 0.01, Mann-Whitney U, Fig. 1D). This is consistent with former studies (26). Previous studies also suggested that lncRNAs can be expressed at high levels in individual cells (26), so we calculated the expression level ratio of lncRNA:PCG in bulk tissues, single cells, and normal kidney cortex. Our result demonstrated that the expression level of lncRNAs was much lower than that of PCGs (lncRNA:PCG < 1) for most single cells, but there was still few single cells with lncRNA:PCG > 1 (Supplementary Fig. S2B). All of these may serve as evidences that the expression of lncRNAs tends to have lower level and higher difference than PCGs, whereas a small-cell subpopulation tends to express lncRNAs at high levels.

Altered genes expression during RCC metastasis
As shown in Supplementary Fig. S3A, a hierarchical clustering of single cells was performed on a set of previously published metastatic tumor-specific genes (27,28). The mRCC and pRCC cells were segregated into distinct clusters, but Pt-mRCC (green) and PDX-mRCC (yellow) cells cannot be distinguished by clustering, indicating their high proximity. These results provided preliminary evidence on the alteration of the expression profile during RCC metastasis. Although these metastasis-specific genes were identified from different cancer types and individuals, some of them have also exhibited differential expression in our ccRCC single-cell samples. For example, TGFB2 can induce EMT in cancer progression and metastatic dissemination (29). In previous study, metastatic tumor-specific genes were found to be involved in the DNA damage response, chromatin modification, differentiation, apoptosis, and the cell cycle (28). However, metastasis-related genes for ccRCC are still lacking, especially for lncRNAs. Therefore, we aim to find more ccRCC metastasis-related genes and lncRNAs.
We identified DEGs using a recently developed Bayesian approach called SCDE (ref. 30 ). This may serve as a candidate set for finding the genes involved in the ccRCC metastasis. Further, we investigated the biological function of those differential expressed PCGs by functional enrichment analysis. We found both upregulated genes and downregulated genes were enriched for regulating cell death, apoptosis, and cell adhesion ( Fig. 2A). To explore whether there is any functional difference between upregulated PCGs and downregulated PCGs, we separately performed biological functional analysis on each PCG set. Strikingly, upregulated PCGs were significantly associated with mRNA synthesis (Fig. 2A, light pink background), whereas downregulated PCGs were associated with membrane process (Fig. 2A, light blue background). Cellcell adhesion is mediated by membrane process, and diminished intercellular adhesion is commonly believed as an important sign for the initiation of metastasis. Similarly, as a basic functional process, the increase of mRNA synthesis is always required for tumor cell invasion. Together, these evidences provide a possible interpretation for the metastasis mechanisms. The downregulation of the cell adhesion-related genes would facilitate the generation of free cancer cells, whereas the upregulation of the mRNA synthesis-related genes would promote the proliferation of the cancer cells. Similar results were obtained using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Fig. 2B), and the consistency suggests the reliability of our results.
We identified 173 differentially expressed lncRNAs from metastatic ccRCC (Supplementary Fig. S4). These lncRNAs may be involved in ccRCC metastasis, so we called them CMALs. From our result, there is no differentially expressed lncRNAs on chromosome 13. According to the fold changes of expression level between mRCC and pRCC, we counted the CMALs and altered PCGs for each fold-change condition. We found that the values of lncRNAs' fold changes detected within single cells were higher than PCGs' (Supplementary Fig. S4, inner circle). Intriguingly, some cancer-related lncRNAs, such as HIF1A-AS1, MALAT1, and PVT1 (33), are not present in our differentially expressed lncRNA set. This observation may suggest that these cancer-related lncRNAs are not metastasis-related. But, MEG3 (34), previously known as kidney cancer-related lncRNA, also showed low expression in metastatic RCC cells.

Validation of CMALs in the ccRCC metastasis
In order to prove that CMALs are indeed associated with ccRCC metastasis, we derived a set of 448 TCGA RCC tumors (22) according to CMALs' expression levels. The most widely used features for predicting tumor aggressiveness are the cancer stage, grade, and degree of metastasis. Besides this clinical information, mRNA clusters and miRNA clusters obtained from the previous work (22) were also used for measurement patient status. Next, unsupervised consensus clustering was used to subdivide the patient samples into six robust molecular subpopulations, and a nonclustered group was detected (Fig. 3A). We found that CMAL-based classification was highly concordant with the TCGA clinical features by x 2 tests. Although CMAL-cluster2 and CMAL-cluster4 were classified as two small subclusters, they tend to include higher stage, higher grade, higher metastasis rate, and more frequent lymph node involvement (P < 10 -3 , x 2 test; Supplementary Table S3). This indicated that CMAL-cluster2, 4 tend to be more readily metastasize to other organ and more aggressive. CMAL-clusters were also significantly associated with previously studied cluster results in mRNA and miRNA. Especially, CMAL-cluster2 and CMAL-cluster4 were significantly enriched in the mRNA-cluster2, 3 and miRNA-cluster2 (P value < 10 -5 , x 2 test) which showed high mortality in previous study (22). As aggressive clusters, the proportion of CMAL-cluster2 and CMAL-cluster4 increases with the increase of malignancy grade (Supplementary Fig. S5).Next, we performed Kaplan-Meier survival analyses to compare survival rate of CMAL-clusters (Fig. 3B). Patients with CMAL-cluster2 and CMAL-cluster4 showed lower probability of overall survival. This is consistent with the above conclusion, indicating that the detected CMALs are significantly correlated with ccRCC metastasis and do play roles in predicting patient outcome. Interestingly, patients within the nonclustered group showed relatively the most aggressive and worst overall survival rate (P < 10 -6 ).

Activated lncRNAs during metastasis
To identify the putative roles of the activated lncRNAs in ccRCC metastasis, we performed a "guilt-by-association" analysis to assess the close correlation between PCGs and CMALs, and constructed a coexpressed network between the CMALs and all PCGs in the single cells (Fig. 4A). Subsequently, the functional analysis of these correlated PCGs indicated that downregulated CMALs are clearly functionally related to cell adhesion, immune response, and cell proliferation, which are all known to be associated with cancer metastasis (Supplementary Fig. S6; refs. 35, 36). The similar result was found in downregulated PCGs.
The result of pathway analysis showed that CMALs were activated in cell adhesion, p53, and AMPK signaling pathway that were associated with cancer development. Both up-and downregulated CMALs were activated in the TNF signaling pathway (Fig. 4B). This pathway can induce MAPK and PI3K-AKt signaling pathway, thus regulating the process of apoptosis and cell adhesion (37). Previous work has indicated that TNFa induces EMT and promotes tumorigenicity of RCC by repressing E-cadherin (38). In the EMT process, epithelial cells could break free from the primary tumor site and metastasize at distant sites by attaining a mesenchymal phenotype (29). We extracted 27 lncRNAs involved in TNF signaling pathway (Supplementary Table S4). Overall survival analysis was performed to further investigate the prognostic potential of these TNF-related CMALs. We found that low expression of MEG3 and CTA-392C11 in the metastasis ccRCC cells was correlated with poor survival. Other 7 upregulated CMALs, including CTB-78F1, RP11-66B24, CTD-3131K8, DGUOK-AS1, RP11-323J4, RP11-563N4, and CTD-2020K17, have also indicated poor survival. Expression profiles of the 9 TNFrelated CMAL were displayed in Fig. 4C. The hypoxia-inducible factor (HIF)-1 signaling pathway activation by overexpression of HIF-1a is significant pathway that promotes metastasis and EMT (39). We obtained nine HIF-1related CMALs, including five upregulated CMALs (RP11-167N5, CTB-78F1, RP11-527J8, RP11-563N4, and RP11-66B24) and four downregulated CMALs (LINC01611, CTA-392C11, MEG3, and RP11-380J14; Fig. 4C; Supplementary Table S4). In addition, except for RP11-527J8, LINC01611, and RP11-380J14, other CMALs are also TNF-related. Previous work showed that TNF1a can induce HIF-1a expression, suggesting a complex system of crosstalk between hypoxia and inflammation (40). We assume that these overlapping CMALs may play a regulatory role between these two pathways.
To validate the results of single-cell RNA-seq, we performed qRT-PCR to examine the aberrant expression of CMALs in the metastasis ccRCC cell line (CAKI-1) compared with the primary ccRCC cell line (786-O) and the normal kidney cell line (HK-2). Because of the restriction of PCR primer design and the low expression of lncRNAs, we finally measured four TNF-and HIF-1-related CMALs. In addition, we examined two PCGs (SMAD3 and SERPINE1) mentioned above. Indeed, lncRNA DGUOK-AS1 and protein gene SERPINE1 were confirmed to be highly expressed in the CAKI-1 cell line, lncRNA LINC01611, and RP11-380J14, and protein gene SMAD3 showed lower expressed in the 786-O cell line, whereas lncRNA RP11-323J4 showed no significant difference in expression (Fig. 4D).

Heterogeneous expression of lncRNAs
Cancer cells in a solid tumor behave as communities. To obtain an overview of single-cell heterogeneity of lncRNAs and PCGs, we calculated the Spearman correlation coefficients between single cells and bulk tissues based on the expression data of lncRNAs and PCGs, respectively. We found that the correlation of cells in the same tissue was higher than that between different tissues, but there were still differences between the cells within the same tissue ( Supplementary Fig. S7). We extracted the correlation value between individual tissues. Correlation coefficients derived from lncRNAs profile data (0.494) are significantly lower than those derived from PCGs profile data (0.612; Fig. 5A). From the result, we found that the intratumoral heterogeneity of lncRNAs is higher than PCGs. This finding indicated that lncRNAs contribute more to intratumoral heterogeneity than PCGs.
To investigate whether intratumoral heterogeneity exist in our ccRCC metastasis tumor samples, we performed hierarchical clustering of mRCC (PDX-mRCC and Pt-mRCC) single cells with genes with high loading scores calculated from the PCA. Our algorithm identified four clusters within mRCC single cells (Fig.  5B). To identify the signatures of these subtypes of mRCC cells, we determined the specific expressed lncRNAs and PCGs in each cluster (Fig. 5C). The set of Cluster-1-specific genes contained the greatest number of specific genes compared with other subtypes. Of the top 45 specific PCGs (15 in each cluster), 6 were not annotated in GO biological function. Several ribosomal proteins, such as RPL12A, RPL27A, and RPS25, had high expression in Cluster3, and GO biological function enrichment showed that they participate in SRP-dependent cotranslational protein targeting to membrane. HEG1, HRSP12, and ARID5B from Cluster4 function in kidney development or lung development.

Discussion
LncRNAs are emerging as key players in tumorigenesis. Plenty of previous works demonstrated that these molecules have functional roles in tumor metastasis via associating with apoptosis,   cell cycle, and cellular differentiation (9,44). In this study, we reanalyzed the public paired pRCC and mRCC single-cell RNA-seq data to profile lncRNA expressions and investigated for the differences between them. We demonstrated that the expression of lncRNAs was heterogeneous across all ccRCC cells both intertumor and intratumor. To our knowledge, this is the first study systematically analyzing the differentially expressed lncRNAs and PCGs in human ccRCC metastasis. The altered expression of specific lncRNAs between metastasis and primary tumors may work as independent biomarkers for diagnosis and prognosis (45). We identified 173 CMALs, most of which were not previously reported to be implicated in cancer invasion and metastasis. CMAL clustering with TCGA expression data was significantly associated with clinical manifestations. Notably, the nonclustered patients with unique aberrant expression patterns had the worst overall survival rate. Compared with clustered patients, the CMAL expression patterns of these nonclustered patients are more complex and diverse, making it difficult to classify and analyze. Thus, the treatment and diagnosis of these patients may require more precise personalized care.
Heterogeneous tumors are composed of differentiated cancer cells, stromal cells, cancer stem cells (CSC), and so on (46). CSC is a subset of cancer cells that possess the driving force of tumorigenesis and metastasis (47). Several molecular markers, such as CD44 (48), CXCR4 (48), and ALDH1A1 (49), have been applied to identify CSCs in renal tumor. However, in our study, the SDs of these three stem-cell markers were less than 1, suggesting four subtypes do not tend to be "stem-like." One possible explanation for this finding is the small number of cell samples in our research, and the possibility of exiting CSCs in the 80 mRCC cells is slim. Another explanation is that CSCs did not exist in our ccRCC tumor sample.
Furthermore, we examined intratumor heterogeneity properties of mRCC in the context of specific genes in each subtype. Previous studies showed that tumor heterogeneity manifested in many aspects, such as genetic mutation, plastic gene expression, and protein modification (50). Our study suggests that due to tumor heterogeneity in gene expression profile, there are different tumor cell populations addicted to different oncogene pathways. These subpopulations may be developed in the differentiation of tumor stem cells or primitive cells, and have formed a complete   system to resist the resistance of human immunology and drug therapy. Compared with single genes, common pathways with subtype-specific genes may represent more robust prognostic signatures and deserve further investigation. Our study reveals a prevalent and complex lncRNA-mediated ccRCC metastasis. Those CMALs might have oncogenic or tumorsuppressive functions, and the perturbation of the lncRNAmediated tumor metastasis might be exploited for cancer prognosis and therapy.