Subsets of preoperative sex hormones in testicular germ cell cancer: a retrospective multicenter study

Preoperative homeostasis of sex hormones in testicular germ cell tumor (TGCT) patients is scarcely characterized. We aimed to explore regulation of sex hormones and their implications for histopathological parameters and prognosis in TGCT using a data-driven explorative approach. Pre-surgery serum concentrations of luteinizing hormone (LH), follicle-stimulating hormone (FSH), testosterone (T), estradiol (E2) and prolactin were measured in a retrospective multicenter TGCT cohort (n = 518). Clusters of patients were defined by latent class analysis. Clinical, pathologic and survival parameters were compared between the clusters by statistical hypothesis testing, Random Forest modeling and Peto-Peto test. Cancer tissue expression of sex hormone-related genes was explored in the publicly available TCGA cohort (n = 149). We included 354 patients with pure seminoma and 164 patients with non-seminomatous germ cell tumors (NSGCT), with a median age of 36 years. Three hormonal clusters were defined: ‘neutral’ (n = 228) with normal sex hormone homeostasis, ‘testicle’ (n = 91) with elevated T and E2, low pituitary hormones, and finally ‘pituitary’ subset (n = 103) with increased FSH and LH paralleled by low-to-normal levels of the gonadal hormones. Relapse-free survival in the hormonal subsets was comparable (p = 0.64). Cancer tissue expression of luteinizing hormone- and follicle-stimulating hormone-coding genes was significantly higher in seminomas, while genes of T and E2 biosynthesis enzymes were strongly upregulated in NSGCT. Substantial percentages of TGCT patients are at increased risk of sex hormone dysfunction at primary diagnosis before orchiectomy. TGCT may directly influence systemic hormonal homeostasis by in-situ synthesis of sex hormones.

The retrospective cohort study dataset was provided as an Excel sheet.Data import in R was accomplished with the read_xslx() R function (package readxl) (1).The dataset was formatted with an in-house developed script.
Predominantly embryonic cancers, predominantly teratoma, predominantly yolk sac cancers and predominantly chorionic cancers were defined by the cutoff of ≥ 75% of the respective histology in pathological examination.Tumor stages and Lugano clinical stages were expressed in simplified form as 'I', 'II', 'III'.Laboratory parameter values beyond the detection range were set to their respective lower and upper detection limit.
In the current analysis, n = 518 testis cancer patients of the retrospective cohort with ≤ 50% missing variables were included (Supplementary Figure S1).
The testicle cancer dataset (clinical information, expression) of the TCGA pan-cancer project (2) with n = 149 specimens was fetched from cBioPortal.The R import and formatting was accomplished with an in-house developed script.

Soware
Data analysis was done with R version 4.2.3.
Descriptive statistics were computed with DescTools (6), rstatix (7) and ExDA.Statistical hypothesis testing was done for numeric and categorical variables with the packages rstatix (7) and ExDA, and with the survival (8) and survminer (9) packages for survival data.
Principal component analysis (PCA) was performed with the development package clustTools employing the pcaPP algorithm by Croux et al. (10,11).Multi-dimensional correspondence analysis (MCA) was done with the MASS package (12).Hormonal subsets of testicle cancers in the retrospective cohort were defined by maximum-likelihood latent class analysis (LCA) done with poLCA (13,14).
Multi-parameter classifiers of hormonal subsets of testicular cancer were established for the retrospective cohort with the conditional Random Forest algorithm (15,16) provided by the package party (17).For cross-validation, computation of fit statistics and quality control of the Random Forest models, the packages caret (18) and the development package caretExtra.
Supplementary Material and parts of the manuscript were written in the rmarkdown environment (22) with the package bookdown (23).Management of figures, tables, links and references in the rmarkdown document was accomplished with figur.The report was rendered as a html document with knitr (24) and rmdformats (25).

Stascal signicance and mulple tesng adjustment
If not indicated otherwise, statistical testing p values were corrected for multiple testing with the false discovery rate method (FDR) (26) within each analysis task.Effects with FDR-corrected p < 0.05 are considered statistically significant.Effect sizes for categorical variables were defined with the Cramer V effect size statistic with the following cutoffs: ≤ 0.3: weak, 0.3 -0.6: moderate, > 0.6: strong.Effect sizes for numeric variables analyzed with Mann-Whitney test were assessed by r effect size statistic with the following cutoffs: ≤ 0.3: small, 0.3 -0.5: moderate; > 0.5: strong.Effect sizes for numeric variables analyzed with Kruskal-Wallis test were investigated by η 2 effect size metrics with the following cutoffs: ≤ 0.13: weak, 0.13 -0.26: moderate, > 0.26: strong.Accuracy was assessed by Cohen's κ inter-rater reliability measure with the following effect size cutoffs: ≤ 0.4: weak, 0.4 -0.6: moderate, 0.6 -0.8: good, > 0.8: very good.

Explorave data analysis
Numeric variables are presented in the text and tables as medians with interquartile ranges (IQR), ranges and numbers of complete observations.Qualitative variables are presented as percentages and counts of the categories within the complete observation set.Descriptive characteristic of the retrospective and TCGA cohorts was generated with the explore(what = 'table', pub_styled = TRUE) function (package ExDA) and is presented in Table 1, Table 3 and Supplementary Table S2.
Normality of numeric variables was assessed with Shapiro-Wilk test (function explore(), package ExDA and quantile-quantile plots.Since multiple numeric variables were notnormally distributed, non-parametric statistical hypothesis tests were used.

Comparison of seminoma and non-seminomatous germ cell tumors in the retrospecve cohort
Demographic, clinical, pathological, therapy-related and endocrine parameters were compared between the histology types with Mann-Whitney test and r effect size statistic (numeric variables) or χ 2 test with Cramer V effect size statistic (function compare_variables(), package ExDA).The comparison results are presented in Table 2, Figures 1 -3 and Supplementary Figure S2.Differences in survival between the histology types were analyzed with Peto-Peto test (function surv_pvalue(), package survminer, Figure 1).

Co-regulaon of pre-surgery sex hormones in the retrospecve cohort
Co-regulation of pre-surgery sex hormones treated as normalized, median-centered absolute concentrations were investigated by 4-dimensional PCA (10,11) (function reduce_data(kdim = 4, red_fun = 'pca'), package clustTools).The first four components of the PCA explained >90% of the entire hormone dataset variance as investigated by visual analysis of the scree plot (not shown).PCA loadings for the first two major principal components were visualized with the plot(type = 'loadings') method (package clustTools).
Trends or overlaps in strata of sex hormone levels (concentrations stratified by limits of reference ranges) were explored by two-dimensional MCA performed with the mca() function from the MASS package (12).Column factors were visualized with a custom script.
A total of n = 422 observations were available for these analyses.The results of component and correspondence analysis are presented in Supplementary Figure S3.
The optimal number of hormonal subsets was determined by comparing values of Bayesian Information Criterion (BIC) for LCA models with varying number of classes.The LCA solution with k = 3 classes/hormonal subsets displayed the minimal BIC value suggestive of the optimal model fit.Convergence of this model was achieved in n = 1393 iterations out of 5000 iterations in total (Figure 4A).Conditional probabilities for the hormonal subset assignment estimates for the hormone concentration strata by the final LCA model are presented in Supplementary Figure S4.
Differences in distribution of sex hormone strata between the hormone subsets were assessed by χ 2 test with Cramer V effect size statistic.Differences in blood hormone concentrations expressed as numeric variables between the hormone subsets were investigated by Kruskal-Wallis test with η 2 effect size statistic (both comparison types done with compare_variables(), package ExDA).Differences in sex hormone levels between the hormonal subsets are shown in Figure 4 and Supplementary Table 1.

Demographic and clinical characterisc of the hormonal subsets of the retrospecve cohort
Demographic, clinical, pathological and therapy-related parameters were compared between the hormonal subsets with Kruskal-Wallis test and η 2  Additionally, neutral and pituitary hormonal subset patients were split according to their AFP and HCG status (negative for both AFP/HCG versus positive for AFP or HCG) and differences between the status strata compared with Mann-Whitney test with r effect size statistic or χ 2 test with Cramer V effect size statistic, as appropriate.Significant differences identified in this analysis are presented in Supplementary Figure S6.

Mul-paramater discriminaon between the hormonal subsets of the retrospecve cohort
A multi-parameter model allowing for discrimination between the hormone subsets with solely non-endocrine explanatory variables (available observations: n = 305, Supplementary Figure S7A) was established with the conditional Random Forest algorithm (15-17) using the train(method = 'cforest') wrapper function from the package caret (18).The optimal mtry parameter value obtained by cross-validation-based tuning with Cohen's κ as the tuning metric was 18.The remaining model parameters such as split rule (set to the maximum test statistic), splitting criterion (set to 1.28) and number of random trees (1000) were specified in the control objects returned by the cforest_classical() function from the package party.
Model performance in the training dataset and 10-fold cross-validation was assessed by the accuracy and Cohen's κ statistics computed with the summary() method (package caretExtra) and are presented in Supplementary Figure S7B.Confusion matrices of the predicted and actual hormone subset assignment were visualized with the plot(type = 'confusion') method provided by the caretExtra package.Permutation importance for explanatory variables was computed as described for the histology types and is presented in Supplementary Figure S7A.

Expression of sex hormone-related genes in the TCGA cohort
Expression of 14 sex hormone-related genes in the tumor tissue of the TCGA cohort specimens was investigated in the current report (GNRH1, GNRH2, PRL, CGA, FSHB, LHB, CYP11A1, CYP17A1, HSD17B3, HSD3B1, HSD3B2, CYP19A1, HSD17B1, SHBG).Differences in gene expression between seminoma and NSGCT samples were investigated by Mann-Whitney U test with r effect size statistic (function compare_variables(), package ExDA).

Figure S2. Preoperative absolute serum concentrations of cancer markers in seminoma and non-seminomatous germ cell tumors.
Whitney test with r effect size statistic.P values corrected for multiple testing with the false disocvery rate method.Blood concentrations of sex hormones (T: testosterone, E2: estradiol, FSH: follicle stimulating hormone, LH: luteinizing hormone; PRL: prolactin) as normalized median-centered numeric variables and as clinical strata were subjected to 4-dimensional principal component analysis (PCA, A) and multidimensional correspondence analysis (MCA, B), respectively.Loadings (PCA) and column factors (MCA) for the first two dimensions are presented in scatter plotsplots.Numbers of observations are indicated in the plot captions.(B) Hormonal subset conditional response probabilities estimated by the final LCA model presented as heat maps.Numbers of observations in the subsets and percentage within the hormone dataset are indicated in the plot captions.shuffled.Accuracy differences are displayed in a bar plot.Numbers of samples in the hormonal subsets are shown in the plot caption.(B) Discriminant performance of the Random Forest model in the original data set and 10fold CV was visualized as heat maps of confusion matrices.The overall subset assignment accuracy and Cohen's κ are displayed in the plot captions.Fractions of observations correctly assigned to the neutral, testicle and pituitary subsets are indicated in X axes.HCG: human chorionic gonadotropin; NSGCT: non-seminomatous germ cell tumors; T replacement: testosterone replacement; AFP: alpha fetoprotein; LDH: lactate dehydrogenase; LVI: lymphovascular invasion; SemiCa: percentage of seminoma histology; CS Lugano: Lugano class; RLA: retroperitoneal lymphadenectomy; YSCa: percentage of yolk sac histology; ECa: percentage of embryonic histology; RTI: rete testis invasion.
Differences in serum concentrations of alpha fetoprotein (AFP), human chorionic gonadotropin (HCG) and activity of lactate dehydrogenase (LDH) between seminoma and non-seminomatous germ cell tumors (NSGCT) were assessed by Mann-Whitney test with r effect size statistic.P values were corrected for multiple testing with the false discovery rate method.Concentrations are presented in violin plots wit single observations depicted as points and medians with interquartile ranges represented by red diamonds and whiskers.Effect sizes and p values are displayed in the plot captions.Numbers of complete observations are indicated in the X axes.