Detection of Copy Number Alterations by Shallow Whole-Genome Sequencing of Formalin-Fixed, Paraffin-Embedded Tumor Tissue

Patients qualified for inclusion if both FFPE (nondecalcified; fixation time between 6 and 48 hours) and fresh material were available for sequencing. Notwithstanding that selection was random, we aimed at combining a variety of disease entities, including tumor types routinely tested by clinical laboratories. In total, 21 patients were included, whereof 17 had available aCGH (60K array, Agilent Technologies, Santa Clara, California) profiles (generated in accordance with Kumps and colleagues¹⁴ ) (Supplemental Tables 1 and 2; see supplemental digital content at www.archivesofpathology.org in the August 2020 table of contents, containing 3 tables and 18 figures). Nine more nontumoral FFPE samples, assumed to contain no CNAs, were sequenced as negative controls. These samples comprised lymph node, colon, appendix, ureter, tonsil, skin, stomach, lung, and brain tissue.

This study was conducted according to the guidelines of the institutional ethics committee at Ghent University Hospital (Ghent, Belgium), project number 2004/094.

Although DNA from fresh material was directly available for shearing through routine clinical practice, 7 sections were taken from the paraffin blocks to initiate FFPE DNA extraction. The first and last section (1.5–2 μm) were stained with hematoxylin-eosin and served as references to guide dissection and to determine tumor cell percentage, whereas the others sections (10 μm) were dissected and deparaffinized by deparaffinization solution (Qiagen, Hilden, Germany). Extraction was executed using the QIAamp DNA FFPE tissue kit (Qiagen), according to the manufacturer's instructions, followed by quantification using a Qubit fluorometer and the Qubit dsDNA BR assay kit (Life Technologies, Bleiswijk, The Netherlands).

Size distribution of the FFPE-derived DNA was determined using the Fragment Analyzer, the High Sensitivity Small DNA Fragment Analysis Kit, and PROSize data analysis software, according to the manufacturer's instructions (Agilent).

DNA shearing to 200–base pair (bp) fragments was completed by Covaris' adaptive focused acoustics shearing process and an M220 focused ultrasonicator (Covaris, Woburn, Massachusetts). Different shearing treatment times were evaluated. Using 200 ng/50 μL DNA as input, the following treatment parameters were set: peak incident power = 50 W; duty factor = 20%; cycles per burst = 200; treatment time = 160 seconds for fresh DNA, 100 seconds for FFPE-derived DNA; and temperature = 20°C. The fragments' size distributions were evaluated using the Fragment Analyzer, the dsDNA 910 reagent kit, and PROSize data analysis software (Agilent).

DNA libraries were prepared using the NEXTflex rapid DNA-seq kit (following option 2, with bead-based size selection prior to polymerase chain reaction amplification) and NEXTflex DNA bar codes (PerkinElmer, Austin, Texas), starting from 128 ng of DNA. Libraries were examined using the Fragment Analyzer, the dsDNA 910 reagent kit (fresh DNA) or the High Sensitivity Small DNA Fragment Analysis Kit (FFPE DNA), and PROSize data analysis software, and by quantitative PCR using the KAPA library quantification kit (Roche, Pleasanton, California).

Libraries were finally equimolarly pooled (4–6 nM). Cluster generation was completed using a cBot 2 system (Illumina, San Diego, California). Sequencing was performed by Illumina's HiSeq 3000 (Illumina) in a single-read 50-cycle run mode. The minimal number of reads per sample was set at 15 million (mean coverage of 0.25x).

Demultiplexing was executed by the bcl2fastq software (v2.19.1, Illumina). To allow fair comparison, paired fresh and FFPE sample data were read normalized by randomly removing reads from the highest coverage sample until equal read depths were obtained. Raw reads were mapped by Bowtie 2 (v2.2.8) onto human reference genome GRCh37, using the fast-local flag,¹⁵  where the resulting BAM files were indexed by SAMtools (v1.5).¹⁶  Reads with a mapping quality lower than 40 were removed. To deduce normalized binwise log₂ ratios between the observed and expected copy number, the reference-free QDNAseq (v1.18.0) was favored.⁵  A bin size of 50 kilobases (kb) was selected, resembling the resolution of the aCGH 60K arrays (Agilent). Profile segmentation was executed by circular binary segmentation, implemented in the R package DNAcopy (v1.56.0).¹⁷  At some loci, copy number is undeterminable by short-read shallow whole-genome sequencing (eg, at centromeres): QDNAseq masks these regions by applying a blacklist. Therefore, a unified blacklist, defined as the union of all obtained blacklists, was applied to all resulting profiles (including aCGH), enabling unbiased concordance analyses.

Assuming our negative control DNA had no structural aberrations, the variance across these flat FFPE-derived profiles illustrated the minimal required tumor fraction (TF) to significantly call aberrations. These profiles were sequentially divided into segments of length 100 kb (2 bins) until 10 Mb (200 bins), where their log₂ ratio equaled the interpolated average of its bins'. For each segmental locus, the standard deviation (sd) and mean (x) of the log₂ ratio across nontumoral profiles was derived—this for every evaluated segment length. Because genomic loci differ in their level of normal variance, the median segmental locus was used to calculate the required log₂ ratio (r) to reach significance (z > 3.29 for gains, z < −3.29 for losses, corresponding to the Bonferroni corrected 5% significance level, assuming 100 segments), using r = ± (3.29 · sd − x). This value corresponds to duplications (3n) and heterozygous deletions (1n) for samples with a TF of ±(2^r · 2) − 2. These formulas allow us, given an aberration size lower detection limit (1 Mb in this study), to calculate the minimal required TF.

The TF was estimated by ichorCNA (Broad Institute, Cambridge, Massachusetts), a software that predicts the portion of tumor fragments based on observed structural aberrations.¹⁸  Although ichorCNA is benchmarked on cell-free DNA, the authors ensured its applicability to FF and FFPE material (Viktor Adalsteinsson, PhD, written communication, December 21, 2017). The tool generates elementary copy number profiles (bin size 500 kb), which were normalized using the provided guanine and cytosine content and mappability tracks. Five different starting values, representing the nontumoral contamination, were passed to the optimization algorithm: 10%, 25%, 50%, 75%. and 90%. The obtained tumor content corresponding to the highest likelihood was interpreted as the final TF.

Aberrations were called once the log₂ ratio of profile segments surpassed a gain or loss cutoff valued as log₂\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\bf{\alpha}}\)\(\def\bupbeta{\bf{\beta}}\)\(\def\bupgamma{\bf{\gamma}}\)\(\def\bupdelta{\bf{\delta}}\)\(\def\bupvarepsilon{\bf{\varepsilon}}\)\(\def\bupzeta{\bf{\zeta}}\)\(\def\bupeta{\bf{\eta}}\)\(\def\buptheta{\bf{\theta}}\)\(\def\bupiota{\bf{\iota}}\)\(\def\bupkappa{\bf{\kappa}}\)\(\def\buplambda{\bf{\lambda}}\)\(\def\bupmu{\bf{\mu}}\)\(\def\bupnu{\bf{\nu}}\)\(\def\bupxi{\bf{\xi}}\)\(\def\bupomicron{\bf{\micron}}\)\(\def\buppi{\bf{\pi}}\)\(\def\buprho{\bf{\rho}}\)\(\def\bupsigma{\bf{\sigma}}\)\(\def\buptau{\bf{\tau}}\)\(\def\bupupsilon{\bf{\upsilon}}\)\(\def\bupphi{\bf{\phi}}\)\(\def\bupchi{\bf{\chi}}\)\(\def\buppsy{\bf{\psy}}\)\(\def\bupomega{\bf{\omega}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\left( {{{2 \pm \left( {tumor \ fraction/2} \right)} \over 2}} \right)\)⁠. If the TF did not exceed the previously established lower detection limit, no copy number states were called.

All comparative analyses were executed by custom scripting in R (v3.5.1). The median segment variance (MSV), a measure that reflects the overall vertical spread of log₂ ratios within segments, was adopted to perform profile noise comparison.¹⁹ 

Profile correlation was visualized by hierarchical clustering, using the unweighted pair group method with arithmetic mean as an agglomerative clustering procedure and the Pearson correlation as a concordance measure. The latter was applied to a predefined set of 1000 equally sized windows, obtained by merging segment-interpolated binwise log₂ ratios—this to reduce computational complexness and to increase robustness.

Finally, performance was established in the form of a receiver operating characteristic analysis. The CNVseq/FFPE combination was compared with the CNVseq/fresh and aCGH/fresh setups, using aCGH/fresh and CNVseq/fresh as gold standards, respectively. Patients without aCGH test or with insufficient TF (according to the lower detection limit) were excluded. Specificity and sensitivity were measured while sliding the aberration cutoff from liberal (ratio of 1) to conservative (0.5 for deletions; 1.5 for gains). The mean area under the curve defined the overall performance.

Extraction of FFPE DNA yielded on average 4.56 ± 0.71 μg DNA (Supplemental Table 2). Fragment analysis of this DNA showed a wide size distribution with both highly degraded and less-degraded fragments (Supplemental Figures 1 through 3). Four FFPE samples needed no further fragmentation (cases 5, 12, 14, and 17). After DNA shearing, fragmented FFPE-derived DNA showed an overall wider size distribution when compared with fragmented fresh tumor tissue DNA (Supplemental Figures 4 through 7). The DNA quality number, a quality metric provided by the PROSize data analysis software, gives an indication of the fraction of fragments with sizes above a specified threshold, set at 200 bp. The average DNA quality number was 5.9 ± 0.1 (Supplemental Table 3) for fresh DNA and 6.4 ± 0.1 (Supplemental Table 3) for FFPE-derived DNA, supporting different size distributions (paired t test; P < .001).

Library construction resulted in on average 39.90 ± 4.79–nM libraries, starting from fresh samples (Supplemental Table 1), and 9.18 ± 2.32–nM libraries starting from FFPE samples (Supplemental Table 2). Library validation showed different traces for FFPE-derived libraries compared with fresh DNA libraries. Fresh samples showed a peak around the expected library size of 300 bp, but had a right-skewed distribution. Formalin-fixed, paraffin-embedded samples displayed 2 peaks, one around 300 bp and one of higher molecular weight (Supplemental Figures 8 through 12).

Very similar copy number profiles were obtained across paired fresh and FFPE samples (Supplemental Figures 13 through 18), as illustrated, for example, by case 14 (Figure 1 and Supplemental Figure 16, n). Here, copy number predictions correspond to each other, with the exception of chromosomes 5 and 11, presumably caused by tumor heterogeneity.

Measuring copy number profile noise for CNVseq, FFPE material showed more variance in comparison with fresh material, measured by means of MSV (paired t test; P = .001; Figure 2, a). The mean MSV for aCGH was 0.025 ± 0.003, where for CNVseq this measure amounted to 0.019 ± 0.001 for fresh material and 0.026 ± 0.002 for FFPE-derived DNA (Supplemental Table 3). For 2 cases, the aCGH profile (eg, case 16; arrow in Figure 2, a; Supplemental Figure 16, p) displayed an abundance of false outlier bins.

The increased noise levels observed in FFPE were partly caused by the number of reads that were available for CNVseq (Figure 2, b). The elevated variance levels in FFPE-derived copy number profiles were partly caused by the number of available reads for CNVseq, which is expected to be inversely proportional to noise: although both FFPE and fresh sample read files were sampled towards the same number of raw reads, mapping and quality filtering resulted in fewer final reads for FFPE material (paired t test; P = .003). The percentage of mapped reads was approximately 70% for fresh samples (13.7 million reads mapped of 19.4 million raw reads; Supplemental Table 1) and 65% for FFPE samples (12.4 million reads mapped of 19.3 million raw reads; Supplemental Table 2), resulting in a mean coverage of 0.21x ± 0.07x for fresh DNA and 0.19x ± 0.06x for FFPE-derived DNA. As expected, there was an inverse relation between the MSV and number of mapped reads (Figure 2, c): MSV levels were higher for FFPE samples, with the exception of 3 cases (cases 1, 18, and 19; Supplemental Table 3). These observations were not related to DNA quality number (r = 0.11, P = .54). The use of FFPE tissue for CNA detection by CNVseq would theoretically require more sequencing reads compared with fresh material when aiming at obtaining equal noise levels (Figure 2, d).

Concordance was supported by unsupervised clustering applied to the fresh and FFPE-derived CNVseq profiles, which resulted in an almost perfectly paired pattern (Figure 3). Matched samples clustered together for all but 2 cases, which showed no detectable CNAs (cases 7 and 9; Supplemental Figures 14, g, and 15, i).

In this study, we aimed at detecting CNAs of at least 1 Mb.^3,20,21  Assuming a coverage of 0.25x (the negative control samples had a mean coverage of 0.25x ± 0.01x, where coverage impacts noise) and 100 segments (the case samples had on average 96 ± 13 segments, where the number of segments impacts Bonferroni correction), the minimal required TF to detect CNAs of at least 1 Mb was measured to be 20% (Figure 4, a). Three cases (cases 4, 7, and 9) had a TF below the lower detection limit (Supplemental Tables 1 and 2). This naturally explains why matched samples do not cluster for cases 7 and 9 (Figure 3).

Although the variance in the profiles of FFPE material was often higher, this had little impact on the performance of aberration calling (Figure 4, b and c). Receiver operating characteristic analysis using aCGH as a gold standard showed a mean area under the curve of 0.960 for CNVseq of fresh material and 0.949 for CNVseq of FFPE material (Figure 4, b). For one case (case 17), CNVseq showed poor concordance with aCGH. The aCGH profile (Supplemental Figure 17, q) displayed an abundance of, probably, false-positive bins.

Using CNVseq in combination with fresh material as the gold standard, the mean area under the curve for FFPE material amounted to 0.973 (Figure 4, c).

In one case (case 12), a 6.2-Mb gain on the short arm of chromosome 2 (amp(2)(2p16.1p14)) could be detected in both fresh and FFPE-derived DNA samples, using aCGH or/and CNVseq (Supplemental Figure 15, l), illustrating that small CNAs could be detected using FFPE tumor tissue. In this region, the REL gene is located; amplification of this gene is a recurrent event in B-cell lymphomas.

In this study, we demonstrate the feasibility of CNA detection by CNVseq using FFPE-derived DNA. In total, 21 patients were selected during a recent time frame (2016 for cases with aCGH data)—this to minimize the effects of fixation and storage time.^6,10  Nondecalcified FFPE material was chosen as decalcification introduces additional DNA damage.²²  We observed that the DNA integrity of FFPE tissue was less affected than expected. Like Hostetter et al⁷  and Munchel et al,¹²  but different from Mc Sherry et al,¹  we saw a wide size distribution of FFPE-derived DNA, indicating only partial DNA degradation. Higher-molecular-weight DNA was more present than expected.

Very similar CNA profiles were obtained across paired fresh and FFPE samples. Differences possibly were caused by sampling errors or resulted from tumor heterogeneity, as previously cited by Little et al.²³ 

Although we observed less degradation, poorer DNA quality from FFPE samples was obvious. The higher variance seen in our FFPE-derived copy number profiles is partly caused by the lower-quality FFPE-derived DNA or the higher error load in FFPE-derived sequencing data. This was reflected mostly during the mapping phase of sequencing. In contrast to Schweiger et al,¹⁰  who saw variations in read mappability independent of the DNA source, we observed that FFPE samples had a lower percentage of mapped reads; thus, the lower percentage of mapped reads was probably due to DNA quality.

Using nontumoral FFPE samples, we determined that the minimum TF necessary for accurate aberrations calling on FFPE material is 20%. This is similar to the sensitivity of aCGH, debated to be 20% as well,⁷  whereas with targeted sequencing, where coverage is regionally increased, a lower detection limit is possible.²⁴ 

Besides detection of large chromosomal alterations, the determination of single-gene copy number status is also important in cancer diagnosis, prognosis, and treatment. We showed that detection of even small CNAs, such as the gain at 2p (containing the REL gene) that we could readily detect in our study, is possible using CNVseq on FFPE material. Classical techniques such as conventional karyotyping are limited to a resolution of 5 to 10 Mb because of chromosome banding.²⁵ 

We showed that the higher variance seen in FFPE-derived profiles has little impact on the performance of statistically calling an aberration if TF is at least 20%.

Although several studies have reported massive parallel sequencing applied on FFPE samples, especially in the context of mutation detection, CNA detection on FFPE material is not yet commonly executed. Although our study is limited to recently archived samples, it shows the possibility of performing CNA detection on routinely collected FFPE material. The use of FFPE material for DNA analysis has several important advantages. Formalin fixation and paraffin embedding is the standard procedure for examination of surgical biopsies, as it provides optimal morphology. Handling and storage of FFPE material is much easier and cheaper than that of FF or unfixed material. In addition, selection of tissue blocks or dissection of areas with a maximal tumor cell percentage makes the process of tumor tissue selection for DNA analysis much more controlled.

Compared with aCGH, CNVseq can be easily implemented as routine clinical practice, because a very similar workflow, as for noninvasive prenatal testing, is often already implemented in numerous genetic laboratories. Moreover, pooling enables CNVseq for multiple types of analyses at once (eg, in combination with noninvasive prenatal testing), resulting in overall turnaround times of less than 10 days. In addition, using low-coverage sequencing, this test can be performed at a low cost relative to aCGH. When aiming at detecting smaller insertions or deletions, the sequencing depth can be easily increased, but will result in an increase of the price tag.

In contrast to targeted sequencing, shallow whole-genome sequencing cannot be used for single-nucleotide variant detection, as it is low coverage and does not need predefined targets.

We believe that CNVseq for CNA detection is likely to become important for diagnosis, prognosis, and therapy response prediction for a wide range of different solid tumor entities. Copy number variation sequencing can be used for detection of losses or gains at single loci, like HER2 amplification in breast cancer,²⁶  but is especially interesting for tumor types with complex CNA profiles. For instance, bladder carcinoma and melanoma harbor tumor-specific CNAs.^27–29  The current standard diagnostic procedure to detect these CNAs is via fluorescent in situ hybridization (FISH) using different probes. However, this technique is labor-intensive and highly dependent on the quality of the samples. The advantages of CNVseq over FISH are that it is less observer dependent and it is cheaper when multiple probes are needed. Even more important, CNVseq is genome wide, and thus can be used without prior knowledge of the target region and can be applied on every tumor type. A time-consuming FISH test that is currently offered by many molecular pathology centers detects for 1p/19q codeletion in gliomas. Codeletion of 1p and 19q is characteristic for oligodendroglioma, a glioma subtype associated with better prognosis and response to chemotherapy.³⁰  Testing with FISH for oligodendrogliomas requires evaluation of 4 signals in many tumor cells (at least 50–100). Moreover, interpretation of this test is difficult, as exemplified by the many different cutoffs and varying guidelines that are used by different laboratories and in clinical studies.^31,32  Given the advantages of CNVseq, FISH may, in time, be replaced by CNVseq on FFPE tumor biopsies. Copy number variation sequencing on FFPE material has also great potential for the diagnosis of high-grade serous ovarian carcinoma (a tumor with an exceptionally high CNA burden),³³  the evaluation of BRCAness in breast carcinoma, the detection of CNAs in neuroblastoma,³⁴  and the study of CNAs in many other benign and malignant tumor entities.

Formalin-fixed, paraffin-embedded tissue enables harvesting of tumor material in a more controlled manner. The higher variance seen in FFPE-derived copy number profiles has little impact on the performance of statistically calling an aberration if the TF exceeds the established lower detection limit of 20%. Importantly, archival FFPE tissue is available for many tumor entities in laboratories for surgical pathology and thus represents a large biobank for future clinical research.

We are thankful for the technical assistance of Lynn Supply, Tamara De Clercq, and the Fragment Analyzer and CNVseq team.