Research Paper Volume 13, Issue 24 pp 26046—26062

An immune-related lncRNA risk coefficient model to predict the outcomes in clear cell renal cell carcinoma

Cheng Tang1, *, , GenYi Qu1, *, , Yong Xu1, , Guang Yang1, , Jiawei Wang1, , Maolin Xiang1, ,

  • 1 Department of Urology, The Affiliated Zhuzhou Hospital XiangYa Medical College CSU, Zhuzhou 412007, China
* Equal contribution

Received: September 17, 2021       Accepted: December 8, 2021       Published: December 26, 2021      

https://doi.org/10.18632/aging.203797
How to Cite

Copyright: © 2021 Tang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Objective: Using model algorithms, we constructed an immune-related long non-coding RNAs (lncRNAs) risk coefficient model to predict outcomes for patients with clear cell renal cell carcinoma (ccRCC) to understand the infiltration of tumor immune cells and the sensitivity to immune-targeted drugs.

Methods: Open genes data were downloaded from The Cancer Genome Atlas and The Immunology Database and Analysis Portal, and immune-related lncRNAs were obtained through Pearson correlation analysis. R language software was used to obtain differentially expressed immune-related lncRNAs and immune-related lncRNA pairs. The model was constructed using least absolute shrinkage and selector operation regression analysis, and receiver operator characteristic curves were drawn. The Akaike information criterion was used to distinguish the high-risk from the low-risk group. We also conducted correlation analysis for the high- and low-risk subgroups.

Results: We identified 27 immune-related lncRNAs pairs, 16 of which were included in the model construction. After merging clinical data, the areas under the curve of 1 -year, 3-year, and 5-year survival times of ccRCC patients were 0.867, 0.832, and 0.838, respectively. Subgroup analyses were conducted according to the cut-off value. We found that the high-risk group was associated with poor outcomes. The risk score and tumor stage were independent predictors of the outcome of ccRCC. The risk model predicted specific immune cell infiltration, immune checkpoint gene expression levels, and high-risk groups more sensitive to sunitinib targeted therapy.

Conclusion: We obtained prognostic-related novel ccRCC markers and risk model that predicts the outcome of patients with ccRCC and helps identify those who can benefit from sunitinib.

Introduction

Renal cancer is among the top ten most common cancers worldwide; it accounts for 5% of new cases each year. Approximately 45,520 people were diagnosed with renal cancer in the United States in 2020 [1]. Clear cell renal cell carcinoma (ccRCC) is the most common subtype of renal cancer, accounting for about 85% of the total cases, with a male-to-female ratio of 1.7:1. Most patients are middle-aged and elderly, with an average age of 64 [2].

Surgery is the standard treatment for localized ccRCC; however, this cancer carries a high risk of metastasis, poor outcome and is insensitive to conventional chemotherapies and radiotherapies [3]. In the most recent decade, with the rapid development of molecular biological tools, drug treatment of ccRCC ranged from the initial non-specific immunotherapy to targeted therapy and then to the immune checkpoint inhibitors [4]. Among them, PD-1-based blocking therapy has become the first-line therapy for patients with advanced or platinum-intolerant ccRCC [57]. In the phase III trial CheckMate 025, nivolumab was more effective than everolimus in treating patients with advanced ccRCC who had previously received treatment (the 5-year survival rate was about 26% vs. 18%). The (Food and Drug Administration) FDA approved Nivolumab for the treatment of ccRCC [8]. Although some of these drugs can improve outcomes, their effectiveness is limited [9]. Therefore, we need to explore the mechanisms of ccRCC occurrence and development and identify new biomarkers to predict related drug sensitivity.

The tumor microenvironment of ccRCC is heterogeneous, and several studies showed that the degree of immune cell infiltration in the tumor microenvironment was related to outcomes. Studies showed that high levels of CD8+ T cell infiltration in ccRCC was associated with poor outcome. Macrophages were also an essential part of the tumor microenvironment; high degrees of infiltration of M2-type macrophages are associated with tumor invasion and poor outcome of ccRCC [5]. Therefore, it is necessary to identify immune-related tumor prognostic markers.

Long non-coding RNA (lncRNA) does not code for protein. The length is more than 200 bp. It was initially considered a by-product of RNA polymerase II transcription and did not have biological functions. Many studies showed that lncRNA has a conserved secondary structure that can interact with proteins, DNA and RNA, and participate in regulating various biological processes, especially in tumors. It plays critical regulatory roles in tumors, including chromatin modification, transcription activation and inhibition, post-transcriptional mediation, and miRNA-induced molecular interference with gene expression [68]. Many studies showed that changes in molecular biology are closely related to the occurrence and development of ccRCC. Transcriptome studies described the abnormal expression of specific long non-coding RNAs, and the occurrence and progression of ccRCC have a close relationship. The expression level of lncRNA HOTAIRM1 in ccRCC decreased and inhibited the hypoxic pathway of tumor development [10]. Another study showed that lncRNA URRCC and EGFL7/P-AKT/FOXO3 signaling was related to poor outcomes and promoted the proliferation and invasion of ccRCC [11]. Another study showed that lncARSR transported by exosomes promoted the expression of AXL and c-MET in ccRCC cells by competitively binding miR-34/miR-449, rendering ccRCC patients resistant to sunitinib [12]. The binding lncRNA-LET and miR-373-3p induced the up-regulation of DKK1 and TIMP2 levels and reduced the anti-tumor effect of lncRNA-LET-mediated by ccRCC cells [13]. Other studies demonstrated that immune-related lncRNA has vital clinical significance in predicting outcomes in patients with ccRCC and as a target for targeted therapy [1419]. Therefore, the present study aimed to construct an immune-related lncRNA risk coefficient model using a model algorithm, lncRNA pairing, and iteration to predict outcomes in patients with ccRCC, understanding the tumor immune cell infiltration and the sensitivity of targeted drugs.

Results

The illustration of summary highlight was provided in Figure 1.

Summary flowchart of this study.

Figure 1. Summary flowchart of this study.

Analysis of differential expression of immune-related lncRNAs in ccRCC

The transcriptome and immune gene-related data of ccRCC were obtained from The Cancer Genome Atlas (TCGA) database and The Immunology Database and Analysis Portal (ImmPort). The Ensembl database was used to annotate and distinguish transcriptome data. Using Pearson correlation analysis, with co-expression correlation coefficient >0.4 and P < 0.001 as the identifying criteria, 433 immune-related lncRNAs were identified. We used differential expression analysis, with |log Fold Change| >1.5 and false discovery rate (FDR) <0.05 as the identifying criteria. We obtained 90 differentially expressed immune lncRNAs, and the gene heatmap (Figure 2A) was generated using R software. Sixteen lncRNAs expressions were downregulated, and 74 lncRNAs were upregulated in ccRCC (Figure 2B and Table 1).

Heatmap and differential expression analysis of immune-related lncRNA in ccRCC. (A) Heatmap of immune-related lncRNA genes between clear cell renal cell carcinoma and normal tissues. Red indicates upregulated, and blue indicates downregulated. (B) Volcano map of immune-related lncRNA between clear cell renal cell carcinoma and normal tissues. Red dots: upregulation with significant differential expression, green dots: downregulation with significant differential expression, black dots indicate no significant difference.

Figure 2. Heatmap and differential expression analysis of immune-related lncRNA in ccRCC. (A) Heatmap of immune-related lncRNA genes between clear cell renal cell carcinoma and normal tissues. Red indicates upregulated, and blue indicates downregulated. (B) Volcano map of immune-related lncRNA between clear cell renal cell carcinoma and normal tissues. Red dots: upregulation with significant differential expression, green dots: downregulation with significant differential expression, black dots indicate no significant difference.

Table 1. Immune-related lncRNAs of ccRCC obtained after differential expression analysis.

lncRNANormal-meanTumor-meanlogFCP-valueFDR
AC104984.53.6825740.462925451−2.991863.54E-382.04E-36
AC015911.30.1004970.8217785063.03167.92E-318.18E-30
PTOV1-AS20.6428682.7714832372.1080631.53E-185.01E-18
AC003070.10.2839991.3904841242.2916312.55E-231.27E-22
AC020913.30.0740730.6179734113.0605241.27E-132.62E-13
AD001527.10.1275940.5760038832.1745214.91E-139.69E-13
AC004687.10.0901070.9867543283.4529745.85E-295.24E-28
AC007098.10.1061510.5725073982.4311752.41E-281.87E-27
AC093001.10.0111440.9659695286.4376662.11E-103.40E-10
AC103563.75.1736960.597977578−3.113031.51E-364.93E-35
HOXB-AS35.5869560.94096728−2.569856.02E-361.62E-34
AL021707.60.4579541.9203098142.0680651.75E-185.68E-18
AL590764.10.153930.7534716962.2912773.26E-313.46E-30
SMIM250.257031.5188049322.5629276.13E-338.52E-32
TRG-AS10.1102140.7128909892.693381.82E-354.08E-34
PRDM16-DT13.26860.316471521−5.38982.97E-403.99E-38
TRBV11-20.1397080.86309822.6271065.48E-191.87E-18
AL049555.13.5404250.370257526−3.257325.29E-361.52E-34
SEMA6A-AS10.1447240.7403733422.3549451.72E-271.24E-26
AC091185.10.1719090.8467357212.3002657.20E-243.72E-23
AL021707.70.123280.5766570082.2257725.35E-171.43E-16
AL513327.10.2538271.0529852682.052571.35E-247.65E-24
LINC008610.0898580.7289560243.020111.93E-312.18E-30
AC079015.10.0698130.7320730123.3904131.14E-331.77E-32
AC084876.10.0829490.5723143412.786511.30E-247.46E-24
AC127024.40.2504251.0557929942.0758763.21E-178.87E-17
AC010973.20.1707490.7750833662.1824734.60E-242.44E-23
AC023421.110.335120.109043989−6.56652.25E-343.94E-33
ARHGAP27P1-BPTFP1-KPNA2P30.3173141.3676950282.1077639.89E-234.69E-22
AL513218.10.1426180.5959384022.0630124.51E-151.02E-14
CR936218.10.2592811.0939315752.0769369.57E-224.10E-21
PCED1B-AS10.4775393.1261440122.7106947.88E-361.92E-34
FOXC2-AS10.0993790.7241794662.8653389.22E-131.79E-12
LINC008930.1688030.8009610542.2463929.65E-172.51E-16
LINC020840.1474120.8327380922.4980139.70E-276.51E-26
AC011462.40.3843461.666939172.1167242.50E-187.88E-18
AL031710.110.234241.697657813−2.591795.68E-338.17E-32
AC004921.10.137650.7042398812.3550672.66E-344.28E-33
AC008735.20.6093163.2372100572.4094897.61E-202.71E-19
MCF2L-AS12.4118780.498237863−2.275251.59E-364.93E-35
AL662844.30.1150570.9892064623.1039222.20E-281.74E-27
IGFL2-AS10.0214510.8274399575.2695342.01E-072.79E-07
AC023669.20.2776181.9037757962.7776910.0007370.000869
AL135999.10.2112460.8930686722.0798477.64E-172.00E-16
AP000757.18.8711181.183177983−2.906451.56E-364.93E-35
AL022322.10.1963931.376384692.8090668.38E-244.28E-23
AC104564.30.1773130.7809168282.1388727.16E-182.15E-17
AC004253.10.1480810.8312917332.4889645.64E-222.47E-21
USP30-AS10.3942992.0726222372.3940941.15E-289.68E-28
LINC013550.2460171.1261336732.194551.08E-204.17E-20
AC004923.40.1578510.6490310362.0397243.03E-221.41E-21
AC019197.12.1445480.331330829−2.694332.06E-354.31E-34
FAM13A-AS10.2165051.2417919462.5199491.42E-291.33E-28
AC093788.10.1383160.5849699892.0803993.62E-191.25E-18
HCG270.2022451.5784507132.9643313.72E-324.68E-31
AC016773.20.1078340.6779468272.6523651.23E-204.70E-20
AC092119.20.1532250.659499182.1057178.33E-151.86E-14
AC148477.45.9740090.803677607−2.894012.37E-343.98E-33
AC005785.10.1210530.6104570232.3342519.01E-297.89E-28
AC027796.40.2242491.0666274412.249881.74E-174.93E-17
AC245884.80.2221971.3057499852.5549661.03E-214.38E-21
AC005104.10.2909041.2309542112.0811621.63E-195.75E-19
AC009549.10.2738161.4045156162.3587921.63E-311.99E-30
AP006621.20.4413351.8541613752.070823.74E-148.01E-14
AC100830.20.2363980.9966194942.0758251.22E-173.56E-17
LINC016123.1505610.338202929−3.219652.61E-461.05E-43
AL353152.12.9912480.17047375−4.133132.50E-391.75E-37
AL035661.129.664282.555782522−3.536891.54E-364.93E-35
AC012615.60.1805080.822616012.1881537.62E-233.65E-22
MALAT115.8350171.623804482.1773223.77E-105.96E-10
AL031123.110.588340.906174056−3.546541.09E-391.09E-37
AC135050.30.1606021.4822466493.2062231.68E-301.65E-29
AC006435.20.2013040.9377882142.2198842.41E-155.55E-15
LINC010940.277251.9636312572.8242626.62E-372.96E-35
AC008870.20.1404590.608695822.115575.25E-222.33E-21
AC009704.20.1502050.9407645752.6468971.69E-206.36E-20
AC243960.10.1907920.905843612.247262.19E-281.74E-27
AC008105.30.0517660.8257179213.9955774.20E-336.28E-32
AC130469.10.0574920.5849189943.3468042.14E-197.51E-19
AC116914.20.236941.3745030752.5363151.69E-227.93E-22
AL031714.10.1885080.7984534272.0825841.83E-241.03E-23
AC020907.40.1343781.0945451793.0259671.54E-251.00E-24
MMP25-AS10.3153791.6807449652.413948.09E-361.92E-34
AC063965.20.1110370.7058994322.6684232.30E-176.49E-17
AL157392.40.1340440.5661377112.0784431.63E-153.77E-15
ITGB2-AS10.2230361.6652707362.9004072.01E-291.84E-28
AC015819.20.0906530.6453245432.8316068.37E-265.53E-25
LINC003420.5915093.1522182352.4138966.71E-274.58E-26
AC004585.10.0898460.9221310893.3594422.00E-312.18E-30
MIR200CHG2.7605870.476473591−2.534511.27E-412.56E-39
Abbreviations: logFC: log fold change; FDR: false discovery rate.

Establishment of differentially expressed immune lncRNA pairs and risk coefficient scoring model

By matching the 90 differentially expressed immune lncRNA pairs for multiple cycles, a total of 2663 differentially expressed immune lncRNA pairs were obtained. Next, 27 immune lncRNA pairs were identified using least absolute shrinkage and selector operation (LASSO) regression analysis and Cox univariate regression analysis (Figure 3A). Then, a Cox multivariate regression analysis was performed based on these 27 immune lncRNA pairs, 16 of which can participate in constructing the risk coefficient scoring model (Figure 3B) and the risk coefficient of each immune lncRNA pair was obtained (Table 2).

Cox regression analysis was performed on 27 immune lncRNA pairs related to clear cell renal cell carcinoma outcome. (A) Cox univariate regression analysis forest plot of 27 immune lncRNA pairs related to the outcome of clear cell renal cell carcinoma. (B) Cox multivariate regression analysis forest plot of 27 immune lncRNA pairs related to clear cell renal cell carcinoma outcome. Red indicates risk factors, and green indicates protective factors.

Figure 3. Cox regression analysis was performed on 27 immune lncRNA pairs related to clear cell renal cell carcinoma outcome. (A) Cox univariate regression analysis forest plot of 27 immune lncRNA pairs related to the outcome of clear cell renal cell carcinoma. (B) Cox multivariate regression analysis forest plot of 27 immune lncRNA pairs related to clear cell renal cell carcinoma outcome. Red indicates risk factors, and green indicates protective factors.

Table 2. Analysis of regression coefficients of 27 pairs of immune-related lncRNA to Cox related to outcome.

lncRNA pairsCoefficientHRHR.95LHR.95HP-value
AC003070.1|LINC01355−0.32880.7197850.5009281.0342610.075434
AC007098.1|AL513218.1−0.308490.7345550.4922591.096110.13089
AC007098.1|AC093788.1−0.526590.5906160.3916290.8907080.012002
AC093001.1|MCF2L-AS10.6771761.9683121.4041972.7590528.49E-05
AC103563.7|AL031123.10.3390061.4035521.0210371.9293710.03678
HOXB-AS3|AC027796.4−0.478820.6195110.3824771.0034440.051656
SMIM25|AC008105.3−0.472750.6232840.400960.9688810.035692
SEMA6A-AS1|AC084876.10.4061831.5010780.9873182.2821780.0574
SEMA6A-AS1|CR936218.1−0.417780.6585090.4250951.0200890.06136
SEMA6A-AS1|AC104564.3−0.380170.6837440.4644031.0066840.054078
LINC00861|AC084876.1−0.314410.7302170.4889261.0905880.124475
AC079015.1|AC093788.1−0.423860.6545170.4515170.9487840.025254
AC084876.1|AC100830.20.5757641.7784891.1615222.7231710.008078
AC084876.1|AC009704.20.674631.9633061.3909422.7711940.000125
ARHGAP27P1-BPTFP1-KPNA2P3|AC116914.2−0.414810.6604650.4740510.9201820.014223
LINC00893|AC027796.4−0.36580.6936420.4507631.067390.096231
AC011462.4|MMP25-AS10.3205171.377840.9757171.9456890.068706
AL031710.1|MCF2L-AS1−0.366670.6930350.4967670.9668480.030895
AL662844.3|LINC010940.535781.708781.0555442.7662790.029265
AL662844.3|ITGB2-AS1−1.046230.3512580.217930.5661571.74E-05
AC023669.2|AC063965.2−0.317320.7280960.5160441.0272850.070809
AL022322.1|AC020907.4−0.301210.7399230.5311661.0307240.074908
AC005785.1|AC063965.20.5817871.7892331.1940422.6811070.004812
AC005785.1|AC004585.10.6569071.9288171.3090242.8420670.000895
AC005104.1|AL031714.10.5676011.764031.0812792.8778910.023033
AC100830.2|AC006435.21.1674583.2138142.0756414.9761021.66E-07
AC012615.6|AC008870.2−0.568780.5662170.3839010.8351160.004121
Abbreviations: HR: hazard ratio; HR.95L: 95% CI lower limit; HR.95H: 95% CI upper limit.

Evaluation of the prognostic predictive power of the risk model

Above 27 prognostic-related immune lncRNA pairs were used to construct the 1-year, 3-year, and 5-year receiver operator characteristic (ROC) curves of patients (Figure 4A), and the 1-year area under the curve (AUC) was calculated to be the largest AUC of 0.867 (Figure 4B). In addition, the 3-year and 5-year AUC obtained were 0.832 and 0.838, respectively, which also had predictive power. Through the best fit, the cut-off value for distinguishing between high-and low-risk groups of ccRCC patients was 2.822. We included 190 patients in the low-risk group and 360 patients in the high-risk group.

The ROC curves by establishing the risk coefficient model through the immune lncRNA pairs of ccRCC. (A) The 1, 3, and 5-year ROC curves were obtained using model construction. The AUC values were all higher than 0.83. (B) One-year ROC curve with maximum AUC value obtained by the model. (C) The cut-off value of 2.822 that distinguishes between high- and low-risk patients was obtained using the best fit.

Figure 4. The ROC curves by establishing the risk coefficient model through the immune lncRNA pairs of ccRCC. (A) The 1, 3, and 5-year ROC curves were obtained using model construction. The AUC values were all higher than 0.83. (B) One-year ROC curve with maximum AUC value obtained by the model. (C) The cut-off value of 2.822 that distinguishes between high- and low-risk patients was obtained using the best fit.

Analysis of the correlation of clinical indicators by the risk model

The relationship between the risk factor score and the risk subgroup patients ccRCC was analyzed via R language software (Figure 5A). According to the time process, the relationship between the patient’s survival status and the risk coefficient score were obtained (Figure 5B), and the Kaplan-Meier curve was constructed according to the survival status of the high- and low-risk groups (Figure 5C). The result was that the survival rate of patients in the low-risk group was significantly higher than that in the high-risk group (P < 0.001).

The risk coefficient model of ccRCC predicted outcome. (A) The risk score was divided into high- and low-risk groups. (B) Scatter plot of risk score and outcome for each patient. (C) A Kaplan-Meier curve was constructed based on the survival status of the high- and low-risk groups.

Figure 5. The risk coefficient model of ccRCC predicted outcome. (A) The risk score was divided into high- and low-risk groups. (B) Scatter plot of risk score and outcome for each patient. (C) A Kaplan-Meier curve was constructed based on the survival status of the high- and low-risk groups.

The heatmap in Figure 6A described the relationship between the level of risk scores and clinically relevant indicators. We found that the survival status of patients with ccRCC (P < 0.001), tumor grade (P < 0.001), tumor stage (P < 0.001), T stage (P < 0.001), M stage (P < 0.001), and N stage (P < 0.01) were related to risk coefficient score significantly. It can be seen from the box plot we constructed that ccRCC patients with a higher risk factor had a higher chance of death (Figure 6B). Furthermore, tumor grade (Figure 6C), tumor clinical stage (Figure 6D), T stage (Figure 6E), N stage (Figure 6F), and M stage (Figure 6G) were also higher.

ccRCC risk coefficient model for clinical correlation analysis. The clinical correlation heatmap (A) illustrating that survival (B), tumor grade (C), tumor clinical stage (D), T stage (E), N stage (F), and M stage (G) were closely related to risk factor scores.

Figure 6. ccRCC risk coefficient model for clinical correlation analysis. The clinical correlation heatmap (A) illustrating that survival (B), tumor grade (C), tumor clinical stage (D), T stage (E), N stage (F), and M stage (G) were closely related to risk factor scores.

A Cox univariate and a multivariate regression analysis were performed on the risk score and clinical correlation indicators. Then the R language software’s survival package was used to visualize the data, and forest maps were done (Figure 7A and Figure 7B). It was found that the tumor grade, tumor stage, TNM stage, and risk coefficient score were related to the outcome of the Cox univariate analysis, but in the Cox multivariate analysis, age, gender, and risk coefficient score were independent predictors of outcome. The ROC curve of clinical-related indicators and the 1-year risk coefficient score were compared in the same chart (Figure 7C). The result was that the patient’s risk coefficient score (AUC = 0.867) and tumor stage (AUC = 0.868) had the highest predictive efficacy.

Cox regression analysis of clinical correlation indicators and integrated ROC curves. (A) Clinical-related indicators Cox univariate regression analysis showing that tumor grade, clinical stage, TMN stage, and risk score were related to outcome. (B) Cox multivariate analysis showing that risk scores are independent predictors of outcome. (C) The comparison of risk coefficient score and clinical-related indicators showing that risk coefficient score (AUC = 0.867) and tumor stage (AUC = 0.868) had the highest predictive efficacy.

Figure 7. Cox regression analysis of clinical correlation indicators and integrated ROC curves. (A) Clinical-related indicators Cox univariate regression analysis showing that tumor grade, clinical stage, TMN stage, and risk score were related to outcome. (B) Cox multivariate analysis showing that risk scores are independent predictors of outcome. (C) The comparison of risk coefficient score and clinical-related indicators showing that risk coefficient score (AUC = 0.867) and tumor stage (AUC = 0.868) had the highest predictive efficacy.

Correlation analysis between risk coefficient model and immune cell infiltration

XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT were used to estimate the proportion of immune cells in these samples of ccRCC patients based on marker gene and deconvolution algorithm. Pearson correlation test was used to analyze the correlation between the risk coefficient model and tumor immune infiltrating cells with screening criteria P < 0.05, and R language software was used for data visualization (Figure 8). We found that the samples of the high-risk group were positively correlated with the infiltration of NK cells, regulatory T cells, and M1 macrophages in ccRCC and negatively correlated with the infiltration of neutrophils in ccRCC.

Correlation analysis of immune infiltrating cells in ccRCC.

Figure 8. Correlation analysis of immune infiltrating cells in ccRCC.

Correlation analysis of risk coefficient model with genes

Immune-targeted therapy is one of the most popular drugs for the treatment of renal clear cell carcinoma. We further explored the relationship between the risk coefficient model and genes, and found that among the high-risk patients, the expression levels of CTLA4 (P < 0.001; Figure 9A), LAG3 (P < 0.001; Figure 9B), PDCD1 (P < 0.001; Figure 9C), GAL9 (P < 0.001; Figure 9D), and TIGIT (P < 0.001; Figure 9E) increased. The expression level of PDCD1LG2 increased but not significantly (P > 0.05; Figure 9F). The expression levels of CD274 (P < 0.01; Figure 9G), HAVCR2 (P < 0.01; Figure 9H) decreased in this model. These genes are the potential therapeutic targets for ccRCC.

Correlation analysis of genes in patients with ccRCC. In the high-risk group, the expression levels of CTLA4 (A), LAG3 (B), PDCD1 (C), GAL9 (D), and TIGIT (E) increased. Although the expression level of PDCD1LG2 (F) increased, it was not statistically significant. CD274 (G) The expression level of HAVCR2 (H) decreased.

Figure 9. Correlation analysis of genes in patients with ccRCC. In the high-risk group, the expression levels of CTLA4 (A), LAG3 (B), PDCD1 (C), GAL9 (D), and TIGIT (E) increased. Although the expression level of PDCD1LG2 (F) increased, it was not statistically significant. CD274 (G) The expression level of HAVCR2 (H) decreased.

Correlation analysis between risk coefficient model and targeted therapy drugs

Targeted therapy drugs are the first-line therapy for patients with advanced ccRCC. We also analyzed the relationship between the risk coefficient scoring model and the sensitivity of targeted therapy drugs. The IC50 was used to evaluate the efficacy of drugs. Lower IC50 suggests higher sensitivity. We found that the high-risk group was associated with higher sensitivity of sunitinib (Figure 10A), which was statistically significant (P = 3 e-08); axitinib (Figure 10B), bevacizumab (Figure 10C), pazopanib (Figure 10D), and sorafenib (Figure 10E) were not significantly different in the high- and low-risk groups.

Correlation analysis of immune-targeted drugs in patients with ccRCC. The risk factor score was used as a potential predictor. Compared with the low-risk group, the sunitinib IC50 value of high-risk patients was lower (A), which was statistically significant (P = 3e-08). Axitinib (B), bevacizumab (C), and pazopanib (D), and sorafenib (E) were not significantly different between the high- and low-risk groups.

Figure 10. Correlation analysis of immune-targeted drugs in patients with ccRCC. The risk factor score was used as a potential predictor. Compared with the low-risk group, the sunitinib IC50 value of high-risk patients was lower (A), which was statistically significant (P = 3e-08). Axitinib (B), bevacizumab (C), and pazopanib (D), and sorafenib (E) were not significantly different between the high- and low-risk groups.

Discussion

Many studies found that transcriptome RNA expression levels are related to the outcomes of malignant tumors [2022]. Recent research on the role of non-coding RNA in ccRCC is also a focus of research [10, 11, 13, 23]. In previous studies, lncRNA-related models of ccRCC were constructed based on the expression level of transcriptome data [14, 16, 24]. In the present study, we used immune-related lncRNAs pairs to construct risk coefficient models to assess the outcome of patients with ccRCC, not based on expression levels of lncRNA. We first used TCGA and ImmPort to obtain the lncRNA and immune-related gene data of patients with ccRCC and then used R software to identify immune-related lncRNAs. Then, the differential expression of ccRCC and normal samples adjacent to the cancer were analyzed, and immune-related lncRNA pairs were obtained. We obtained the risk coefficient of each sample of ccRCC patients and established a risk coefficient model using Cox univariate factors, multivariate regression analysis, and LASSO regression analysis. By generating ROC curves, we found that the AUC for one-year outcome was the largest, and the Akaike information criterion (AIC) optimal fitting was used to obtain the critical value for distinguishing high from low-risk groups. The survival analysis of the high- and low-risk groups showed that the survival rate of patients in the low-risk group was significantly higher than that of the high-risk group (P < 0.001).

We also calculated the correlation between the risk coefficient score of each ccRCC sample and various clinical indicators. We found that age, gender, and risk coefficient score were independent predictors of outcome through Cox multivariate regression analysis. We also constructed the ROC curve of clinical-related indicators, which compared the ROC curve of the 1-year risk coefficient score in the same chart. We found that the one-year outcome risk coefficient score and tumor stage were the best predictors of ccRCC outcome, suggesting the reliability of the risk coefficient model.

To analyze the relationship between risk factor score and immune cell infiltration, the immune cell infiltration data of patients with ccRCC, we used XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT and correlation analysis. We found that the level of infiltration of NK cells, regulatory T cells, and M1 macrophages in the high-risk group was high.

Sierra et al. found that tumor-infiltrating PD-L1+ NK cells were highly expressed in renal clear cell carcinoma patients. In in vitro experiments, NK cells inhibited the proliferation of CD8+ T cells, suggesting that NK tumors infiltrating cells weakened immune regulatory functions [25]. A study showed that high levels of infiltration of dendritic cell quiescence, dendritic cell activation, mast cell quiescence, mast cell activation, and eosinophils were associated with a good outcome in ccRCC patients, while B cell memory, T cell follicular helper cells, and T cell regulation associated with poor outcome in ccRCC [26]. Xu et al. found that HK3 promoted the infiltration of monocytes and macrophages that present cell surface antigens and regulated the critical genes PD-1 and CTLA-4 of debilitating T cells, thereby affecting the immune escape process [27].

We also performed correlation analysis on the risk model for immune checkpoint genes and targeted therapy drugs and found that expression levels of CTLA4, LAG3, PDCD1, GAL9, and TIGIT increased, while expression levels of CD274 and HAVCR2 decreased in samples from patients in the high-risk group; these can be used as immune targeted therapy with potential therapeutic targets.

ccRCC discovered early and mid-term can be removed surgically, while patients with advanced ccRCC experience poor outcomes because of metastasis and missing the optimal time for surgery. Immune-targeted drug therapy has become the first-line treatment for patients with advanced ccRCC. Of these, vascular endothelial growth factor monoclonal antibody and tyrosine kinase inhibitors are the primary drugs used for anti-tumor angiogenesis therapy [2830]. However, changes in the tumor microenvironment may be associated with the emergence of resistance of ccRCC to immune-targeted drugs. Therefore, identifying sensitive drugs may reduce treatment costs and reduce the side effects of immune-related drugs. In the risk coefficient model, we included sunitinib, axitinib, bevacizumab, pazopanib, and sorafenib and found that patients with ccRCC in the high-risk group were more sensitive to sunitinib than the low-risk group.

This study has some limitations, although we adopted rigorous methods and algorithms to build the model. It is necessary to validate the reliability of our risk coefficient model using external data. In a future study, we will collect more clinical data and expand the sample size.

Outcome-related novel ccRCC markers and risk coefficient model were obtained by constructing immune-related lncRNA pairs, which can predict outcomes of ccRCC and help distinguish those who can benefit from sunitinib.

Materials and Methods

Data acquirement

The GDC Data Transfer Tool was used to download open transcriptome data from TCGA (https://cancergenome.nih.gov/) [31] of ccRCC and normal tissues adjacent to the cancer, which included 539 cases of ccRCC samples and 72 cases of normal tissue samples adjacent to cancer. Then gene transfer format (GTF) files were downloaded using Ensembl (http://asia.ensembl.org) [32] to annotate and distinguish the mRNA and lncRNA of the transcriptome data. Immune-related genes were obtained from ImmPort (http://www.immport.org).

Differential expression analysis of immune-related lncRNAs

Immune-related lncRNAs were identified based on the co-expression strategy and Pearson correlation analysis according to the co-expression correlation coefficient > 0.4 and P < 0.001. DEirlncRNA was selected using the limma package of R language software [33] with |log Fold Change|>1.5 and FDR < 0.05. Obtained lncRNAs were visualized using the heatmap package.

Construction of immune-related lncRNA pairs

DEirlncRNAs were identified using multiple rounds of pairing. Parameter values of 0 or 1 were used for definition, and α was defined as the parameter value. If the expression of lncRNA A was greater than that of lncRNA B in an immune-related lncRNA pair in a sample, then the α value of the lncRNA pair was 1; otherwise, the value of α was 0. If the ratio of the α value (either 0 or 1) of the immune-related lncRNA pair in all samples is less than 80%, it means that the immune-related lncRNA pair was an effective match; otherwise, it needs to be re-paired.

Acquisition of clinical data and establishment of model

First, clinical data related to ccRCC were downloaded from TCGA. Then, the limma package of R software was used to match the immune-related lncRNA pairs in the previous step. We then took the intersection and deleted the repeated clinical data with no follow-up time. A single factor regression analysis was performed on the immune-related lncRNA pairs initially obtained, and the immune-related lncRNA pairs related to the survival status were found. The significance screening criterion was P < 0.01.

To prevent over-fitting, the glmnet package of R language was used to perform LASSO regression analysis [34] on the obtained immune lncRNA pairs, run 1000-repeated random cycles, and immune lncRNA pairs with a matching frequency of more than 100 times identified those with P < 0.05 after the second cross-validation. The best pairing combination was selected to obtain immune lncRNA pairs that can participate in constructing the Cox risk coefficient model. By constructing Cox univariate and multivariate analysis models, the risk coefficient of each immune lncRNA pair related to the outcome was obtained, and the risk score of each patient’s tumor sample was determined. The total risk score of each ccRCC patient sample was equal to the sum of the expression amount of each immune lncRNA pair in the sample multiply risk coefficient. The formula is following:

Risk Score=i=1nRisk Coefficienti×IrlncRNA Expressioni

The Cox analysis results were visualized using the survminer and survival packages of R software.

Construction of the ROC curve with risk coefficient model

ROC curves were constructed using the survivalROC package of the R software, which included ROC for 1-, 3-, and 5- years and the AUC values were calculated to determine the value predicted by the model. We found that the 1-year ROC curve had the largest AUC value. According to the AIC best fit [35], it was possible to distinguish low-risk with high-risk patients by finding a critical value with the largest sum of specificity and sensitivity.

Clinical correlation analysis with risk coefficient model

Survival and survminer packages of the R language software were used to compare the survival differences between the high- and low-risk groups. P < 0.001 indicated a significant difference. A Kaplan-Meier curve was constructed to visualize the data. The relationships between the risk score and the previously obtained clinical indicators (survival status, age, gender, tumor grade, tumor stage, and T, N, and M stages) were analyzed using the chi-square test. The relationships between the risk score and the different subgroups of clinical indicators were analyzed using the Wilcoxon rank-sum test. The limma package and ggpubr package of the R language software were used to visualize the data. To determine whether the risk score can be used as an independent predictor related to the outcome of patients with ccRCC, we performed Cox univariate and multivariate regression analysis on the risk score and clinical correlation indicators, which used the hazard ratio to evaluate. P < 0.05 was the identification criterion, and the survival package of R software was used to visualize the data. To compare the accuracy of the risk score and clinically relevant indicators in predicting survival and outcome, we compared the ROC curves obtained for the 1-year follow-up with the ROC curve of clinically relevant indicators in the same graph.

Correlation analysis of immune cells

To analyze the relationship between risk factor score and immune cell infiltration, the immune cell infiltration data of patients with ccRCC in TCGA was calculated based on CIBERSORT (http://cibersort.stanford.edu/) [36], TIMER (version 2.0; http://timer.cistrome.org/) [37], QUANTISEQ (http://icbi.at/quantiseq) [38], Microenvironment Cell Populations-counter [39], EPIC (http://epic.gfellerlab.org) [40], and XCELL (http://xCell.ucsf.edu/) [41]. The correlation between immune cell infiltration data and risk coefficient score was analyzed using the limma, scales, ggplot2, and ggtext packages of R software, which can be visualized to obtain a bubble chart according to P < 0.05 as the identifying criterion.

Gene correlation analysis

We found that CD274, CTLA4, HAVCR2, LAG3, LGALS9, PDCD1, PDCD1LG2, and TIGIT were abundantly expressed in ccRCC samples. To determine whether these genes differed between the high- and low-risk groups of the risk model, the limma and ggpubr packages of the R language software were used to analyze and visualize the data using violin charts.

Correlation analysis of targeted drugs

To determine whether there was a difference in patients’ response in the high- and low-risk groups of ccRCC patients to targeted drugs, the half-inhibition rate (IC50) of the drug was used as an index to measure drug sensitivity. The data were analyzed and visualized using the limma package, ggpubr, ggplot, and pRRophetic package in R software.

Abbreviations

ccRCC: clear cell renal cell carcinoma; TCGA: The Cancer Genome Atlas; ImmPort: The Immunology Database and Analysis Portal; lncRNAs: long noncoding RNAs; GTF: Gene Transfer Format; DEirlncRNAs: Differential Expression immune-related long noncoding RNAs; FDR: false discovery rate; LASSO: Least Absolute Shrinkage and Selector Operation; RS: risk score; ROC: Receiver Operator Characteristic; OS: overall survival; AUC: area under the curve; AIC: Akaike information criterion; MCPcounter: Microenvironment Cell Populations-counter; IC50: half inhibitory concentration; HR: hazard ratio; HR.95L: 95% CI lower limit; HR.95H: 95% CI upper limit; logFC: log fold change.

Author Contributions

C. Tang wrote the main manuscript text and analyzed the data; G. Qu performed experiments and analyzed data; Y. Xu proposed the idea; G. Yang, J. Wang, and M. Xiang collected data. All the authors reviewed the manuscript and discussed the results and edited the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

This work was supported in part by grants from Zhuzhou City Science and Technology Plan Project (#2019-001).

References

  • 1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020; 70:7–30. https://doi.org/10.3322/caac.21590 [PubMed]
  • 2. Barata PC, Rini BI. Treatment of renal cell carcinoma: Current status and future directions. CA Cancer J Clin. 2017; 67:507–24. https://doi.org/10.3322/caac.21411 [PubMed]
  • 3. Pirrotta MT, Bernardeschi P, Fiorentini G. Targeted-therapy in advanced renal cell carcinoma. Curr Med Chem. 2011; 18:1651–7. https://doi.org/10.2174/092986711795471293 [PubMed]
  • 4. Choueiri TK, Motzer RJ. Systemic Therapy for Metastatic Renal-Cell Carcinoma. N Engl J Med. 2017; 376:354–66. https://doi.org/10.1056/NEJMra1601333 [PubMed]
  • 5. Motzer RJ, Penkov K, Haanen J, Rini B, Albiges L, Campbell MT, Venugopal B, Kollmannsberger C, Negrier S, Uemura M, Lee JL, Vasiliev A, Miller WH Jr, et al. Avelumab plus Axitinib versus Sunitinib for Advanced Renal-Cell Carcinoma. N Engl J Med. 2019; 380:1103–15. https://doi.org/10.1056/NEJMoa1816047 [PubMed]
  • 6. St Laurent G, Wahlestedt C, Kapranov P. The Landscape of long noncoding RNA classification. Trends Genet. 2015; 31:239–51. https://doi.org/10.1016/j.tig.2015.03.007 [PubMed]
  • 7. Goodall GJ, Wickramasinghe VO. RNA in cancer. Nat Rev Cancer. 2021; 21:22–36. https://doi.org/10.1038/s41568-020-00306-0 [PubMed]
  • 8. Statello L, Guo CJ, Chen LL, Huarte M. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol. 2021; 22:96–118. https://doi.org/10.1038/s41580-020-00315-9 [PubMed]
  • 9. Zahoor H, Barata PC, Jia X, Martin A, Allman KD, Wood LS, Gilligan TD, Grivas P, Ornstein MC, Garcia JA, Rini BI. Patterns, predictors and subsequent outcomes of disease progression in metastatic renal cell carcinoma patients treated with nivolumab. J Immunother Cancer. 2018; 6:107. https://doi.org/10.1186/s40425-018-0425-8 [PubMed]
  • 10. Hamilton MJ, Young M, Jang K, Sauer S, Neang VE, King AT, Girke T, Martinez E. HOTAIRM1 lncRNA is downregulated in clear cell renal cell carcinoma and inhibits the hypoxia pathway. Cancer Lett. 2020; 472:50–8. https://doi.org/10.1016/j.canlet.2019.12.022 [PubMed]
  • 11. Zhai W, Zhu R, Ma J, Gong D, Zhang H, Zhang J, Chen Y, Huang Y, Zheng J, Xue W. A positive feed-forward loop between LncRNA-URRCC and EGFL7/P-AKT/FOXO3 signaling promotes proliferation and metastasis of clear cell renal cell carcinoma. Mol Cancer. 2019; 18:81. https://doi.org/10.1186/s12943-019-0998-y [PubMed]
  • 12. Qu L, Ding J, Chen C, Wu ZJ, Liu B, Gao Y, Chen W, Liu F, Sun W, Li XF, Wang X, Wang Y, Xu ZY, et al. Exosome-Transmitted lncARSR Promotes Sunitinib Resistance in Renal Cancer by Acting as a Competing Endogenous RNA. Cancer Cell. 2016; 29:653–68. https://doi.org/10.1016/j.ccell.2016.03.004 [PubMed]
  • 13. Ye Z, Duan J, Wang L, Ji Y, Qiao B. LncRNA-LET inhibits cell growth of clear cell renal cell carcinoma by regulating miR-373-3p. Cancer Cell Int. 2019; 19:311. https://doi.org/10.1186/s12935-019-1008-6 [PubMed]
  • 14. Sun Z, Jing C, Xiao C, Li T. Long Non-Coding RNA Profile Study Identifies an Immune-Related lncRNA Prognostic Signature for Kidney Renal Clear Cell Carcinoma. Front Oncol. 2020; 10:1430. https://doi.org/10.3389/fonc.2020.01430 [PubMed]
  • 15. Khadirnaikar S, Kumar P, Pandi SN, Malik R, Dhanasekaran SM, Shukla SK. Immune associated LncRNAs identify novel prognostic subtypes of renal clear cell carcinoma. Mol Carcinog. 2019; 58:544–53. https://doi.org/10.1002/mc.22949 [PubMed]
  • 16. Jiang Y, Gou X, Wei Z, Tan J, Yu H, Zhou X, Li X. Bioinformatics profiling integrating a three immune-related long non-coding RNA signature as a prognostic model for clear cell renal cell carcinoma. Cancer Cell Int. 2020; 20:166. https://doi.org/10.1186/s12935-020-01242-7 [PubMed]
  • 17. Tian P, Wei JX, Li J, Ren JK, Yang JJ. LncRNA SNHG1 regulates immune escape of renal cell carcinoma by targeting miR-129-3p to activate STAT3 and PD-L1. Cell Biol Int. 2021; 45:1546–60. https://doi.org/10.1002/cbin.11595 [PubMed]
  • 18. Su Y, Zhang T, Tang J, Zhang L, Fan S, Zhou J, Liang C. Construction of Competitive Endogenous RNA Network and Verification of 3-Key LncRNA Signature Associated With Distant Metastasis and Poor Prognosis in Patients With Clear Cell Renal Cell Carcinoma. Front Oncol. 2021; 11:640150. https://doi.org/10.3389/fonc.2021.640150 [PubMed]
  • 19. Li CS, Lu ZZ, Fang DL, Zhou WJ, Wei J. Immune-related long non-coding RNAs can serve as prognostic biomarkers for clear cell renal cell carcinoma. Transl Androl Urol. 2021; 10:2478–92. https://doi.org/10.21037/tau-21-445 [PubMed]
  • 20. Shen X, Piao L, Zhang S, Cui Y, Cui Y, Quan X, Sun H. Long non-coding RNA activated by TGF-β expression in cancer prognosis: A meta-analysis. Int J Surg. 2018; 58:37–45. https://doi.org/10.1016/j.ijsu.2018.08.004 [PubMed]
  • 21. El-Ashmawy NE, Al-Ashmawy GM, Hamouda SM. Long non-coding RNA FAM83H-AS1 as an emerging marker for diagnosis, prognosis and therapeutic targeting of cancer. Cell Biochem Funct. 2021; 39:350–6. https://doi.org/10.1002/cbf.3601 [PubMed]
  • 22. Yu M, Lu B, Zhang J, Ding J, Liu P, Lu Y. tRNA-derived RNA fragments in cancer: current status and future perspectives. J Hematol Oncol. 2020; 13:121. https://doi.org/10.1186/s13045-020-00955-6 [PubMed]
  • 23. Zhong W, Zhang F, Huang C, Lin Y, Huang J. Identification of Epithelial-Mesenchymal Transition-Related lncRNA With Prognosis and Molecular Subtypes in Clear Cell Renal Cell Carcinoma. Front Oncol. 2020; 10:591254. https://doi.org/10.3389/fonc.2020.591254 [PubMed]
  • 24. Zhang H, Qin C, Liu HW, Guo X, Gan H. An Effective Hypoxia-Related Long Non-Coding RNAs Assessment Model for Prognosis of Clear Cell Renal Carcinoma. Front Oncol. 2021; 11:616722. https://doi.org/10.3389/fonc.2021.616722 [PubMed]
  • 25. Sierra JM, Secchiari F, Nuñez SY, Iraolagoitia XLR, Ziblat A, Friedrich AD, Regge MV, Santilli MC, Torres NI, Gantov M, Trotta A, Ameri C, Vitagliano G, et al. Tumor-Experienced Human NK Cells Express High Levels of PD-L1 and Inhibit CD8+ T Cell Proliferation. Front Immunol. 2021; 12:745939. https://doi.org/10.3389/fimmu.2021.745939 [PubMed]
  • 26. Pan Q, Wang L, Chai S, Zhang H, Li B. The immune infiltration in clear cell Renal Cell Carcinoma and their clinical implications: A study based on TCGA and GEO databases. J Cancer. 2020; 11:3207–15. https://doi.org/10.7150/jca.37285 [PubMed]
  • 27. Xu W, Liu WR, Xu Y, Tian X, Anwaier A, Su JQ, Zhu WK, Shi GH, Wei GM, Huang YP, Qu YY, Zhang HL, Ye DW. Hexokinase 3 dysfunction promotes tumorigenesis and immune escape by upregulating monocyte/macrophage infiltration into the clear cell renal cell carcinoma microenvironment. Int J Biol Sci. 2021; 17:2205–22. https://doi.org/10.7150/ijbs.58295 [PubMed]
  • 28. Powles T, Albiges L, Staehler M, Bensalah K, Dabestani S, Giles RH, Hofmann F, Hora M, Kuczyk MA, Lam TB, Marconi L, Merseburger AS, Fernández-Pello S, et al. Updated European Association of Urology Guidelines: Recommendations for the Treatment of First-line Metastatic Clear Cell Renal Cancer. Eur Urol. 2018; 73:311–5. https://doi.org/10.1016/j.eururo.2017.11.016 [PubMed]
  • 29. Bedke J, Albiges L, Capitanio U, Giles RH, Hora M, Lam TB, Ljungberg B, Marconi L, Klatte T, Volpe A, Abu-Ghanem Y, Dabestani S, Fernández-Pello S, et al. Updated European Association of Urology Guidelines on Renal Cell Carcinoma: Nivolumab plus Cabozantinib Joins Immune Checkpoint Inhibition Combination Therapies for Treatment-naïve Metastatic Clear-Cell Renal Cell Carcinoma. Eur Urol. 2021; 79:339–42. https://doi.org/10.1016/j.eururo.2020.12.005 [PubMed]
  • 30. Lai Y, Zhao Z, Zeng T, Liang X, Chen D, Duan X, Zeng G, Wu W. Crosstalk between VEGFR and other receptor tyrosine kinases for TKI therapy of metastatic renal cell carcinoma. Cancer Cell Int. 2018; 18:31. https://doi.org/10.1186/s12935-018-0530-2 [PubMed]
  • 31. Tomczak K, Czerwińska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Pozn). 2015; 19:A68–77. https://doi.org/10.5114/wo.2014.47136 [PubMed]
  • 32. Yates AD, Achuthan P, Akanni W, Allen J, Allen J, Alvarez-Jarreta J, Amode MR, Armean IM, Azov AG, Bennett R, Bhai J, Billis K, Boddu S, et al. Ensembl 2020. Nucleic Acids Res. 2020; 48:D682–8. https://doi.org/10.1093/nar/gkz966 [PubMed]
  • 33. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
  • 34. Engebretsen S, Bohlin J. Statistical predictions with glmnet. Clin Epigenetics. 2019; 11:123. https://doi.org/10.1186/s13148-019-0730-1 [PubMed]
  • 35. Bozdogan H. Akaike's Information Criterion and Recent Developments in Information Complexity. J Math Psychol. 2000; 44:62–91. https://doi.org/10.1006/jmps.1999.1277 [PubMed]
  • 36. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–7. https://doi.org/10.1038/nmeth.3337 [PubMed]
  • 37. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017; 77:e108–10. https://doi.org/10.1158/0008-5472.CAN-17-0307 [PubMed]
  • 38. Plattner C, Finotello F, Rieder D. Deconvoluting tumor-infiltrating immune cells from RNA-seq data using quanTIseq. Methods Enzymol. 2020; 636:261–85. https://doi.org/10.1016/bs.mie.2019.05.056 [PubMed]
  • 39. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, de Reyniès A. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016; 17:218. https://doi.org/10.1186/s13059-016-1070-5 [PubMed]
  • 40. Racle J, Gfeller D. EPIC: A Tool to Estimate the Proportions of Different Cell Types from Bulk Gene Expression Data. Methods Mol Biol. 2020; 2120:233–48. https://doi.org/10.1007/978-1-0716-0327-7_17 [PubMed]
  • 41. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017; 18:220. https://doi.org/10.1186/s13059-017-1349-1 [PubMed]