Reconstruction of clone- and haplotype-specific cancer genome karyotypes from bulk tumor samples

Abstract

Many cancer genomes are extensively rearranged with aberrant chromosomal karyotypes. Deriving these karyotypes from high-throughput DNA sequencing of bulk tumor samples is complicated because most tumors are a heterogeneous mixture of normal cells and subpopulations of cancer cells, or clones, that harbor distinct somatic mutations. We introduce a new algorithm, Reconstructing Cancer Karyotypes (RCK), to reconstruct haplotype-specific karyotypes of one or more rearranged cancer genomes from DNA sequencing data from a bulk tumor sample. RCK leverages evolutionary constraints on the somatic mutational process in cancer to reduce ambiguity in the deconvolution of admixed sequencing data into multiple haplotype-specific cancer karyotypes. RCK models mixtures containing an arbitrary number of derived genomes and allows the incorporation of information both from short-read and long-read DNA sequencing technologies. We compare RCK to existing approaches on 17 primary and metastatic prostate cancer samples. We find that RCK infers cancer karyotypes that better explain the DNA sequencing data and conform to a reasonable evolutionary model. RCK's reconstructions of clone- and haplotype-specific karyotypes will aid further studies of the role of intra-tumor heterogeneity in cancer development and response to treatment. RCK is freely available as open source software.

Figure 1.

Overview of the RCK algorithm. The inputs to RCK (white dotted boxes) are clone- and allele-specific copy numbers (top left) and novel adjacencies (top right) from bulk tumor samples that are derived from alignments of DNA sequencing (top) reads using existing tools. The RCK algorithm (blue shaded elements) builds a diploid interval adjacency graph integrating copy number and novel adjacency information (for details, see Methods). RCK then solves a mixed-integer linear program (MILP) to find an optimal assignment of segment copy numbers and novel adjacencies to alleles and clones, subject to copy number balance on segment ends and satisfying evolutionary constraints from a generalized infinite sites model. Constraints on groups of novel adjacencies from the third-generation sequencing technologies may optionally be included. The outputs of RCK are clone- and haplotype-specific cancer genome karyotypes.

Figure 2.

Results of RCK on simulated bulk tumor samples with two clones. (A) False negative rate (FNR) and false positive rate (FPR) of novel adjacencies used by RCK using adjacency utilization parameter P = 0.9 (RCK-0.9) and P = 0.75 (RCK-0.75). (B) Length-weighted segment copy number distances between input copy numbers (${\widetilde{\overline{\mathbf{C}}}}^{\mathrm{\prime }}$) and karyotypes inferred by RCK.

Figure 3.

Comparison of RCK and ReMixT on heterogeneous prostate cancer samples. (A) Length-weighted segment copy number distances between segment copy numbers from HATCHet and segment copy numbers output by ReMixT and RCK. (B) Fractions of novel adjacencies from input that are inferred to be present by ReMixT or RCK for each sample in the heterogeneous group. RCK used segment copy numbers from HATCHet in input and novel adjacency utilization parameter P = 1.0, 0.9, 0.75, 0.5.

Figure 4.

ReMixT karyotypes from heterogeneous prostate cancer samples have numerous violations of the generalized infinite sites constraints. In A, C, and E, solid edges represent segment edges, black-dashed edges represent reference adjacency edges, and red dashed edges represent novel adjacency edges. Integer values indicate copy numbers of corresponding segment and adjacency edges. (A) An intra-genome violation of the homologous-extremity-exclusivity constraint. To achieve copy number balance, both homologous vertices $2_{\mathrm{A}}^h$ and $2_{\mathrm{B}}^h$ from genome G_i must be involved in novel adjacencies. (B) Number of novel adjacencies that violate the intra-genome homologous-extremity-exclusivity constraint in each cancer karyotype inferred by ReMixT in each sample. (C) An inter-genome violation of the homologous-extremity-exclusivity constraint. To achieve copy number balance, both homologous vertices $2_{\mathrm{A}}^h$ and $2_{\mathrm{B}}^h$ (in different genomes) must be involved in novel adjacencies. (D) The fraction x/y, where x is the number of novel adjacencies that violate the inter-genome homologous-extremity-exclusivity constraint (on at least one of the extremities involved in a novel adjacency) in ReMixT karyotypes, and y is the total number of novel adjacencies reported by ReMixT as being present in both genomes. (E) A violation of the intra-genome homologous-reciprocal-extremity-exclusivity constraint. To achieve copy number balance, both homologous-reciprocal vertices $2_{\mathrm{A}}^h$ and $3_{\mathrm{B}}^t$ must be involved in novel adjacencies. Inter-genome violations of the homologous-reciprocal-extremity-exclusivity constraint are also possible (Supplemental Fig. S17). (F) Fraction x/y, where x is the number of reciprocal locations with violations of either intra- or inter-genome (or both) homologous-reciprocal-extremity-exclusivity constraint in ReMixT karyotypes; and y is the total number of reciprocal locations that both have novel adjacencies in ReMixT karyotypes.

Figure 5.

Evidence of complex k-break (k ≥ 3) rearrangements in metastatic prostate cancer. (A) Two complex rearrangements across two genomes in a heterogeneous sample. A 5-break rearrangement that produced four novel adjacencies {a, b, c, d} involving five reference adjacencies (X, R, L, O, and M), with novel adjacency a not present in genome G₂. A 3-break rearrangement that produced three novel adjacencies {e, f, j} involving three reference adjacencies (Y, Z, and T), with novel adjacency j not present G₁. (B, top) A complex 5-break rearrangement on Chromosome 10 in the karyotype inferred by RCK on sample A31a. Only the four novel adjacencies, five reference adjacencies, and incident segments involved in the rearrangement are shown. Copy numbers ≤1 are omitted for clarity, and absent segments/adjacencies are shown as faded. (Bottom) The locations of the corresponding double-stranded DNA breakages for the 5-break on Chromosome 10, indicated as x|y for each reference adjacency {(x)^h, (y)^t}. Three reference adjacencies lie in/near genes: reference adjacency 102,756,[799|800] falls within the promoter region for gene LZTS2; reference adjacency 114,208,50[2|3] falls inside gene VTI1A; and reference adjacency 114,062,94[6|7] falls inside gene TECTB. (C) Number of complex k-break (k ≥ 3) rearrangements reported in RCK-reconstructed karyotypes using HATCHet and Battenberg copy number inputs with novel adjacency utilization parameter P = 0.9. Values of 0 are omitted for clarity.

Figure 6.

Segments, extremities, and copy number profiles for genomes. (A) A diploid reference genome r containing two pairs of homologous chromosomes: A Chromosomes are dark blue and dark green, and the homologous B Chromosomes are light blue and light green. Chromosomes are partitioned into consecutive segments labeled 1 through 12. (B, top) Reference genome r is a collection of concatenations of segments; the “flat” end of segment j corresponds to the tail extremity j^t, whereas the “pointy” end of each segment j corresponds to the head extremity j^h. Dashed lines correspond to reference adjacencies between adjacent extremities. The set $\mathcal{T}\left(\mathrm{R}\right)=\{1_{\mathrm{A}}^t,1_{\mathrm{B}}^t,5_{\mathrm{A}}^h,5_{\mathrm{B}}^h,6_{\mathrm{A}}^t,6_{\mathrm{B}}^t,{12}_{\mathrm{A}}^h,{12}_{\mathrm{B}}^h\}$ of extremities is the telomere set. (Bottom) The diploid segment copy number profile C_R = (a, b) for the genome R with colors (dark/light blue/green) corresponding to A/B labeled segments. (C, top) A derived genome g obtained via multiple large-scale rearrangements from the reference genome R. Red dashed lines correspond to novel adjacencies, for example, $\{3_{\mathrm{A}}^h,7_{\mathrm{B}}^h\}$. (Bottom) The diploid segment copy number profile C_G = (a, b) for the genome g with colors (dark/light blue/green) corresponding to A/B labeled segments. The set $\mathcal{T}(\mathrm{G})$ of telomeres in the derived genome G is identical to the set $\mathcal{T}\left(\mathrm{R}\right)$ of telomeres in the reference genome R.

Figure 7.

Ambiguity and errors in inferring segment copy number (SCN) profiles for a heterogeneous sample S = (G₁, G₂) under different assumptions about the sample composition. (A) A two-genome proper sample S = (G₁, G₂): each genome G_i ∈ S is depicted as collections of adjacent blocks (top), and the corresponding sequences of signed blocks (bottom). (B) The copy number profile c = [c₁, c₂, c₃, c₄] inferred under the assumption that the sample is homogeneous (i.e., comprised of a single derived genome) and the reference genome is haploid (i.e., each segment has only a single haplotype in the reference). Each value c_j is the weighted average of the sums of haplotype-specific (or allele-specific) copy numbers $a_{i,j}+b_{i,j}={\hat{c}}_{i,j}+{č}_{i,j}$ over the genomes G_i ∈ S. (C) Allele-specific copy number profiles $\hat{\mathbf{c}}=\left[{\hat{c}}_1,{\hat{c}}_2,{\hat{c}}_3,{\hat{c}}_4\right]$ and $\mathbf{č}=[{č}_1,{č}_2,{č}_3,{č}_4]$ inferred under the assumption that the sample is homogeneous and the reference genome is diploid (i.e., each segment has two haplotypes labeled A and B). Here, the entries ${\hat{c}}_j$ and ${č}_j$ for segment j are averages $\left({\hat{c}}_{1,j}+{\hat{c}}_{2,j}\right)/2$ and $\left({č}_{1,j}+{č}_{2,j}\right)/2$ of genome- and allele-specific copy number values. Note that the vectors $\hat{\mathbf{c}}$ and $\mathbf{č}$ do not preserve the true A/B label of each allele: dark blue are true counts of allele A and light blue are true counts of allele B. Here, segments 2 and 4 are flipped. (D) Genome-specific copy number profiles c₁ = [c_1,1, c_1,2, c_1,3, c_1,4] and c₂ = [c_2,1, c_2,2, c_2,3, c_2,4] inferred under the assumption that the sample is heterogeneous, but the reference genome is haploid. Here, the entry c_i,j for a segment j and genome G_i is the sum ${\hat{c}}_{i,j}+{č}_{i,j}$ of allele-specific copy number values in a genome G_i. (E) Allele- and genome-specific copy number matrices $\tilde{\mathbf{C}}=\left(\hat{\mathbf{C}}={\left[{\hat{\mathbf{c}}}_1,{\hat{\mathbf{c}}}_2,\dots ,{\hat{\mathbf{c}}}_n\right]}^T,\mathbf{Č}={\left[{\mathbf{č}}_1,{\mathbf{č}}_2,\dots ,{\mathbf{č}}_n\right]}^T\right)$ inferred under the assumption that the sample is heterogeneous and the reference genome is diploid. Segments 2 and 4 are flipped alleles: $\left({č}_{1,2},{\hat{c}}_{2,2}\right)=\left(a_{1,2},b_{2,2}\right)$ and $\left({č}_{1,4},{\hat{c}}_{2,4}\right)=\left(a_{1,4},b_{2,4}\right)$.

Figure 8.

A DIAG $D(\mathrm{R},{\tilde{\mathcal{A}}}_N)=\left(V,E\right)$ constructed on a set $\left\{1,2,\dots ,12\right\}$ of segments, and a set $\mathcal{A}\left(\mathrm{R}\right)\cup \mathcal{H}\left({\tilde{\mathcal{A}}}_N\right)$ of adjacencies. The set $\mathcal{A}\left(\mathrm{R}\right)$ corresponds to reference adjacencies in a diploid reference R shown in Figure 6B, and the set ${\tilde{\mathcal{A}}}_N=\left\{\left\{3^h,7^h\right\},\left\{2^h,9^h\right\},\left\{4^t,8^t\right\},\left\{4^h,4^h\right\},\left\{5^t,8^h\right\},\left\{3^t,{10}^t\right\},\left\{6^h,{11}^t\right\}\right\}$ represents unlabeled novel adjacencies that were measured from a derived genome G shown in Figure 6C. Squares indicate telomere vertices $\mathcal{T}\left(G\right)=\mathcal{T}\left(\mathrm{R}\right)\subseteq V$, and circles are non-telomere vertices. Solid edges correspond to segment edges in E_S, with dark blue/green edges corresponding to segments labeled A, and light blue/green edges corresponding to segments labeled B. Black-dashed edges are reference adjacency edges E_R, and red-dotted edges are novel adjacency edges E_N.

Figure 9.

Derivation of extremities and novel adjacencies for input to RCK and ReMixT. (A) An example of derivation of coordinates that resembles a reciprocal signature in measured unlabeled novel adjacencies on a chromosome a. Positions p₁ = (a, 100, + ) and p₂ = (a, 107, − ) have reciprocal signature (i.e., $|\mathrm{coor}{\mathrm{d}}_1-\mathrm{coor}{\mathrm{d}}_2|=7<50$ and $\mathrm{st}{\mathrm{r}}_1=-\ne \mathrm{st}{\mathrm{r}}_2=+$). Updated pair $\left\{{p^{\mathrm{\prime }}}_1=\left(a,103,+\right),{p^{\mathrm{\prime }}}_2=\left(a,104,-\right)\right\}$ of coordinates constitutes a reciprocal location. (B) An example of partitioning of a set $\mathcal{F}=\left\{f_1,f_2,f_3,f_4\right\}$ of fragments from allele-specific copy number calls into a set S = {s₁, s₂, s₃, s₄, s₅, s₆, s₇, s₈} of segments. Extremities of segments in S correspond to either preprocessed coordinates of unlabeled novel adjacencies (e.g., $s_1^h=p_1^{\mathrm{\prime }}$,$s_2^t=p_2^{\mathrm{\prime }}$) or to the extremities of fragments in $\mathcal{F}$ (e.g., $s_3^h=f_2^h$,$s_4^t=f_3^t$).