Abstract
Purpose: This study aims to construct the stromal immunotype, which could improve the prediction of postsurgical survival and adjuvant platinum-based chemotherapy in muscle-invasive bladder cancer (MIBC).
Experimental Design: A total of 118 patients with MIBC from Shanghai Cancer Center, 140 patients with MIBC from Zhongshan Hospital, and 287 patients with MIBC from TCGA cohort were included in the study. Immune cell infiltration was evaluated by IHC staining or CIBERSORT method. Five immune features were selected out of 22 immune features to construct immunotypes based on the LASSO Cox regression model.
Results: Using the LASSO model, we classified patients with MIBC into stromal immunotype A subgroup (CTLhighNKhighTreglowMacrophagelowMClow) and stromal immunotype B subgroup (CTLlowNKlowTreghighMacrophagehighMChigh). Significant differences were found between immunotype A and immunotype B in the combined cohort with 5-year overall survival (OS, 76.0% vs. 44.0%; P < 0.001) and 5-year disease-free survival (DFS, 62.8% vs. 48.3%; P < 0.001). Stromal immunotype was revealed to be an independent prognostic indicator in multivariate analysis in all cohorts separately. Either OS or DFS was not improved by adjuvant chemotherapy (ACT) in pT2 stage patients or pT3+pT4 patients, but further analysis revealed that OS and disease-free was significantly improved by ACT in pT3+pT4 patients (P = 0.016 and P = 0.006, respectively). Finally, stromal immunotype A showed higher immune checkpoint molecules (PD-L1, PD-1, and CTLA-4) expression.
Conclusions: The stromal immunotypes could effectively predict survival and recurrence of MIBC. Furthermore, the immunotypes might be a practical predictive tool to identify pT3+pT4 patients who would benefit from ACT. Clin Cancer Res; 24(13); 3069–78. ©2018 AACR.
Both malignant cells and stromal cells were known for orchestrating tumor progression in muscle-invasive bladder cancer (MIBC). However, despite the increasing studies of the clinical significance of immune elements in MIBC, no study has evaluated the whole immune cells' landscape based on their location. Here, we thoroughly evaluate immune cells in tumor core and stroma in two independent cohorts, and then identify two major immune phenotypes by using bioinformatics and IHC methods. Immunotype A tumors, characterized by relatively higher levels of CD8+ T cells and CD56+ NK cells, are associated with better prognosis in multivariate analysis and higher immune checkpoint (PD-L1, PD-1, and CTLA-4) mRNA expression. Immunotype B tumors, which contain higher levels of FoxP3+ T regulatory cells, macrophages, and mast cells, derive more benefit from adjuvant chemotherapy in higher state (T3 and T4) patients. Therefore, the immunotype could be a potential prognostic indicator and predictor for adjuvant chemotherapy and immunotherapy.
Introduction
As the fifth most common cancer, bladder cancer is diagnosed in more than 430,000 patients worldwide each year (1). Although representing only about 25% of patients, patients with muscle-invasive bladder cancer (MIBC) have a much worse prognosis than non-muscle–invasive disease and often need radical cystectomy. Besides surgical treatment, platinum-based chemotherapy is considered as a first-line regimen for advanced or metastatic urothelial bladder cancer (UBC; ref. 2). However, although neoadjuvant chemotherapy remains the preferred approach for patients with MIBC, the role of adjuvant chemotherapy still has not been established in fully accrued (3). A statistically significant progression-free survival improvement was observed in the EORTC intergroup trial, whereas an Italian trial accruing over 180 patients had reported conflicting results (4, 5). These conflicting data suggested that besides traditional clinical high-risk feature (pT3+pT4 or N+), more specific markers and subtype classification are required to aid in patient selection.
Recent studies showed that bladder tumors can be grouped into basal and luminal subtypes based on gene expression patterns (6, 7). However, these studies mostly focused on the characteristics of malignant cells. Both malignant cells and stromal cells were known for orchestrating tumor-associated inflammation, tumor progression, and metastasis (8). On the other hand, growing studies suggested that immune response plays a critical role in chemotherapy-induced cytotoxicity and might indicate patients' response to cytotoxic regimens (9–11). By far, although some studies have unveiled the clinical significance of some immune cells like CD3+ tumor-infiltrating lymphocytes (TILs), CD8+ cytotoxic T cells (CTLs), CD68+ tumor-associated macrophages (TAM), and foxp+ regulatory T cells (Treg) in MIBC, no study has evaluated the whole immune cells' landscape based on their location (12–17).
Therefore, in this study, we used CIBERSORT, a bioinformatics method (18), to screen 22 immune cells phenotypes and find out the most relevant prognostic immune cells. Then, we performed IHC to locate these immune cells at tumor core and stroma, which aimed to construct immunotypes to predict survival and benefits of adjuvant chemotherapy in patients with MIBC.
Materials and Methods
Study population
The study design was shown in Fig. 1A. The study was conducted in accordance with the International Ethical Guidelines for Biomedical Research Involving Human Subjects (CIOMS) and performed after approval by the Institutional Review Board (IRB) of Fudan University Shanghai Cancer Center. Three independent cohorts of patients with MIBC were included. The training cohort, including 118 patients with MIBC, underwent radical cystectomy between 2008 and 2012 in Fudan University Shanghai Cancer Center (Shanghai, China). The testing cohort was composed of 140 patients with MIBC who underwent radical cystectomy between 2002 and 2014 in Zhongshan Hospital (Shanghai, China). All the patients signed informed consent given by the investigators. Following criteria were conformed in all the cases: (i) informed consent; (ii) no comorbidities (e.g., had suffered from other malignant tumors.); (iii) complete available follow-up data; (iv) pT≥2; (v) performed cystectomy and confirmed postoperative histopathology diagnosis; and (vi) tumor stage classification was carried out or updated according to 2009 TNM classification. Patients who had metastatic disease or underwent prior chemotherapy or radiotherapy at the time of cystectomy were excluded. Five patients and seven patients were excluded because of metastatic disease in the training cohort and testing cohort, respectively. One patient in the training cohort received prior radiotherapy (XRT) in our cohort. Three patients and four patients were excluded because of the prior chemotherapy. Pathologic data consisted of tumor grade, tumor and nodal stage, and histotype according to the 2004 World Health Organization (WHO) classification. The platinum-based treatment was recorded as four cycles of chemotherapy with gemcitabine/cisplatin (GC) or methotrexate/vinblastine/adriamycin/cisplatin(M-VAC). Gemcitabine/carboplatin or gemcitabine/oxaliplatin were administered when patients did not tolerate cisplatin. Overall survival (OS) was computed as the time from the day of surgery to death from all causes. Disease-free survival (DFS) was computed from the date of surgery to the first recurrence (local recurrence or distant metastases, identified by physical examination, CT, or MRI) or the most recent follow-up. The follow-up protocol was composed of patient history, physical examination, and laboratory measurement every 3 months in the first year, every 6 months in the second year, and annually thereafter routinely. We ended follow-up in June 2016. The median follow-up in the training cohort, testing cohort, and validation cohort were 29, 56, and 20 months, respectively. A total of 47, 77, and 89 patients died in these cohorts, separately.
The validation cohort was composed of 287 patients with MIBC from TCGA provisional database. Biospecimens were collected by either transurethral resection or radical cystectomy. The selection criteria were: (i) available for clinical data and mRNA expression data; (ii) available for OS data and DFS data; (iii) available for tumor grade; and (iv) pT stage ≥ 2. Demographics of patients are presented in Fig. 1B. The Philadelphia neoadjuvant chemotherapy (NAC) cohort data were obtained from the GEO database (GSE48277; ref. 7).
IHC and evaluation of immunostaining
Tissue microarrays were constructed with two 1.0-mm tissue cores from two distinct areas, and at least one out of the two selected tissue cores was available to locate stroma. IHC was performed as described in Supplementary Data S1. The CIBERSORT method was performed to analyze associations between clinical outcomes and 22 leukocyte phenotypes in TCGA cohort. Univariate Cox regression was performed to select six relevant prognostic leukocytes as presented in Supplementary Table S1: Mast cells (tryptase), macrophages (CD68), Tregs (foxp3+), NK cells (CD56), CTLs (CD8), and follicular helper T cells (Bcl-6+ and CD3+). These immune cells were stained with different antibody as summarized in Supplementary Table S2. No follicular helper T cell was found after dual stain with anti-Bcl6 and anti-CD3 antibody. The negative controls were equally treated with the primary antibody omitted. All the slides were analyzed with Leica DM6000 B (Leica Microsystems). The number of immune cells per core was evaluated with Image Pro plus (Media Cybernetics Inc.). To evaluate the density of stained immune cells, three respective areas of stroma and tumor core were evaluated at ×200 magnification and the mean value was adopted. The counts of all cells were defined as the number of nucleated stained cells per field and change into density as cells/mm2. Two pathologists (Y. Zhu and Z. Liu) who were blinded to the clinicopathologic data and scored all samples independently, and the mean count was adopted. The intraclass correlation between the two pathologists from the same slides was 0.935 [95% confidence interval (CI), 0.916–0.942; P < 0.001]. The tumor mutation burden (TMB) was calculated from comprehensive genomic profiling data of TCGA cohort. The three major T-cell categories (i.e., immune desert, immune excluded, and inflamed) were evaluated on the basis of the density of CD8+ CTLs in tumor core and stroma (19). The represented IHC images were shown in Supplementary Fig. S1.
Statistical analysis
The LASSO Cox regression model was applied to determine the ideal coefficient for each prognostic feature and estimate the likelihood deviance. The coefficients and partial likelihood deviance were calculated with “glmnet” package in R. In TCGA cohort, the relative leukocyte fraction was evaluated instead of leukocyte density. GSEA was used to identify the pathways that were significantly enriched between stromal immunotype A and stromal immunotype B (20). Statistical analysis was performed with MedCalc software (version 18; MedCalc), Stata (version 14; StataCorp.), and R software packages, version 3.4.2 (The R foundation for Statistical Computing, https://www.r-project.org/) in the study.
Results
Construction of stromal immune index and definition of stromal immunotype
Tumor core and stromal mast cells, macrophages, Tregs, NK cells, and CTLs stained distinguishable and clear as shown in Fig. 2A in training cohort and testing cohort. All five immune cells were significantly more abundant in the stroma in training and testing cohorts (Fig. 2B). The mean densities (/mm2) of mast cells, macrophages, Tregs, NK cells, and CTLs in tumor core were 3.525, 14.550, 1.330, 0.652, and 4.432, respectively. Correspondingly, the mean densities of these cells in stroma were 11.584, 53.822, 19.161, 3.567, and 23.194, respectively. Stromal immune cells' infiltration was showed to be significantly relevant to MIBC patient prognosis compared with infiltration in tumor core (Supplementary Table S3). Given the abundant infiltration and tight outcome correlation, five stromal immune cells in training cohort were selected to construct stromal immune index (IIstromal) by the LASSO Cox regression method. With the LASSO Cox regression models, the coefficients for these five features were calculated when Logγ = −4 as shown in Fig. 2C. The partial likelihood deviance for the selected lambda was 10.243 as presented in Fig. 2D. A formula was derived on the basis of the five stromal immune cell infiltration feature, where IIstromal = 2(0.0294 × MCstromal+ 0.01151 × Macrophagestromal+ 0.00585 × Tregstromal− NKstromal× 0.00808 − 0.0295 × CTLstromal). A modified formula was applied to the TCGA cohort leukocyte fraction data to minimize bias as: modified IIstromal = 2[n × (0.0294 × MCstromal+ 0.01151 × Macrophagestromal+ 0.00585 × Tregstromal− NKstroma× 0.00808 − 0.0295 × CTLstromal)], n = 378.5 (modified factor). In these formulas, low expression status was equivalent to 0, and high expression status was equivalent to 5.
On the basis of IIstromal, stromal immunotype A was defined as stromal immune index <1 and stromal immunotype B was defined as stromal immune index ≥1. Distribution of clinical characteristics did not vary significantly between stromal immunotype A and stromal immunotype B (Supplementary Table S4).
Association of stromal immunotype with OS and DFS
To validate the rationality of new defined stromal immunotype, hierarchical tree analysis was performed to classify different subgroups of patients and clusters of immune cells in training cohort FUSCC, testing cohort, and validation cohort (Fig. 3A–C). In all cohorts, a CTLhighNKhighTreglowMacrophagelowMClow group and a CTLlowNKlowTreghighMacrophagehighMChigh group could be discriminate, whereas the latter one showed worse 5-year OS and DFS. To evaluate the prognostic value of IIstromal as a linear variable, we performed smooth HR curves of OS and DFS in three cohorts as shown in Fig. 3D–F and Supplementary Fig. S2A–S2C. The cut-off value of IIstromal (1) was used as the reference line. When IIstromal was below 1, 95% CI of Ln (HR) area was lower than the reference line and vice versa. Univariate Cox regression analysis identified IIstromal was a statistically significant clinical factor associated with OS and DFS.
To further evaluate the prognostic value of stromal immunotype, we compare OS and DFS between stromal immunotype subgroups with Kaplan–Meier survival analysis in three cohorts. Patients with stromal immunotype B showed significantly worse OS in the training, testing, and validation cohorts (P < 0.001, P < 0.001, and P < 0.001, respectively) as shown in Fig. 3G–I. Also, the DFS in three sets showed the consistent results (P = 0.008, P < 0.001, and P < 0.001) as presented in Supplementary Fig. S2D–S2F. After the combination of all three cohorts, significant differences were found between immunotype A and immunotype B in the combined cohort in 5-year OS (76.0% vs. 44.0%, respectively; P < 0.001) and 5-year DFS (62.8% vs. 48.3%, respectively; P < 0.001). The statistically significant clinical factors identified by univariate Cox regression analysis (Supplementary Table S5) were included to perform multivariate Cox regression analysis. Table 1 presented that stromal immunotype remained a powerful and independent prognostic indicator for OS and DFS in the training, testing, and validation cohort.
Variables . | HR (95% CI) . | P . |
---|---|---|
OS | ||
Training cohort (n = 118) | ||
pN (pN+ vs. N0+Nx) | 4.088 (2.077–8.048) | 0.001 |
Immunotype (type B vs. type A) | 6.815 (2.419–19.201) | <0.001 |
Testing cohort (n = 140) | ||
Grade (High vs. PUNLMP+low grade) | 2.606 (1.249–5.438) | 0.011 |
Immunotype (type B vs. type A) | 5.029 (2.572–9.836) | <0.001 |
Validation cohort (n = 287) | ||
pT (pT3+pT4 vs. pT2) | 1.791 (1.032–3.107) | 0.039 |
pN (pN+ vs. N0+Nx) | 2.223 (1.423–3.473) | <0.001 |
Immunotype (type B vs. type A) | 2.065 (1.334–3.198) | 0.001 |
DFS | ||
Training cohort (n = 118) | ||
pN (pN+ vs. N0+Nx) | 3.942 (1.836–8.462) | <0.001 |
Immunotype (type B vs. type A) | 4.382 (1.666–11.525) | 0.003 |
Testing cohort (n = 140) | ||
pN (pN+ vs. N0+Nx) | 2.264 (1.029–4.982) | 0.043 |
Grade (High vs. PUNLMP+low grade) | 3.920 (1.572–9.776) | 0.004 |
Immunotype (type B vs. type A) | 3.970 (2.074–7.601) | <0.001 |
Validation cohort (n = 287) | ||
pT (pT3+pT4 vs. pT2) | 2.067 (1.324–3.227) | 0.001 |
pN (pN+ vs. N0+Nx) | 1.948 (1.354–2.802) | <0.001 |
Immunotype (type B vs. type A) | 1.675 (1.176–2.386) | 0.004 |
Variables . | HR (95% CI) . | P . |
---|---|---|
OS | ||
Training cohort (n = 118) | ||
pN (pN+ vs. N0+Nx) | 4.088 (2.077–8.048) | 0.001 |
Immunotype (type B vs. type A) | 6.815 (2.419–19.201) | <0.001 |
Testing cohort (n = 140) | ||
Grade (High vs. PUNLMP+low grade) | 2.606 (1.249–5.438) | 0.011 |
Immunotype (type B vs. type A) | 5.029 (2.572–9.836) | <0.001 |
Validation cohort (n = 287) | ||
pT (pT3+pT4 vs. pT2) | 1.791 (1.032–3.107) | 0.039 |
pN (pN+ vs. N0+Nx) | 2.223 (1.423–3.473) | <0.001 |
Immunotype (type B vs. type A) | 2.065 (1.334–3.198) | 0.001 |
DFS | ||
Training cohort (n = 118) | ||
pN (pN+ vs. N0+Nx) | 3.942 (1.836–8.462) | <0.001 |
Immunotype (type B vs. type A) | 4.382 (1.666–11.525) | 0.003 |
Testing cohort (n = 140) | ||
pN (pN+ vs. N0+Nx) | 2.264 (1.029–4.982) | 0.043 |
Grade (High vs. PUNLMP+low grade) | 3.920 (1.572–9.776) | 0.004 |
Immunotype (type B vs. type A) | 3.970 (2.074–7.601) | <0.001 |
Validation cohort (n = 287) | ||
pT (pT3+pT4 vs. pT2) | 2.067 (1.324–3.227) | 0.001 |
pN (pN+ vs. N0+Nx) | 1.948 (1.354–2.802) | <0.001 |
Immunotype (type B vs. type A) | 1.675 (1.176–2.386) | 0.004 |
Abbreviations: CI, confidence interval; HR, hazard ratio; PUNLMP, papillary urothelial neoplasm of low malignant potential.
Stromal immunotype B predict patients' response to adjuvant chemotherapy
After the combination of all three cohorts together, adjuvant chemotherapy could not enhance OS or DFS, no matter in pT2 stage or pT3+pT4 stage in combined cohort (Fig. 4A–D). Interestingly, after grouping patients into different immunotypes, adjuvant chemotherapy significantly enhanced OS and DFS in stromal immunotype B subgroup patients with pT3+pT4 stage (P = 0.016 and P = 0.006) (Fig. 4E). Also, adjuvant chemotherapy significantly enhanced OS and DFS in stromal immunotype A subgroup with pT2 stage (P = 0.038 and P = 0.036). However, further evaluation of subgroup interaction showed that only patients in pT3+pT4 stage with OS were significantly discriminated by stromal immunotype (P = 0.004). Taken together, these results indicated that stromal immunotype could successfully identify patients with pT3+pT4 stage MIBC who were candidates for adjuvant chemotherapy.
Identification of stromal immunotype associated biological pathways and immune checkpoint molecules
Three important cytokine expression profiles (IL2, IFNγ, and TGFβ) were compared in TCGA cohort between immunotype A and B. IFNγ expression was significantly higher in immunotype A subgroup while TGFβ expression was significantly higher in immunotype B subgroup (Fig. 5A). The gene-set enrichment analysis (GSEA) showed that immunotype A subgroup was highly enriched in T-cell receptor signaling pathway and NK-cell–mediated immunity (Fig. 5B). Evaluation of the Philadelphia neoadjuvant chemotherapy (NAC) cohort data showed that NK cells significantly increased, whereas M2 macrophages significantly decreased after chemotherapy (Fig. 5C). Three immunotype A patients changed into immunotype B, and three immunotype B patients changed into immunotype A after chemotherapy (Fig. 5D).
Furthermore, we found that the expression of several immune checkpoint molecules (PD-L1, PD-1, and CTLA-4) was significantly higher in immunotypes A in TCGA cohort (Fig. 5F; Supplementary Table S6). After we categorized the training and the testing cohort into the three major T-cell categories, we found that immunotype A subgroup matched with higher inflamed tumor and less immune desert tumor (Fig. 5G). The tumor mutation burden (TMB) was not shown to be significantly different between different immunotypes (Fig. 5E).
Discussion
Integrating multiple biomarkers into a single model could improve the prognostic value over that of a single model substantially (21–23). Patients are of great heterogeneity in clinical outcomes even with the identical TNM stage. In this study, new stromal immunotype comprising five selected immune features in the stroma of MIBC was identified, as a prognostic tool independent of other clinical characteristics. In contrast with other studies evaluating immune cells in MIBC, our stromal immunotype was constructed from evaluating immune cells thoroughly based on location. IIstromal was built with LASSO Cox regression models, which significantly improve its predictive value. Compared with traditional TNM staging, the immunotype indicates the bioimmunologic information of MIBC and indicate different information from TNM staging. In comparison with the flow cytometry method, using histopathologic analyses to identify immune cells infiltration could provide location information that gave us a clear preview of immune contexture (24). In our study, immune cell infiltration in different locations showed distinctive density and clinical implications. In the ultimate formula, Tregs, mast cells, and macrophage stood for worse prognostic indicators, whereas NK cells and CTLs stood for favorable prognostic indicators in tumor stroma. Also, we found that the T-cell–mediated pathway and NK-cell functions were increased in the immunotype A patients, which might explain the prognostic difference between different immunotype patients. These findings were also consistent with previous findings (12, 25).
Interestingly, in the IHC cohorts, the macrophages and Tregs appear to be coexpressed in one subset of tumors, and the mast cells, CTLs, and NK cells each appear to occupy their own clusters, while in the TCGA/CIBERSORT analyses, the Tregs and macrophages appear to occupy their own clusters, and a relatively large fraction of the CTLs and NK cells appear to coexist. The inconsistencies might be caused by different evaluation methods between IHC and CIBERSORT analysis. Also, all the patients came from China in the IHC cohorts, whereas patients in the TCGA cohort were of different races. Recent studies have revealed the heterogeneity of immune-related gene in different races (26, 27). The different genetic and environmental factors might affect the elements of cancer immunity in MIBC. However, the distribution of immunocytes of subgroup A and subgroup B in three cohorts support the same conclusion: the subgroup A consisted of a CTLhighNKhighTreglowMacrophagelowMClow group associated with a better 5-year OS and DFS, whereas the subgroup B correlated with a worse clinical outcome.
A durable tumor-targeting immune response was essential for the success of cytotoxic chemotherapy (28). In current clinical practice, adjuvant chemotherapy is provided to those high-risk patients (pT3+pT4 or N+) who do not undergo neoadjuvant chemotherapy, but the debate is still ongoing. In our study, although ACT did not show survival advantage in OS and DFS, stromal immunotype stratified high-risk patients (pT3+pT4) who could benefit from adjuvant chemotherapy. Interestingly, patients with stromal immunotype B, who represented a CTLlowNKlowTreghighMacrophagehighMChigh group, have worse clinical outcomes but better response to platinum-based chemotherapy. A recent study unveiled that chemotherapy reprogrammed tumor-associated macrophage from protumor role to antitumor role (29). We found that NK cells increase and M2 macrophages decrease significantly after chemotherapy. The immunotype B patients have higher infiltration of macrophage and low infiltration of NK cells, which might facilitate the response to the ACT. However, further basic study was required to explore the interaction between stromal cells and platinum-based regimens. By far, we demonstrated that our stromal immunotype is both a prognostic and a predictive tool in patients with pT3+pT4 disease.
In a recent phase II trial of gemcitabine + cisplatin + ipilimumab in patients with metastatic urothelial cancer, the combination of chemotherapy and immunotherapy presented with promising therapeutic potential with response rate as high as 64% (23/36; ref. 30). Besides immune checkpoint molecule expression, the TMB and the major T-cell categories might also indicate the response of immunotherapy in MIBC (31). In this study, we found that stromal immunotype A was associated with higher immune checkpoint molecule expression and more inflamed tumor. These results indicated that, although patients with stromal immunotype A were not sensitive to chemotherapy, immune checkpoint inhibitor might work well in this subgroup. Therefore, stromal immunotype might also be used as a predictor of upcoming popularity of immunotherapy.
Several limitations should be addressed in this study. It was a retrospective study and the sizes of our cohorts were relatively small. Therefore, larger independent and prospective validation is needed in the future. Also, we selected several immune cells from CIBERSORT instead of evaluating them all. Other important immune phenotypes might be missed, such as neutrophil infiltration (32). However, simply evaluating all the immune cells might be costly and wasteful in the clinical practice. Better screening method could be used to identify those important immune markers and fulfill the immunotype classification.
In conclusion, this is the first study evaluating the whole immune landscape of MIBC with three independent cohorts. We identified two immunotypes by IIstromal. Moreover, the immunotypes could be a useful prognostic and predictive tool to identify patients with MIBC who might benefit from adjuvant chemotherapy. Thus, the stromal immunotypes might have crucial implications for the postoperative personalized follow-up and decision-making regarding individualized adjuvant treatment.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Disclaimer
The study sponsors have no roles in the study design, in the collection, analysis, and interpretation of data.
Authors' Contributions
Conception and design: H. Fu, Y. Zhu, D. Ye, J. Xu
Development of methodology: H. Fu, Y. Wang
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): H. Fu, Y. Zhu, Y. Wang
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): H. Fu, Y. Zhu, Y. Wang
Writing, review, and/or revision of the manuscript: H. Fu, Y. Zhu, Y. Wang, J. Xu
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Zhu, Z. Liu, J. Zhang, H. Xie, Q. Fu, B. Dai
Study supervision: D. Ye, J. Xu
Acknowledgments
This study was funded by grants from National Natural Science Foundation of China (81402082, 81471621, 81472227, 81671628, and 31770851), Shanghai Municipal Natural Science Foundation (16ZR1406500), and Shanghai Municipal Commission of Health and Family Planning Program (20144Y0223).
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.