GeneBreak: detection of recurrent DNA copy number aberration-associated chromosomal breakpoints within genes

DNA copy number profiles

The breakpoint detection method we provide is amenable for data from any DNA copy number discovery platform, e.g. array-CGH and (low-pass) WGS. For optimal results, ‘GeneBreak’ takes DNA copy number data that are pre-processed by the R-package ‘CGHcall’¹³ or ‘QDNAseq’¹⁴, both based on the Circular Binary Segmentation algorithm¹⁵, as input. Alternatively, segmented values (log2-ratios) from a different copy number detection algorithm can be used. In addition, it is recommended to provide discrete DNA copy number states (e.g. loss, neutral, gain) that can be used for breakpoint selection. Bioconductor vignette and manual describe commands and workflows in detail (See Supplementary material).

Breakpoint detection and filter options

Breakpoints are defined by the chromosomal locations that separate the contiguous DNA copy number segments pinpointed by a segmentation algorithm. ‘GeneBreak’ identifies chromosomal breakpoint locations for each individual DNA copy number profile. Instead of taking all detected breakpoints, users may want to define more precisely what breakpoints to take into account, based on the two flanking DNA copy number segment characteristics. One of the following three selection options can be applied. A) Copy number-deviation: this selects breakpoints where the shift in log2-ratio between two consecutive DNA copy number segments exceeds the user-defined threshold; B) CNA-associated breakpoints: this selects all breakpoints between consecutive DNA copy number segments, except for breakpoints flanked by two copy number neutral segments; C) CNA-breakpoints: this selects only those breakpoints flanked by segments with dissimilar discrete DNA copy number states.

Breakpoints are defined by the genomic start position of the copy number segments. DNA copy number profiling data is typically granular due to the distance between microarray probes or bin size of WGS copy number data. This means that the genomic location of a breakpoint is not detected at nucleotide resolution but represents a chromosomal interval with a size that is determined by microarray probe density or WGS bin size.

Breakpoint gene identification

For identification of genes affected by chromosomal breakpoints the built-in gene annotations can be used. Alternatively, a user-defined gene annotation file can be provided (see Bioconductor vignette and manual for further details). The implemented mapping approach identifies genes that are associated with one or multiple chromosomal breakpoint intervals.

Cohort-based breakpoint statistics: breakpoint and gene level

Identification of statistically recurrent breakpoint events across all samples can be performed on both chromosomal location- and gene-level. As features, i.e. microarray probes or bins of WGS copy number data, are (nearly) equally distributed over the genome, we assume that the null- probability for breakpoint occurrence is equal for all individual candidate breakpoints (features). It differs per sample, though, and equals p_s = N_s/N, where N is the number of probes, and N_s the total number of breakpoint for samples. The test statistic is T_p is the total number of breakpoints for probe p across all samples. Then, under the null-hypothesis, T_p is simply a sum of independent Bernoulli (p_s) random variables, the null-distribution of which is the same for all probes. It is quickly computed by using probability generating functions, giving also the p-values for any observed value of T_p.

The probe-based statistical analysis uses Benjamini-Hochberg false discovery rate (FDR) correction for multiple testing. For the intended use at gene level, a more advanced statistical null-model is required. For the gene level, the null-probability for a breakpoint to occur within an individual gene, depends on 1) the length of the gene, 2) the number of gene-associated features and 3) the number of breakpoints in the entire tumor profile for the specific sample. Therefore, at gene-level, we apply a linear regression-based correction for covariates. These regression-estimates are then used as gene- and sample-specific breakpoint null-probabilities (p_g,s). The test statistic remains the same, and so does the null-distribution computation, although it has to be repeated for each gene now. Finally, the Gilbert FDR correction that accounts for discreteness in the null-distribution¹⁶ is applied in this analysis to determine significance of recurrent breakpoint genes. Commands and example workflow can be found in Bioconductor vignette and manual.

Identification of recurrent breakpoint genes in advanced colorectal cancers

We applied our method to 352 high-resolution array-CGH samples from a series of advanced colorectal cancers¹⁷ following CNA detection using ‘CGHcall’¹³. Array-CGH data are available in the Gene Expression Omnibus database under accession number GSE63216 (www.ncbi.nlm.nih.gov/projects/geo/). We selected for the CNA-associated breakpoints (setting: ‘CNA-associated’), used gene annotations from ensembl (human genome NCBI build36/hg18, release 54) and applied the dedicated Benjamini-Hochberg-type FDR correction (setting: ‘Gilbert’), for recurrent breakpoint gene identification. A total of 748 genes appeared to be recurrently affected by chromosomal breaks (FDR<0.1)⁵. Breakpoint frequencies of chromosome 20p are visualized with the built-in plot function (Figure 1; see Bioconductor vignette and manual for further details about this function). Interestingly, patient stratification based on recurrent gene breakpoints and well-known point mutations by propagation to the predefined STRING human protein interaction network revealed one CRC subtype with very poor prognosis, which supported clinical relevance of this class of somatic aberrations in advanced colorectal cancers⁵.