Article Text
PDF
Abstract
Background The risk of developing advanced neoplasia (AN; colorectal cancer and/or high-grade dysplasia) in ulcerative colitis (UC) patients with a low-grade dysplasia (LGD) lesion is variable and difficult to predict. This is a major challenge for effective clinical management.
Objective We aimed to provide accurate AN risk stratification in UC patients with LGD. We hypothesised that the pattern and burden of somatic genomic copy number alterations (CNAs) in LGD lesions could predict future AN risk.
Design We performed a retrospective multicentre validated case–control study using n=270 LGD samples from n=122 patients with UC. Patients were designated progressors (n=40) if they had a diagnosis of AN in the ~5 years following LGD diagnosis or non-progressors (n=82) if they remained AN-free during follow-up. DNA was extracted from the baseline LGD lesion, low-coverage whole genome sequencing performed and data processed to detect CNAs. Survival analysis was used to evaluate CNAs as predictors of future AN risk.
Results CNA burden was significantly higher in progressors than non-progressors (p=2×10−6 in discovery cohort) and was a very significant predictor of AN risk in univariate analysis (OR=36; p=9×10−7), outperforming existing clinical risk factors such as lesion size, shape and focality. Optimal risk prediction was achieved with a multivariate model combining CNA burden with the known clinical risk factor of incomplete LGD resection. Within-LGD lesion genetic heterogeneity did not confound risk prediction.
Conclusion Measurement of CNAs in LGD is an accurate predictor of AN risk in inflammatory bowel disease and is likely to support clinical management.
- INFLAMMATORY BOWEL DISEASE
- COLORECTAL CANCER
- ULCERATIVE COLITIS
- CANCER PREVENTION
- GENETIC INSTABILITY
Data availability statement
All data relevant to the study are included in the article or uploaded as supplementary information. Processed data to recreate all modelling results is included in online supplemental file 2 (online supplemental table S2). Copy number segmentations and CNA calls for all 270 samples are available on Mendeley Data (doi:10.17632/vm7wjjb36y.1). Raw sequencing data are not available due to the retrospective and anonymous nature of the study. Code for reproducing main results is publicly available at https://github.com/yosoykit/C4L_CNA.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
Patients with ulcerative colitis (UC) are at increased risk of colorectal cancer (CRC), and this risk increases dramatically in those who develop low-grade dysplasia (LGD). However, there is currently no accurate way to risk-stratify UC patients with LGD, leading to both overtreatment and undertreatment.
WHAT THIS STUDY ADDS
We show that the burden of somatic copy number alterations (CNAs) within resected LGD lesions from UC patients strongly predicts future cancer and/or high-grade dysplasia (HGD) development.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY
The measurement of CNAs in LGD lesions arising in colitis is a robust and low-cost biomarker of HGD/CRC risk in inflammatory bowel disease that can be readily translated into clinical practice. Future clinical use of the biomarker should enable more effective management of LGD in colitis, preventing overtreatment and undertreatment of cancer risk.
Introduction
The management of colitis-associated colorectal cancer (CA-CRC) risk suffers from both chronic underdiagnosis and overdiagnosis, leading to uncertainty in how to best monitor and treat patients. Patients with colitis are offered enrolment into endoscopic surveillance programmes to detect early signs of cancer, namely the pathological diagnosis of dysplasia. Patients with low-grade dysplasia (LGD), itself a contentious definition,1–3 have an approximately 30% chance of progressing to CRC at 10 years4 and this level of risk is considered sufficient to offer prophylactic colon resection. Major surgery has many risks, and subsequent detrimental impact on quality of life due to an ileal-pouch anal anastomosis or a stoma. Alternatively, patients can be offered watchful waiting, where frequent colonoscopy and worry about CRC risk also severely impact quality of life. Patients and clinicians can have disparate thresholds for tolerable cancer risk,5 6 and so the joint decision-making process of which route to follow can prove challenging. Ultimately, these decisions are made against the backdrop where more than two-thirds of colitis patients with resected LGD are not actually at elevated risk of developing CRC, and so any treatment for these patients is overtreatment. Consequently, there is a major unmet need for accurate cancer risk stratification in colitis patients with LGD.
Some level of cancer risk stratification of colitis-associated LGD lesions is achieved by assessing clinical and endoscopic factors. Lesion size greater than 10 mm, non-polypoid lesion shape or endoscopically invisible and a history of indefinite dysplasia are significantly associated with LGD progression to advanced neoplasia (namely high-grade dysplasia (HGD) and CRC).7 We recently found that four clinicopathological variables, endoscopically visible LGD>1 cm, incomplete or impossible endoscopic resection, moderate/severe histological inflammation within 5 years of LGD diagnosis and multifocal dysplasia provided excellent negative predictive value (NPV) of cancer risk.8 Other clinical risk factors known to be associated with progression to CRC in the general inflammatory bowel disease (IBD) population include a concomitant diagnosis of primary sclerosing cholangitis (PSC), cumulative inflammation burden, stricturing, the presence of inflammatory pseudopolyps and colonic scarring.9
The inflammatory microenvironment generated by colitis drives carcinogenesis through genomic and molecular pathway alterations that differ from those seen in sporadic CRC.10–13 One key difference is the earlier onset of aneuploidy in colitis, which can be detected even in non-neoplastic colitic epithelium13–15 and has long been recognised to be associated with progression risk.16 Studies that aim to exploit these molecular differences for risk stratification of LGD remain limited by the small sample sizes,17–19 their use of ‘high risk’ flat or invisible dysplasia lesions, and their reliance on flow cytometry and comparative genomic hybridisation, with further limitations in terms of the large quantities of preserved cells or DNA required.
Low-pass whole genome sequencing (lpWGS) is a cost-efficient approach that uses next-generation sequencing to reliably detect copy number alterations (CNAs). The input can be minimal amounts of DNA (as low as 500 pg), it is compatible with the highly degraded DNA that is typically derived from formalin-fixed paraffin-embedded (FFPE) blocks in hospital archives,20 and it does not require a matched patient-specific reference genome for copy number estimation. In this study, we show that CNAs detected in lpWGS data from LGD biopsies of IBD patients very strongly predict future risk of HGD and/or CRC, and that this genomic measurement adds significant prognostic value to the currently used clinicopathological variables.
Results
Copy number profiles can be derived from archival LGD samples by lpWGS
We conducted a retrospective case–control study (figure 1) in which ‘progressors’ (cases) were defined as patients with IBD where LGD was the highest histological grade detected at colonoscopy and who subsequently developed HGD or CRC within 5 years. The ‘non-progressors (NP)’ (controls) were defined as patients with IBD where LGD was the highest histological grade detected at colonoscopy, and who subsequently remained free from both HGD and CRC for approx. 5 years or more.
Study design. Schematic showing the detection and removal of low-grade dysplasia during routine IBD surveillance, with progressors developing high-grade dysplasia and/or colorectal cancer within 5 years (top), and non-progressors remaining free from these (bottom). CRC, colorectal cancer; IBD, inflammatory bowel disease.
We identified a discovery cohort of 67 patients (22 ‘progressors’ who developed HGD/CRC a median 427 days after LGD diagnosis and 45 ‘NP’) bearing 78 LGD lesions (see table 1 and Methods for details). While progressors and NP patients in the discovery cohort were well matched by age, gender, duration of disease, PSC status and lesion location and morphology, the progressor cohort is statistically over-represented by larger lesions. In the progressor patients, LGD was found at the same anatomical location as the future HGD/CRC in 20 out of 22 cases, strongly suggesting that resection of these LGD lesions left residual ‘cancerised’ cells.10 21 A validation cohort of 55 patients was collected from University of Newcastle (n=15), University of Belfast (n=4) and Oxford University (n=36) that included a total of 59 LGD lesions (see table 2).
Summary of patients and LGD samples used as the discovery data by progression status, clinical features and endoscopic findings
Summary of patients and LGD samples used as the validation data by progression status, clinical features and endoscopic findings
Where patients had multiple LGD lesions, we extracted DNA separately from each using epithelial enrichment by laser capture microdissection or needle macrodissection. To examine intralesion heterogeneity from the largest LGD lesions, we extracted DNA from multiple regions within the same lesion. In total, we extracted DNA from 164 regions of LGD in the discovery cohort and from 106 regions of LGD in the validation cohort. We successfully performed lpWGS on all regions, with median coverage of 0.11× (IQR 0.09–0.13×). For each region, we derived a whole-genome copy number profile (see the Methods section).
Progressor LGD has a higher CNA burden than NP LGD
We used lpWGS data to quantify the burden of somatic CNAs in each LGD region. We noted that lesions which appear to be endoscopically and histologically similar can have strikingly different CNA profiles (figure 2).
Representative endoscopic and histological images with copy number profiles. Endoscopy images (upper right) and H&E (upper left) with genome-wide copy number profile below for representative non-progressor (A) and progressor (B) patients.
We compared CNA burden between progressor (P) and NP lesions and observed a large difference in the number of CNA events between the two groups in the discovery dataset (figure 3A), and this was highly statistically significant (figure 3B, median 4 CNAs in NPs vs 16 CNAs in Ps, p=2×10−6, Mann-Whitney U test). Similarly, significant differences were present in the validation cohort (online supplemental figure S1). We note that for patients who had multiple regions or multiple LGD lesions analysed, we considered only the region with highest CNA burden for this comparison. We compared the burden of CNAs detected after epithelial enrichment by laser capture microdissection versus needle macrodissection and found no consequential difference (online supplemental figure S2). There was no correlation between CNA burden and calendar year of LGD resection (online supplemental figure S3).
Supplemental material
Genomic alterations in non-progressor (NP) vs progressor LGD in discovery cohort. (A) Heatmap of genome-wide copy number alterations for lesions in the discovery LGD cohort, sorted by percent genome altered. Sporadic adenomas included below. (B) Violin plots showing the number of altered genomic segments in progressor and NP lesions. (C) Genome-wide CNA frequency for NP (top) and progressor (P, bottom) patients. Stars indicate significant arm level differences at adjusted p<0.01 level. For patients with multiregion analysis, the most highly altered sample per patient was included. LGD, low-grade dysplasia; MSI, microsatellite instability; PGA, percent genome altered; WGD, whole genome doubling.
We next compared the CNA profiles of our LGD lesions to those of 19 sporadic adenomas from patients without IBD.22 We found that although there were commonly altered regions (such as gains of chromosome 7, 13 and 20), there were alterations unique or highly enriched in IBD-LGD (figure 3A).
To determine whether there were specific CNAs associated with progressor LGD, we reviewed the distribution of individual copy number gains and losses across all 67 patients in the discovery dataset (figure 3C). Evidence of selection is indicated by the non-random distribution of CNA frequencies (assuming an equal genome-wide propensity to generate CNAs)—there were CNAs which were never observed in this cohort (eg, chromosome 4 gains and 2q losses), and there were CNAs that are repeatedly detected (such as loss of chr5q and gain of chr7). Many of the common CNAs found across the progressor cohort were also seen in the NP cohort, although at far lower frequency; however, one of the most striking differences between the cohorts was in the frequency of losses involving chromosome 17 (the chromosome where the tumour suppressor gene TP53 is located) which were strongly associated with progressor LGD (Benjamini-Hochberg adjusted p<0.001 for Fisher’s exact test). However, the strongest prognostic signal on chr17 was for loss of the q arm (odds ratio (OR)=40.96, adjusted p<0.001), not for loss of the p arm where p53 is located (OR=18.96, adjusted p<0.001). Therefore, interestingly, the prognostic signal associated with chr17 loss is not likely to simply be a readout of p53 loss. There were also other specific chromosomal CNAs that were found in a significant proportion of progressor LGDs but not in NP LGDs (online supplemental table S1). Interestingly, we noted CNA losses involving subtelomeric segments of chromosomes 4, 5, 6 and 8 in both NPs and progressors, either in isolation (as seen in the former) or as part of a chromosomal arm/whole chromosome loss (as seen in latter, figure 3A).
For 66 lesions in our study that we analysed histopathologically previously,3 we had sequencing data to investigate potential correlation between CNA burden with non-conventional forms of dysplasia. Progressor status was not associated with LGD subtype (Fisher’s exact test p=0.58). CNAs burden was similar in conventional and non-conventional LGD subtypes (online supplemental figure S4). We detected higher numbers of CNAs in invisible versus visible lesions (67 lesions with data included in figure 3), which is consistent with the known increased CRC risk associated with invisible LGD lesions (online supplemental figure S5).
Genomic features of LGD can predict future HGD/CRC
We further interrogated the lpWGS data by calling whole genome doubling and microsatellite instability (MSI) status in each region (see figure 3A). We used lpWGS-derived genomic features (CNA profiles, MSI status) to derive a binary ‘genomic CNA score’ for risk of progression, in which patients are classified as high-risk (score=1) or low-risk (score=0). A patient was considered as ‘high risk’ if (1) the number of CNA segments was >14.5 (75th percentile of CNA segment number in all discovery data) in any LGD sample sequenced, (2) loss on chromosome 17 in any LGD sample sequenced and/or (3) any LGD sample was classified as MSI. If none of these applied, then the ‘genomic CNA score’ for that region was zero and the patient would be considered ‘low risk’.
Univariate logistic regression analysis considering the genomic CNA score and the clinical and endoscopic variables used in the ulcerative colitis cancer risk estimator (UC-CaRE),8 along with a diagnosis of PSC, found that the genomic CNA score was significantly associated with increased risk of progression to HGD/CRC (OR=36.0, p=9e−07; figure 4A). In the discovery cohort, the univariate model with genomics alone (CNA score) had an area under the receiver operating characteristic curve (AUC)=0.85 for accuracy in identifying progressor versus NP patients. When the same criteria for ‘high risk’ were applied to data from the independent validation cohort, we found that a univariate CNA score model had an AUC=0.84 at 5 years post-LGD resection.
Model performance for prediction of future HGD/CRC. (A) ORs for univariate models considering genomic score, UC-CaRE and PSC variables. ‘Unresected LGD’ is defined as either endoscopically invisible LGD (which was detected on histological examination of biopsy only and so by definition cannot be completely resected) or LGD specifically noted to have unclear resection margins. (B) AUC values for genomic score alone at 5 years post-LGD resection and PPV/NPV over time in validation data (all 55 patients included at baseline) with 95% bootstrap CIs (see the Methods section ‘Prediction accuracy’ for more details). Kaplan-Meier progression-free survival for the discovery (C) and validation (D) cohorts using a univariate genomic score model. (E) ORs for multivariate model chosen using stepwise selection. (F) AUC values for multivariate model at 5 years post-LGD resection and PPV/NPV over time in validation data (n=37 patients included at baseline with data on incomplete LGD resection). (G) Kaplan-Meier curves for validation data; high-risk (>0.5 risk) patients determined by multivariate model prediction. AUC, area under the curve; CRC, colorectal cancer; HGD, high-grade dysplasia; LGD, low-grade dysplasia; MSI, microsatellite instability; NPV, negative predictive value; PPV, positive predictive value; PSC, primary sclerosing cholangitis; UC, ulcerative colitis.
Using the genomic CNA score model in the external validation data, we also computed positive predictive value (PPV) and NPV over time since the LGD was resected (see the Methods section). Usually, this corresponded to the time of LGD diagnosis, but when this was not available, we used the earliest available time point with LGD found. We found that PPV=0.27 at 6 months, PPV=0.45 at 1 year and PPV=0.91 at 5 years post index time point (figure 4B). For patients classified as low risk, we found that NPV=1.0 at 6 months, NPV=0.98 at 1 year and NPV=0.90 at 5 years post index time point (figure 4B, see also online supplemental figures S6 and S7 for corresponding results of sensitivity, specificity and AUC across follow-up time).
To further investigate the utility of lpWGS as a clinical risk stratification tool, we used Kaplan-Meier (KM) survival analysis to compare HGD/CRC progression-free survival curves between low-risk and high-risk patients as determined by the CNA score. In the discovery cohort, patients deemed ‘high risk’ were significantly more likely to progress to HGD/CRC than patients deemed ‘low risk’, with an HR of 13.4 (log-rank p<0.0001, figure 4C). Similarly, in the validation cohort, we found that our genomic CNA score effectively stratifies high-risk versus low-risk patients with HR=14.7 (log-rank p<0.0001, figure 4D). We also evaluated per cent genome altered (PGA) as an alternative genomic biomarker of cancer risk but found slightly inferior biomarker performance (online supplemental figure S8).
We evaluated the genomic CNA score model performance in patients with greater than 1 year of follow-up. In this subset of patients, the genomic biomarker still performed extremely well (online supplemental figure S9). In the validation data, we found that PPV=0.33 at 2 years and PPV=0.83 at 5 years post index time point. For patients classified as low risk, we found that NPV=0.98 at 2 years and NPV=0.92 at 5 years post index time point.
A multivariate model combining genomic CNA score and clinical information increases predictive power
To determine a parsimonious multivariate model that combined genomics with candidate clinicopathological variables, we performed stepwise selection using all variables shown in figure 4A for the discovery data. The best-performing model included our genomic CNA score and the UC-CaRE binary variable of incomplete LGD resection (yes or no). Both variables remained significant in multivariate logistic regression although the OR for our genomic CNA score was six times higher than the OR for incomplete resection (figure 4E). In the discovery data, this multivariate model had an AUC=0.89 for accuracy in classifying progressor vs NP patients. Here, we used predicted risk from the multivariate model and considered a patient to be ‘high risk’ if their predicted risk was greater than 0.5 and ‘low risk’ otherwise. When the exact same risk model was applied to the validation data, we found excellent model accuracy with an AUC=0.95 at 5 years post-LGD resection (n=37 patients included with data on incomplete LGD resection). Using this risk prediction model in the external validation data, we found that PPV=0.4 at 6 months, PPV=0.8 at 1 year and PPV=1.0 at 5 years post index time point (figure 4F). For patients classified as low risk, we found that NPV=1.0 at 6 months, NPV=1.0 at 1 year and NPV=0.90 at 5 years post index time point (figure 4F).
Finally, we used KM survival analysis to compare HGD/CRC progression-free survival curves between low-risk and high-risk patients as determined by the combined risk model. The risk stratification in the discovery data was identical to the univariate model results shown in figure 4C. However, when applying to the independent validation data, we found that the combined risk model improves the stratification of high-risk versus low-risk patients in their progression risk with HR=27.0 (log-rank p<0.0001, figure 4G).
Increased intralesion genomic heterogeneity in progressor LGD lesions compared with NP lesions
To characterise intralesion heterogeneity, we compared the CNA profiles of multiregion sequenced LGD lesions. We derived a CNA-based phylogenetic tree for each multiregion sequenced lesion (median 3 regions per lesion, range 2–8 regions per lesion). As previously described, we saw greater CNA burden in progressor LGD versus NP LGD (figure 5A,B). CNA profiles for all regions sequenced can be found in online supplemental figure S10 (discovery data) and online supplemental figure S11 (validation data).
Phylogenetic analysis of multiregion data. Phylogenetic trees (left) and matched copy number profiles (right) for (A) a representative non-progressor (NP) lesion and (B) a representative progressor lesion. Tree statistics for lesions in the discovery cohort that had at least two regions sequenced were maximum tree length (C), minimum tree length (D), clonal CNA (E) and average subclonal CNA (F). CNA, copy number alteration; LGD, low-grade dysplasia.
We measured four tree statistics to compare all progressor (n=17) and NP (n=22) lesion trees (figure 5C–F). Maximum and minimum tree length (the number of events from a diploid state to the most and least aberrant sample) was significantly higher in progressors compared with NPs (figure 5C,D). The maximum and minimum tree length was 15–2 for the NP example (figure 5A) and 57–39 for the progressor example (figure 5B). Average subclonal CNA (ie, average branch length from the trunk to the tips) was also significantly higher in progressor legions implying greater ongoing genomic ‘evolvability’ (figure 5F).
Importantly, the clonal CNA (ie, trunk length) was much higher in progressor lesions (figure 5E and seen when comparing figure 5A,B). Therefore, the detection of an elevated CNA burden in progressor lesions was expected to be robust to intra-lesion heterogeneity, despite the elevated heterogeneity in progressor lesions. To investigate this further, we repeated our analysis of the predictive power of our genomic CNA score, this time choosing a random sample per patient in the validation data rather than the most highly altered region. We calculated similar PPV and NPV values over time since the initial time point, confirming that the predictor is robust to lesion heterogeneity (online supplemental figure S12).
Discussion
CRC risk is chronically underdiagnosed and overdiagnosed in colitis patients, due to lack of a sufficiently accurate predictor of future risk. Here, we have found that the pattern and burden of somatic CNAs within LGD lesions, as measured by standard short-read shallow whole genome sequencing, very strongly predicts the future risk of CRC, with PPV and NPVs at 5 years above 90%. Further, genomic analysis by lpWGS is the single most significant predictor of progression, exceeding the predictive value of any currently used single clinical-pathological feature. We consider this level of predictive power to be clinically useful to support risk-informed decision-making in this patient cohort.
Aneuploidy, and underlying TP53 alteration, have long been recognised as having a central role in the pathogenesis of CA-CRC,10 11 15 23–28 and previous work has emphasised the enrichment of aneuploidy in dysplastic lesions and progressive disease.16 29–32 Our data provide additional evidence of a likely causative role for aneuploidy in CA-CRC. Our whole genome sequencing data indicate that distinct CNAs are enriched in CA-CRC and precursor dysplasias compared with sporadic CRC,11 27 33–36 presumably because the distinct CNAs provide a fitness benefit in the highly inflamed colitic microenvironment.10 For example, losses involving 17p and gains in 7p are common in both sporadic adenomas and IBD-LGD, whereas 4q loss and 8q gain are significantly over-represented in IBD dysplasia, and chr17 losses are enriched in progressive colitis lesions.
We and others have previously established that CA-CRC arises following a process of ‘field cancerisation’ whereby cells carrying cancer-associated mutations expand through the intestinal mucosa without necessarily causing morphological change37–39 (ie, without necessarily causing dysplasia). We suggest that field cancerisation explains why CNAs detected in LGD lesions predict the future risk of CA-CRC development, despite often complete resection of the LGD lesion: we expect that the CNA-bearing LGD lesions arise within a larger clonal field of morphologically unremarkable cells all carrying the same cancer-risk conferring CNAs, which are left in the bowel after LGD resection. We note that the presence of a large clonal field means that it is likely not necessary to carefully enrich for dysplastic cells when performing genomic analysis of LGD lesions, as the non-dysplastic cells at the dysplasia margins are expected to carry the same alterations as the dysplasia itself. In our study, we used an ad hoc mixture of laser microdissection and needle microdissection, with no consequential decrease in genomic signal for the latter coarser method. Our future work will investigate CNA burden as predictor of advanced neoplasia risk in non-dysplastic mucosa.
Our study had several limitations. Our discovery cohort consisted of patients exclusively from St. Mark’s Hospital, a specialist tertiary referral centre in the UK, where the patient cohort is enriched for complex and high-risk cases. Our validation cohort was drawn from other major centres across the UK and so may also be enriched for high-risk cases. Therefore, translation of the CNA biomarker into the general (lower-risk) surveillance population may require further refinement. Our dataset included no information on family history, a recognised (although weak) predictor of CA-CRC; future work to include this factor may further improve risk prediction. We note that LGD lesions were annotated by expert pathologists prior to microdissection for genomic analysis, but as discussed above, we think that enrichment of dysplastic cells will not be necessary for practice because dysplasia tends to arise from a cancerised field bearing the same genetic alterations.
In conclusion, we have derived a novel genomic biomarker of future CRC risk in patients with IBD with LGD. The biomarker uses standard short-read very low-coverage whole genome sequencing data and consequently we expect minimal barriers to subsequent clinical implementation. This biomarker has major translational potential, and we anticipate it will improve the clinical management of IBD, reduce CRC diagnoses and improve quality of life for patients with IBD.
Methods
FFPE human tissue samples, with corresponding anonymised clinical data, were obtained from St Mark’s Hospital from patients with extensive ulcerative colitis (Montreal Classification E3). All pathology was confirmed by two blinded GI pathologists (from MoM, MJ, MR-J and NAW) using the corresponding H&E slide. Where two pathologists were in disagreement on dysplasia grading, the review of a third pathologist was sought, and the majority grading severity was chosen.
Where a patient had multifocal LGD, all dysplastic lesions were taken forward for analysis. In our St. Mark’s cohort (called the discovery dataset throughout), we aimed for 1:2 patient matching, with one progressor (P) for every two NP roughly corresponding to progression rate in LGD. Clinically relevant risk factors such as patient age, gender, duration of IBD diagnosis and the presence of PSC were matched between P and NP. Moreover, endoscopic factors such as lesion size, shape and location were noted. Censor time is defined as the number of days from the procedure at which the LGD tissue sample is obtained until: (1) the date of the colonoscopy or operation at which either HGD or CRC is first detected (for P) or (2) the date of the final colonoscopy confirming that the patient remains HGD and CRC-free (for NP).
For the Newcastle and Oxford validation cohorts, we replicated this as closely as possible, however, these cohorts had only 3 years of follow-up data at the censoring time.
DNA extraction from FFPE tissue
Lesions with low epithelial content were sectioned onto UV-treated, PEN-membrane covered slides (Zeiss) to allow for epithelial enrichment by laser-capture microdissection, or they were sectioned onto glass slides for epithelial enrichment by needle macrodissection. In each case six sequential 10 µm sections were used. Sections were dewaxed in xylene for 5 min and rehydrated using 5 min incubations in dilutions of ethanol (100%, 90% and 70%) before staining with methyl green (Sigma Aldrich). Sections were briefly washed in water. Sections on PEN-membrane slides were dehydrated in 100% ethanol for 2 min. A sequential annotated H&E was used as a guide for epithelial enrichment, either through targeted needle macrodissection or laser capture microdissection with the PALM MicroBeam (Carl Zeiss Microscopy).
DNA was extracted using a modified protocol of the High Pure FFPET DNA Isolation Kit (Roche Life Science). Briefly, the isolated epithelium-enriched tissue samples were digested in 23 µL proteinase K and 60 µL lysis buffer for 16 hours at 56°C, followed by enzyme denaturation using a 1-hour incubation at 90°C. After cooling to room temperature, 67 µL of DNA binding buffer is added and incubated for 10 min at room temperature. 33 µL of isopropranol was then added and the solution passed through DNA-binding spin columns twice by centrifugation. After three washes of the spin columns using the supplied wash buffers, DNA was eluted in 30 µL nuclease-free water with an average yield of about 5 ng.
Library preparation
Between 2 and 3 ng of input DNA was used for library preparation with the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs). Minor modifications to the manufacturer’s protocol included the halving of all recommended reagent volumes for cost-effectiveness, the elimination of DNA fragment size selection step prior to addition of adaptor and optimised dilution of Illumina adaptor (dependent on input DNA). Libraries were barcoded using the NEBNext Multiplex Oligos for Illumina kit (Dual index Primers Set 1) and subjected to 11–13 cycles of library amplification. The final library quantity was confirmed with the Qubit Fluorometer, and quality (fragment size distribution and the presence of unbound adaptor) was verified using the Agilent High Sensitivity D1000 DNA TapeStation system. Libraries were sequenced if the peak fragment size was greater than 184 base pairs. Samples were equimolar pooled and sequenced to a target 0.1× coverage using the Illumina NextSeq 500.
lpWGS: preprocessing bioinformatics pipeline
Bioinformatic analysis was performed on the QMUL High Performance Computing Cluster (Apocrita).40 After merging flow cell FASTQ output for every sample, reads in the generated files are subjected to adaptor trimming using Skewer,41 with trimmed reads less than 35 base pairs long filtered out. A quality control report (including per base sequence quality and estimated duplicate rate) is generated for each FASTQ file using FASTQC.42 The reads are aligned against the reference genome hg19 using the Burrows-Wheeler Aligner package (BWA-MEM)43 to produce corresponding SAM files. These generated files are then converted to compressed BAM files, with filtering out of poorly mapped reads (defined as a MAPQ score of less than 10−37/10, which is approximately less than 1 in 5000) using SAMtools.44 Duplicate reads were subsequently filtered out using Picard Tools,45 and BAM index files were created for more rapid data extraction including read visualisation on Integrated Genome Viewer.46 A quality control report relating to mapping quality is generated for each BAM file using BAMQC,47 and mean genome coverage is estimated using GATK.48
Next, the QDNAseq49 package in R was used to arrange BAM reads into bins that are 500 kilobase pairs (kbp) in size, followed by correction for GC content, the filtering out of ‘blacklisted’ bin clusters identified by ENCODE50 (poorly mapped genomic regions rich in nucleotide repeats such as centromeres, telomeres and microsatellites) and by the QDNAseq package developers and our own team (other recurrently deleted bin clusters seen in both normal and neoplastic tissue, often located in pericentromeric and subtelomeric regions), and then the smoothing of outliers. Sex chromosomes are also excluded from downstream analysis. The remaining 4401 bins are subsequently normalised against the mean followed by log2 transformation and segmented using R package DNAcopy.51 To allow for differences in cellularity and ploidy, we then used R package ACE52 to obtain a best-fit model of the segmented copy numbers to determine the cellularity and ploidy for each sample with the following function call: squaremodel (template, penalty=0.5, penploidy=0.5), where template is the segmented copy numbers output from QDNAseq. Chromosomal CNAs were then identified using CGHcall53 where per sample cellularity was given as input from ACE output, with each bin assigned a probability for the different types of CNA (losses, gains, complete deletions and amplifications of chromosomal regions). Bins with a CNA probability of greater than or equal to 50% probability are ‘called’ as per CGHcall recommendations and taken forward for downstream analysis in R (see figure 3 for representative copy number profiles). CNAs were also derived from sporadic (non-IBD) adenomas as previously described.22
lpWGS: downstream analysis of copy number
The following measurements were made to compare the tissue types of interest:
Genomic aneuploidy: The proportion of the 500kbp bins with a called CNA.
The number of chromosomal CNA events: The number of individual segments where the constituent bins have a called CNA (eg, in figure 2, part b there are three CNA events on chromosome 17).
A chromosomal arm loss or gain of interest: A chromosomal arm gain or loss of interest was considered to be present if a segment(s) within that arm have been annotated to have a gain, loss, deletion or amplification by CGHcall.
Statistical comparisons of these measurements were made in R using Mann-Whitney U tests (for measurements 1–2) or Fisher’s exact tests (measurement 3), with a p value of 0.05 used as a cut-off for significance. Where multiple comparisons are made, a Benjamini-Hochberg correction was used to generate adjusted p values. Violin plot generation was performed in R using the package ggplot2.
MSI classification pipeline using lpWGS
In this study, we employed our in-house developed method MILO (Microsatellite Instability in Low-quality biopsies) to investigate the MSI status in all IBD biopsies.54 MILO is specifically designed for the assessment of MSI status in low-quality samples, including low-coverage FFPE biospecimens. We identified highly sensitive genomic features for MSI classification from publicly available MSI cohorts with low-quality data and subsequently used this information to create MSI classifiers, referred to as MILO. We demonstrated the performance of MILO using previously labelled datasets, achieving an accuracy of 0.99. For more comprehensive information, please refer to the original manuscript.
Whole genome doubling classification
Using International Cancer Genome Consortium (ICGC) segment calls,55 we binned the genome into 4 Mb bins to mimic data expected in low-coverage WGS. Here, we took the median calls for each bin when comparing overlapping segments. For each sample, we calculated the number of segments after binning by detecting changes in call status in adjacent bins across the genome. We calculated PGA by measuring the fraction of bins across the genome that did not equal the median copy number of the bins of the sample (baseline copy number). Taking the called genome-doubled (GD) status per sample from the ICGC resource, we trained a support vector machine using the e1071 package in R (no scaling, a linear kernel and a cost of constraints violation of 5). The classifier applied to these data calls GD as present in a sample (figure 3) if: 0<0.000879NS+PGA-0.347, where NS denotes the number of segments and PGA denotes the percentage genome altered (fraction of genome with non-baseline copy number).22 The classifier shows very good performance when applied to Pan-Cancer Analysis of Whole Genomes data.
Logistic regression modelling and survival analyses
For patients with samples taken from multiple regions and/or at multiple time points, the region with maximal number of chromosomal CNA segments at the first time point was used in the calculation of the number of CNA segments and taken forward for logistic regression and survival analysis. For the case–control study design, univariate and multivariate logistic regression analyses were performed for the discovery data in R using the glm function. The final multivariate model was chosen by performing backward stepwise selection in R using the step function and initially including all UC-CaRE variables and genomic score. For classification of ‘low-risk’ and ‘high-risk’ groups using the multivariate model, the predicted probabilities from function glm.predict in R were used for each patient with complete data and a threshold of 0.50 was used to designate a high-risk patient. To use time-to-event information, patients with LGD were censored at time ‘c’ when HGD and/or CRC was first detected, or at the final recorded clinical report if the patient did not progress to HGD or CRC, as described above. KM survival estimation analyses with right censoring were performed in R using the survival package. HRs and CIs were computed using output from survdiff function, and log-rank tests were used for statistical comparison of the two KM curves. Plots of KM curves were created using the survminer package.
Prediction accuracy
For calculations of NPV and PPV over time, we employed a non-parametric approach that uses the values of the KM estimator directly for predicted high-risk (PPV) and low-risk (NPV) groups (see ref56 for statistical details). For time-dependent calculations of sensitivity, specificity and AUC, we similarly considered a range of censoring times c (typically c was 6 months, 1 year, 2 years, etc). We then relabelled patients as ‘progressors’ or ‘non-progressors’ according to their progression status by censoring time c. If the follow-up for an NP was shorter than the censoring time being considered, then this patient was excluded from the calculation for this time point (as it was possible that the patient had become a progressor, but this was not captured in our data). Confusion matrices for validation data are shown in online supplemental figure S7 which clarify patient numbers at each censoring time point considered. AUC was calculated using the ROCR57 package in R. For each time point, 95% bootstrap CIs were calculated by simulating 1000 runs with ¼ of validation data held-out.
Multiregion analysis and phylogenetic tree reconstruction
To assess intralesion heterogeneity, MEDICC258 was used to create phylogenetic trees for lesions with multiregion samples, where branch lengths are proportional to copy number ‘events’. We transformed the integer relative copy number changes (−2 to +2) from CGHcall into pseudo absolute copy number calls (0 to 4) by normalising so that a copy number change of zero was considered equivalent to diploid state of CN=2. The number of ‘clonal’ CNAs can be visually interpreted as the ‘trunk length’ of the phylogenetic tree, while subclonal CNAs per sample are the average distance of each tip from the trunk.
Supplemental material
Supplemental material
Data availability statement
All data relevant to the study are included in the article or uploaded as supplementary information. Processed data to recreate all modelling results is included in online supplemental file 2 (online supplemental table S2). Copy number segmentations and CNA calls for all 270 samples are available on Mendeley Data (doi:10.17632/vm7wjjb36y.1). Raw sequencing data are not available due to the retrospective and anonymous nature of the study. Code for reproducing main results is publicly available at https://github.com/yosoykit/C4L_CNA.
Ethics statements
Patient consent for publication
Ethics approval
This was a retrospective anonymised analysis of archival patient samples. The study was approved by the UK REC approval number 18/LO/2051.
References
Footnotes
IAB and KC are joint first authors.
SJL, ALH and TAG are joint senior authors.
X @DrPipDunne, @trevoragraham
Contributors Study conception and design: TAG, ALH, SJL, IAB and KC. Bioinformatics: KC, GDC, HEG, SN and QG. Tissue processing and genomics: IAB, NN, KS, HD, TC, CK, MaM and A-MB. Provision of samples: PDD, MBL, AS, JEE, MoM and ALH. Pathology and endoscopy data: JF, ALH, SJL, NAW, MR-J, MJ and MoM. Paper writing: TAG, ALH, A-MB, IAB and KC. TAG is responsible for the overall content (as guarantor).
Funding The study was principally funded by Cancer Research UK (Early Detection project award 25901) to TAG and SJL, and the Barts Charity (large project grant 472-2300) to TAG. We acknowledge a contribution from Bowel Research UK to TAG via their PhD student funding scheme and pilot funding from the 40tude charity to TAG and AH. TAG received additional funding from Cancer Research UK (CRUK A19771 and DRCNPG-May21_100001). SJL was supported by CRUK Program Grant (DRCNPG-Jun22\100002). IAB was funded by an MRC clinical studentship. Support for this research was provided by the NIHR Imperial BRC. NAW received funding from Cancer Research UK Program Grant A2187. MJ is supported by a Cancer Research UK Clinician Scientist Fellowship (A22745). KC received funding from an MRC HDR-UK programme (UKRI Rutherford Fund Fellowship). This work was supported by AGA Research Foundation (AGA Research Scholar Award AGA2022-13-05; KC) and NIH grants (R01 CA270235, P30 CA023100; KC). This work was supported in part by a Merit Review Award I01 BX005958 (KC) from the US Department of Veterans Affairs Biomedical Laboratory Research and Development Service. The study was supported in part by the NIDDK-funded San Diego Digestive Diseases Research Center (P30 DK120515). JEE is funded by the National Institute for Health Research Oxford Biomedical Research Centre. The authors would like to thank Dr Adam Brentnall and Dr Ronghui Xu for helpful statistical discussions.
Disclaimer The contents do not represent the views of the US Department of Veterans Affairs or the United States Government. The views expressed are those of the authors and not necessarily those of the National Health Service, the National Institute for Health Research, or the Department of Health.
Competing interests The authors are in discussions about potential commercialisation and clinical translation of the findings described here. AH has served as consultant, advisory board member or speaker for AbbVie, Arena, Atlantic, Bristol-Myers Squibb, Celgene, Celltrion, Falk, Galapogos, Lilly, Janssen, MSD, Napp Pharmaceuticals, Pfizer, Pharmacosmos, Roche, Shire and Takeda. KC has an investigator-led research grant from Phathom Pharmaceuticals. TAG and A-MB are named as coinventors on a patent application that describe a method for TCR sequencing (GB2305655.9), and TAG is named on a method to measure evolutionary dynamics in cancers using DNA methylation (GB2317139.0). TAG has received honorarium from Genentech and DAiNA therapeutics.
Patient and public involvement Patients and/or the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.
Provenance and peer review Not commissioned; externally peer reviewed.
Author note GDC, HEG, NN and KS are joint second authors.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.