The International Cancer Genome Consortium (ICGC) Data Access Committee granted us access to the CLL sequencing data (Ramsay et al., 2013) deposited in the European Genome-Phenome Database (EGA) under DACO-1040945. For this analysis, we used exome-seq data from matched control and tumor samples of CLL cases under the accession code EGAD00001001464. Patient characteristics can be consulted in Table 1 .

Samples were processed by Puente et al. (2015) as described in their original paper (Puente et al., 2015). Briefly, 3 μg of genomic DNA was used for paired-end sequencing library construction, followed by enrichment in exomic sequences using the SureSelect Human All Exon 50Mb kit (Agilent Technologies). Next, DNA was pulled down using magnetic beads with streptavidin, followed by 18 cycles of amplification. Sequencing was performed on an Illumina GAIIx or on a HiSeq2000 sequencer (2 × 76 bp). Exome-seq data were aligned to the reference genome (GRCh37.75) using bwa (Li and Durbin, 2009). Duplicate read removal, sorting, and indexing were done using samtools (Li et al., 2009). Base quality score recalibration was made with BamUtil (Breese and Liu, 2013) using a logistic regression model.

We analyzed paired tumor-normal exome-sequencing data with Control-FREEC version 11.3 in order to identify somatic CNA and CNN-LOH regions (Boeva et al., 2012). Control-FREEC uses aligned reads to construct and normalize a copy number profile and a B-allele frequency (BAF) profile. Then, it performs profile segmentation and infers genotype status for each segment using both copy number and allelic frequency information. Finally, genomic aberrations are identified and annotated.

The following specifications were used: “window = 0,” “ploidy = 2,” “breakPointThreshold = 1.2,” “noisyData = TRUE,” “readCountThreshold = 50,” “forceGContentNormalization = 1,” “contaminationAdjustment = TRUE,” “telocentromeric = 50000,” and “mateOrientation = 0.” BAFs were estimated using all variants reported in dbSNP (version 150) with a minimal coverage per variant position of 5 reads and a minimal sequencing quality per position of 20 Phred. Variant calling was limited to regions covered in the SureSelect Exome Capture 50Mb version 4 kit.

CNN-LOH were called with p-values < 0.05 (Kolmogorov–Smirnov test) and an uncertainty upper threshold of 5%. Regions with lowmappability according to UCSC 75-bp mappability tracks (score below 0.5) were filtered out. As expected, we observed regions that seemed prone to CNA erroneous detection. Thus, we decided to apply a hard filter and discard those CNAs significantly enriched in both amplifications and deletions, as well as those located near a telomeric or centromeric region. CNA events were detected using GISTIC2.0 (Beroukhim et al., 2007). GITIC2.0 was run with default parameters plus the arm peel-off filter. Focal recurrent CNAs were defined as those spanning less than 50% of a chromosome arm with a residual q-value < 0.05 and a wide peak size below 10 megabases. A 1 − log2 tumor/normal ratio above 0.3 was used to define amplifications, and a ratio below −0.3 was used to define deletions.

Variables associated with time to treatment (TTT) and overall survival (OS) were analyzed using Cox regression as implemented in the survival R package (Therneau and Grambsch, 2000; R Development Core Team, 2011; Therneau, 2015). Assessment of the proportional hazards assumption was performed using the cox.zph function. We created two different models: a univariate model that only includes CNA status for each gene and an adjusted model that included variables associated with clinical outcome at a p-value < 0.2 in a univariate model. The combined model (CM) for TTT analysis included as covariates donor sex, stage at diagnosis (monoclonal B-cell lymphocytosis (MBL), Binet A, Binet B, and Binet C), and IGHV mutational status, while the CM for OS analysis included stage at diagnosis, IGHV mutation status, and age at diagnosis. The Benjamini–Hochberg (BH) method was used to adjust for multiple testing.

Two hundred twenty patients had matched RNAseq data of CLL-purified cells (accession IDs EGAD00001001443 and EGAD00001000258). Illumina adapters were removed using cutadapt (Martin, 2011), and alignment to the human reference genome (GRCh37) was performed using Hisat2 (Kim et al., 2015) with default specifications. We used the Hisat2-provided Hierarchical Graph FM index for GRCh37 with SNP and Ensembl transcript information. Bam files were sorted and indexed using samtools (Li et al., 2009). Bam files were processed in R (R Development Core Team, 2011) according to the RNAseq gene expression protocol developed by Love et al. (2015). Briefly, bam files were read using Rsamtools (Morgan et al., 2017), followed by gene-level expression estimation using the SummarizeOverlaps function from the GenomicAlignments package (Lawrence et al., 2013). Gene models in GTF format were downloaded from Ensembl (GRCh37.75 version) (Yates et al., 2016). A log2-transformation on normalized frames per kilobase counts was performed. Focal CNAs were classified according to Gistic into low-range events (tumor/normal log2 ratio > 0.3 and <0.9 for amplifications and less than −0.9 and more than −0.3 for deletions) and higher-range events (tumor/normal log2 ratio > 0.9 for amplifications and less than −0.9 for deletions). Correlation between CNA status and gene expression was performed using Spearman’s correlation. A minimum of 5 CNA events with matched transcriptomic data was set for analysis. Furthermore, immunoglobulin and T-cell receptor gene rearrangements were not included. p-values were adjusted for multiple testing using the BH method.

We identified 54 recurrent focal CNAs in the CLL genome (residual q-value < 0.05, Supplemental Figure 1 , Table 2 ). Among them, there were 31 recurrent amplifications with a wide peak size of 75.1 kb ( Figure 1 ) and 23 recurrent deletions with a median wide peak size of 405 kb ( Figure 2 ).

The most frequently amplified regions were found in 11p13 (adjacent to the WT1 locus), 8q24.13 (FAM91A1), 3q25.1 (PFN2 and RNF13), 18q22.1 (TMX3), 3q29 (ACAP2), 17q11.2 (SUZ12), 14q32.33 (adjacent to the AKT1 locus), 13q13.3 (NBEA/BCL8B), 1q25.2 (RFWD2), 5q35.3 (FLT4), 19p13.3 (adjacent to the APC2 locus), 1p22.2 (HFM1), and 1p31.1. Similarly, the most recurrent focal deletions were detected in 13q14.2 (DLEU1 locus), 11q22.3 (ATM locus), 14q21.1 (MIA2), 11q14.3 (NOX4), 14q32.33 (IGH locus), and 4q13.2 (UGT2B15 and UGT2B17 loci). Moreover, we observed frequent deletions of immunoglobulin loci and recurrent amplified regions at TCR genes, likely reflecting deletions present in the T lymphocytes within control samples.

Three focal deletions were associated with TTT (BH q-value < 0.05): 11q22.3 (ATM locus), 14q32.33 (IGH locus), and 7q21.2 (CDK6 locus). 11q22.3 loss was associated with shorter survival too. No event was associated with TTT or OS after adjusting for IGHV status, sex, and disease stage (BH q-value < 0.05). No recurrent gain was associated with treatment-free survival or OS. The association results can be consulted in Supplementary Table 1 .

Furthermore, SETD2 deletions (nine cases) and IRF4 gains (five cases) were nearly significant (GISTIC residual q-values of 0.08 and 0.06, respectively) and associated with clinical evolution. SETD2 deletion was associated with shorter treatment free survival (p-value 1.3 × 10⁻⁸) independent of IGHV status, sex, and disease stage (p-value 2.3 × 10⁻³); and it was also associated with shorter survival (p-value 9.9 × 10⁻³) but not independently of IGHV status (p-value 0.24). Amplifications in IRF4 were associated with shorter survival (p-value 6.4 × 10⁻⁵) in a partially IGHV-independent manner (p-value 0.025). Nevertheless, this finding must be interpreted with caution due to the position of IRF4 near the telomeric end of the short arm of chromosome 6.

We detected significant correlations between some recurrent CNAs and their correspondingly encoded genes ( Supplementary Table 2 ). As expected, deletions in 11q22.3 and 13q14.2 and the expression of their respective genes (ATM, FDX1, and RDX in the first case, and DLEU1, DLEU2, and KCNRG in the second case). Deletions in 6q21 were correlated with lower expression of CDK19, and so did those in 14q21.1 with the expression of MIA2. Furthermore, we detected significant correlations between amplifications in 3q25.1 and 3q29 and expression of their target genes PFN2 and ACAP2, respectively. Surprisingly, an inverse correlation was observed between 11q14.3 deletions and FOLH1B expression.

Similarly, we detected significant correlations between 19 of these CNAs and the expression of 926 protein-coding genes (q-value < 0.01; Supplementary Table 3 ). The CNAs with more correlated genes were deletions in 13q14.2 (389 genes), 12q24.33 (ULK1 locus, 170 genes), 6q21 (CDK19 locus, 79 genes), 11q22.3 (ATM locus, 39 genes), and 3q25.33 (IL12A locus, 26 genes), as well as amplifications in 18q22.1 involving TMX3 (135 genes).

Nonetheless, this study is probably underpowered to detect CNA–transcript correlations due to the fact that only 50% of the exome-seq samples had matched RNAseq data of purified CLL cells.

Recurrent broad amplifications and deletions were detected in the CLL genome ( Table 3 , Supplementary Figure 1 ). Among the amplifications, the most frequent were found in chromosomes 12 (60 cases, 13.6% of patients), 2p (10 patients), 8q (6 patients), 18q (6 patients), 18p (5 patients), and 3q (5 patients). Similarly, the most recurrent large deletions were detected in 17p (8 patients), 18p (7 patients), and 8p (6 patients).

Amplifications in chromosomes 12, 2p, and 8q, as well as deletion of 17p, were significantly associated with shorter TTT ( Supplementary Tables 4 and 5 ). Furthermore, deletion of 17p was significantly associated with TTT independently of IGHV status, sex, and Binet staging. Amplifications in chromosomes 12 and 8q were associated with shorter OS, but no event was significant after adjusting for IGHV status, age, and Binet staging. Nonetheless, we detected a significant difference in OS and TTT between IGHV-mutated cases with and with and without trisomy 12 ( Figures 3A, B respectively).

Control-FREEC identified 63 regions of CNN-LOH with a genotype uncertainty below 5% ( Supplementary Table 6 ). Ten events were located at 1p, six of which affected ARID1A. By comparing with mutation data published by Puente et al. (2015), none of these samples bore concurrent non-synonymous mutations in ARID1A. Other three events were detected at 1q, with a minimally affected region on 1q21.3, which holds likely driver genes such as PI4KB and IL6R/CD126. Four CNN-LOH events with a minimally involved region in 9q34.13–q34.3 involved the NOTCH1 gene. One of these cases also had a frameshift deletion in NOTCH1. Three events at 11p15.5–15.4 affected the imprinted locus of IGF2 and CDKN1C. Eight CNN-LOH events affected the 11p11.2–q13.2 region. Two different samples had events at 11q involving the ATM locus, both of which also had non-synonymous ATM mutations. Three events were located in 16p13.3–p13.11, which involve the CREBBP gene. None of these had mutations in CREBBP. Ten CNN-LOH events affected chromosome 17, three of which included the TP53 locus. Among the latter, two had non-synonymous mutations in TP53. Four samples had CNN-LOH events at chromosome 20, all of which affect the ASXL1 gene but without any concurrent mutation in it. It is interesting to mention that only two of the CNN-LOH events overlapped with those reported by Puente et al. (2015) using SNP array technology.