Radiogenomics Consortium Genome-Wide Association Study Meta-Analysis of Late Toxicity After Prostate Cancer Radiotherapy

Abstract Background A total of 10%–20% of patients develop long-term toxicity following radiotherapy for prostate cancer. Identification of common genetic variants associated with susceptibility to radiotoxicity might improve risk prediction and inform functional mechanistic studies. Methods We conducted an individual patient data meta-analysis of six genome-wide association studies (n = 3871) in men of European ancestry who underwent radiotherapy for prostate cancer. Radiotoxicities (increased urinary frequency, decreased urinary stream, hematuria, rectal bleeding) were graded prospectively. We used grouped relative risk models to test associations with approximately 6 million genotyped or imputed variants (time to first grade 2 or higher toxicity event). Variants with two-sided Pmeta less than 5 × 10−8 were considered statistically significant. Bayesian false discovery probability provided an additional measure of confidence. Statistically significant variants were evaluated in three Japanese cohorts (n = 962). All statistical tests were two-sided. Results Meta-analysis of the European ancestry cohorts identified three genomic signals: single nucleotide polymorphism rs17055178 with rectal bleeding (Pmeta = 6.2 × 10−10), rs10969913 with decreased urinary stream (Pmeta = 2.9 × 10−10), and rs11122573 with hematuria (Pmeta = 1.8 × 10−8). Fine-scale mapping of these three regions was used to identify another independent signal (rs147121532) associated with hematuria (Pconditional = 4.7 × 10−6). Credible causal variants at these four signals lie in gene-regulatory regions, some modulating expression of nearby genes. Previously identified variants showed consistent associations (rs17599026 with increased urinary frequency, rs7720298 with decreased urinary stream, rs1801516 with overall toxicity) in new cohorts. rs10969913 and rs17599026 had similar effects in the photon-treated Japanese cohorts. Conclusions This study increases the understanding of the architecture of common genetic variants affecting radiotoxicity, points to novel radio-pathogenic mechanisms, and develops risk models for testing in clinical studies. Further multinational radiogenomics studies in larger cohorts are worthwhile.

Associations between pairs of toxicities were assessed by hazard ratios, considering each toxicity as a time-dependent covariate in a Cox model for each other toxicity, unadjusted for any other predictor. If the explanatory toxicity was censored before the dependent toxicity, the dependent toxicity was artificially censored at the same earlier time.

Genotype Imputation
Genetic data were imputed using, as reference haplotypes, the1000 Genomes Project Phase 3 (Haplotype release date October 2014) for chromosomes 1 to 22 and the 1000 Genomes Project Phase 1 (Haplotype ChrX release date Aug 2012) for chromosome X, since the phased data for Chr X from 1000GP Phase 3 was not available. A two-stage procedure used SHAPEIT (shapeit.v2.r790.Ubuntu_12.04.4.static) to derive phased genotypes (default parameters with the following modifications: 10 burn-in iterations, 10 pruning iterations, and 50 iterations to compute transition probabilities) and IMPUTEv2 (impute_v2.3.2_x86_64_static) to perform imputation of the phased data (default parameters with the following modifications: 5Mb nonoverlapping intervals, 800 reference haplotypes to use as templates when imputing missing genotypes, and 500kb buffer region). 1000 Genomes Project variants whose minor allele frequency in Europeans and East Asians was lower than 0.001 were excluded from imputation. All OncoArray datasets were imputed jointly; the Affymetrix and Illumina CytoSNP12 datasets were imputed separately following the same procedure.

Fine-scale mapping
Genomic regions were defined as the 1Mb interval surrounding each statistically significant independent association. We re-imputed genotypes for the non-directly-genotyped variants using IMPUTE2 [15] and a reference panel [16] using the standard IMPUTE2 MCMC algorithm for follow-up imputation (see [17] for detailed description of the parameters used) to improve accuracy at low frequency variants. Variants with imputation info score ≥0.3 in all cohorts and MAF ≥0.02 in at least one cohort were included. 4,190 variants across the chr1:230337180-231337180 region; 3,776 at chr5:156903410_157903410 and 3,987 at chr9:30366808-31366808 were evaluated for hematuria, rectal bleeding or decreased urinary stream risk, respectively.
Then, the most statistically significant variant (index variant at signal 1) was used to perform conditional analysis in each cohort independently. To define the cumulative posterior probability of the credible set, we estimated the empirical Bayes Factor [18]. with more than one independent signal Bayes Factor was estimated using the summary statistics from the conditional analysis, after adjusting for other index variants at the region.

Credible causal variant (CCV) annotation
Variants were annotated with Variant Effect Predictor [19] to determine their effect on genes, transcripts, and protein sequences. To evaluate whether CCVs were located at regulatory regions, we overlapped our CCVs with Encode enhancer-like and promoter-like regions for 73 tissues and cells (primary, immortalized, in vitro differentiated) with available data for both enhancer-and promoter-like regions ( [19][20][21] and DCC accession: ENCSR037HRJ; GEO accession: GSE30567). In order to evaluate whether the CCVs could drive the expression of local genes, we accessed the GTEx Portal on 04/19/2018 to retrieve the metasoft results for all tissues in the V7 release. LocusZoom [22] was used to visualize associations for regions containing CCVs. Linkage disequilibrium was estimated using as reference the European ancestry populations from the 1000 Genomes Project (Phase 3, release 20130502; [16]).

Pathway Analysis
Gene-and pathway-based analysis was performed using Pascal (Pathway scoring algorithm) [23].
Gene-based scores were computed using the default "sum" option, which calculates the sum of chi-squared statistics and measures the strongest association signal per gene, respectively. SNPs were mapped to genes using a 100kb window surrounding each gene. Pathway-based scores are computed using a modified Fisher method, which improves statistical power compared with enrichment-based analysis while maintaining rigorous type I error control. The KEGG, Biocarta, and Reactome databases were queried for the pathway-based analysis.

Multivariable Modeling
Clinical variables were combined with genetic variants (identified via GWAS meta-analysis and from prior studies) using cohort-stratified grouped relative risk models, assuming an additive model for each variant, resulting in per allele hazard ratios (HR). Such grouped relative risk models estimate hazard ratios based on grouped survival data, assuming proportional hazards for the latent continuous survival times within each cohort stratum. The discrete monitoring times need not be equally spaced, nor need they be the same across cohort strata, but it is necessary that they be on the same temporal grid for all subjects within each cohort stratum. Confidence intervals and p-values were likelihood based and two-sided, with p-values ≤ 0.05 considered statistically significant.
Stepwise model selection was used to identify a parsimonious multivariable model for each toxicity outcome. For the multivariable cohort-stratified grouped relative risk models presented in Table 4, an alpha to enter of 0.10 and an alpha to stay of 0.05 were used as model selection parameters. Genetic variants were forced into the model in advance of the inclusion of any clinical variables. The large number of tied follow-up times were handled using the exact method, equivalent to marginal likelihood (that Efron's method approximates) -not the discrete time method that assumes proportional odds rather than proportional hazards.
The log 2 transformation was used to symmetrize the distribution of strongly positively skewed continuous variables, thus reducing the influence of the most extreme observed covariate values and resulting in a hazard ratio per doubling of the predictor. The functional form of each continuous variable was chosen via model selection from the following options: linear; piecewise constant histospline with knots at one or more quartiles; or piecewise linear spline with knots at one or more quartiles, with the option to force the slope to be 0 to the left of the first knot (as a reference group, similar to that of a histospline).
Missing data were imputed within cohorts as follows. If a variable had ≤25% of values missing within a cohort, within-cohort mean imputation was used to impute the missing values. If a variable had >25% of values missing within a cohort, the variable was set to a constant of 0 within the cohort, allowing the hazard ratio for that variable to be estimated based only on cohorts with no more than 25% missing data, without requiring subjects missing data on a subset of variables within some cohorts to be excluded entirely from the analysis.
This novel approach allowed us to at least partially adjust for variables that were available in some cohorts but not others, where the adjustment would be complete for variables that truly did not vary within cohorts in which they were missing -irrespective of the true constant value within each such cohort.
In addition to grouped relative risk regression, two other multivariable modeling methods were applied to derive separate predictive models for each of the four toxicity endpoints: Polygenic Risk Score (PRS) and a machine-learning method. In both methods, a "training" set was used to derive the model and a "test" set was used to evaluate model performance such that the training set was independent of the test set. The training set included a randomly selected 50% of the RAPPER study participants and all other cohorts; the testing set included the 50% of the RAPPER study participants not included in training data.
The first method, PRS, is a linear combination of risks of multiple SNPs identified by GWAS. The risk SNPs comprising the PRS were identified via GWAS meta-analysis of results from the cohorts comprising the training set for each of the four toxicity endpoints (rectal bleeding, increased urinary frequency, decreased urinary stream, and hematuria). On the training data, we tested several P meta thresholds (1E-1, 1E-2, …, 1E-8) in selecting SNPs for PRS, followed by LD-pruned using 1000 Genomes Project EUR panel. Afterward, we computed PRS for each individual of the testing data, and examine the association between PRS and toxicity endpoint.
The machine learning-based method, which was previously developed and applied to RNAseq-derived gene expression data [24], was used to derive multi-SNP models for predicting the toxicity outcomes considered in this study. Since this method was designed for binary classification tasks, we focused on the binarized versions of these endpoints (grade 2 or worse toxicity vs. grade 0 or 1 toxicity) in these experiments, as in the GWAS meta-analysis. We applied this method to the training set, with the constituent SNPs filtered at the same thresholds used for PGS, and evaluated the resultant models on the test set.

Data Management and Analysis
Genomic data were formatted using R (version 3.2.2, R Foundation for Statistical Computing, Vienna, Austria). Analysis was carried out using ProbABEL [25], which employs the coxfit2 function in the R package survival. GWAS results were meta-analyzed using Stata (version 14.2, StataCorp LLC, College Station, TX).
Multivariable modeling was done using SAS (version 9.4, SAS Institute, Cary, NC). Pascal was used to compute gene and pathway scores [23].

Supplementary Tables
Supplementary Table 1. Definitions of toxicity grades. Semicolons indicate "or". CTCAEv3.0 -Urinary retention 0 = No toxicity 1 = Hesitancy or dribbling, no significant residual urine; retention occurring during the immediate postoperative period 2 = Hesitancy requiring medication; or operative bladder atony requiring indwelling catheter beyond immediate postoperative period but for <6 weeks 3 = More than daily catheterization indicated; urological intervention indicated (e.g., TURP, suprapubic tube, urethrotomy) 4 = Life-threatening consequences; organ failure (e.g., bladder rupture); operative intervention requiring organ resection indicated Gene-PARE NTMC American Urological Association Symptom Score Q5 ¶ 0 = Had a weak urinary stream 'not at all' or 'less than 1 time in 5' 1 = Had a weak urinary stream 'less than ½ the time' 2 = Had a weak urinary stream 'about ½ the time' or 'more than ½ the time' 3 = Had a weak urinary stream 'almost always' Hematuria #

RAPPER
LENT Institutional scale 0 = None 2 = Bleeding present * Increased urinary frequency was not analyzed in CCI-BT because pre-radiotherapy assessments were more than one year prior to starting radiotherapy for the majority of participants. Abbreviations: OPD, outpatient department; NA, grade not applicable † UGhent used an institutional-specific scaled based on CTCAEv3.0 ‡ During the past month or so, how often have you had to urinate again less than two hours after you finished urinating? § Over the past month, how many times per night did you most typically get up to urinate from the time you went to bed at night until the time you got up in the morning? || Decreased urinary stream was not analyzed in CCI-BT because pre-radiotherapy assessments were more than one year prior to starting radiotherapy for the majority of participants. ¶ During the past month, how often have you had a weak urinary stream? # Endpoint not available in CCI-BT, UGhent or NTMC ** Rectal bleeding was not assessed in NTMC or PRRG. Rectal bleeding was assigned a single grade in GenePARE using information across all follow up assessments, and so this endpoint was not available for analysis using Cox proportional hazards modeling. 38.1 * Minor allele frequency is from PRACTICAL oncoarray samples of European ancestry. Abbreviations: SNP, single nucleotide polymorphism; MAF, minor allele frequency; HR, hazard ratio; CI, confidence interval; BFDP, Bayesian false discovery probability; NA, not available; INS, insertion. † Imputation info score values are from the oncoarray. ‡ Hazard ratio corresponds to the minor allele with the major allele treated as the reference group. § Two-sided P meta was calculated using a Wald test. || BFPD estimated assuming a prior variance, W = 0.32^2, and prior probability of a non-null association 0.0001 ¶ Hazard ratio is for the major allele with the minor allele treated as the reference group # Beta coefficient from linear regression of STAT score at 2 years after radiotherapy  † Two-sided P meta was calculated using a Wald test. ‡ Two-sided P meta was calculated using a chi-square test. * Two-sided P meta was calculated using a Wald test. 0.05 * Two-sided pvalue was calculated using a chi-square test. † Following methods in [23], a Monte Carlo estimate of the p-value is obtained by sampling random gene sets of size and calculating the fraction of sets reaching a higher score than gene set of the given pathway. ‡ Pathway is associated with both rectal bleeding and increased urinary frequency § Pathway is associated with both rectal bleeding and hematuria || Pathway is associated with both rectal bleeding and decreased urinary stream ¶ Pathway is associated with both increased urinary frequency and hematuria # Pathway is associated with both decreased urinary stream and hematuria