Principles of phasing mutations to haplotypes: mutations are physically linked to SNPs on long reads, SNPs are phased to haplotypes, and the mutation inherits the haplotype of its linked SNP alleles.
phaseTime combines long-read and short-read whole genome sequencing (WGS) of tumour samples to
This tutorial walks through every step of the method on a simulated tumour. Alongside the plots produced here, schematic illustrations explain the principles, and example output from real multi-region long-read sequenced lung tumours shows what the results look like on real data.
Principles of phasing mutations to haplotypes: mutations are physically linked to SNPs on long reads, SNPs are phased to haplotypes, and the mutation inherits the haplotype of its linked SNP alleles.
# install.packages("remotes")
remotes::install_github("kerstinth/phaseTime", build_vignettes = TRUE)phaseTime only depends on ggplot2 and
patchwork.
phaseTime takes the following inputs for each tumour region (sample):
| input | columns | typical source |
|---|---|---|
mutations |
sample, chr, pos,
vaf (or var_count, depth),
optionally cluster |
Mutect2 calls from paired short-read WGS |
cn |
sample, chr, seg_start,
seg_end, major / minor copy number |
Battenberg on short-read WGS |
purity |
named vector, one value per sample | Battenberg |
snps |
region, chr, pos,
baf of heterozygous germline SNPs |
SNP calling / Battenberg allele counts |
linkage |
sample, chr, pos,
snp_pos, mt_snp_alt,
mt_snp_ref |
pileup of Nanopore long reads (Dorado, minimap2 to GRCh38) |
tree |
edges (parent, child) and trunk | CONIPHER on short-read WGS |
linkage contains, for every mutation and nearby
heterozygous SNP, the number of long reads that carry the mutation
and the alternative (mt_snp_alt) or
reference (mt_snp_ref) allele of the SNP.
To make this tutorial self-contained we simulate a tumour with two
regions with simulate_tumour(). Mutations are acquired on
individual DNA copies along the phylogeny 1 -> 2 -> 3
and 1 -> 4 -> 5 (cluster 2 is clonal in region 1,
cluster 4 in region 2), and copy number gains duplicate DNA copies
together with the mutations they carry. The simulated reads mimic a
typical long-read study (~200x short-read and ~17x long-read coverage).
Because the gain times are known, we can check the estimates against the
truth.
sim <- simulate_tumour(seed = 1)
names(sim)
#> [1] "mutations" "cn" "snps" "linkage" "purity" "tree"
#> [7] "truth"
head(sim$mutations[, c("sample", "chr", "pos", "ref", "alt", "var_count", "depth", "vaf", "cluster")])
#> sample chr pos ref alt var_count depth vaf cluster
#> 1 SIM001_R1 chr1 10914802 C G 131 211 0.6208531 1
#> 2 SIM001_R1 chr1 25932362 G A 131 213 0.6150235 1
#> 3 SIM001_R1 chr1 18487166 A T 140 208 0.6730769 1
#> 4 SIM001_R1 chr1 54304392 C G 126 203 0.6206897 1
#> 5 SIM001_R1 chr1 81994754 G C 104 192 0.5416667 1
#> 6 SIM001_R1 chr1 49001200 C G 118 189 0.6243386 1
head(sim$cn[, c("sample", "chr", "seg_start", "seg_end", "nMajor", "nMinor")])
#> sample chr seg_start seg_end nMajor nMinor
#> 1 SIM001_R1 chr1 1 102193308 3 1
#> 2 SIM001_R1 chr1 102193309 248956422 1 0
#> 3 SIM001_R1 chr2 1 242193529 1 1
#> 4 SIM001_R1 chr3 1 57739176 2 0
#> 5 SIM001_R1 chr3 57739177 154447013 4 1
#> 6 SIM001_R1 chr3 154447014 198295559 2 1
head(sim$snps[, c("region", "chr", "pos", "baf")])
#> region chr pos baf
#> 1 SIM001_R1 chr1 1784 0.7272727
#> 2 SIM001_R1 chr1 33879 0.7262570
#> 3 SIM001_R1 chr1 111767 0.7327189
#> 4 SIM001_R1 chr1 119144 0.7524752
#> 5 SIM001_R1 chr1 120435 0.7336683
#> 6 SIM001_R1 chr1 121217 0.7211538
head(sim$linkage)
#> sample chr pos snp_pos mt_snp_alt mt_snp_ref
#> 1 SIM001_R1 chr1 25932362 25933463 0 5
#> 2 SIM001_R1 chr1 25932362 25930705 6 0
#> 3 SIM001_R1 chr1 25932362 25930550 0 6
#> 4 SIM001_R1 chr1 18487166 18489363 0 9
#> 5 SIM001_R1 chr1 54304392 54302251 0 11
#> 6 SIM001_R1 chr1 81994754 81992765 1 8
sim$purity
#> SIM001_R1 SIM001_R2
#> 0.70 0.55
sim$tree
#> $edges
#> parent child
#> [1,] "1" "2"
#> [2,] "2" "3"
#> [3,] "1" "4"
#> [4,] "4" "5"
#>
#> $trunk
#> [1] "1"Long reads can carry both a somatic mutation and a heterozygous SNP. Simulating random positions of 50,000 mutations and 2.25 million heterozygous SNPs shows how much longer reads increase the fraction of mutations that can be phased: about 10% with 150 bp reads, but almost all with 5 kb reads.
phasing_sim <- simulate_snp_phasing(seed = 1)
phasing_sim
#> read_length number_mutations number_phased_muts pct_phased
#> 1 100 50000 3591 7.182
#> 2 150 50000 5391 10.782
#> 3 200 50000 6920 13.840
#> 4 300 50000 9926 19.852
#> 5 500 50000 15223 30.446
#> 6 1000 50000 26030 52.060
#> 7 2000 50000 38612 77.224
#> 8 5000 50000 48762 97.524
#> 9 7500 50000 49801 99.602
#> 10 10000 50000 49976 99.952
plot_snp_phasing_simulation(phasing_sim)
Example from real data: (A) a mutation and a heterozygous SNP on the same long reads; (B) the simulation above; (C-F) percentage of mutations phased to SNPs per tumour region and its relation to read length, coverage and number of heterozygous SNPs.
To reduce run time, only the 16 SNPs closest to each mutation (within
20 kb) are considered. select_phasing_snps() creates this
list of mutation-SNP pairs that is then counted in the long-read BAM
file (e.g. with samtools mpileup), discarding loci with
fewer than 3 reads, to produce the linkage table.
r1_snps <- unique(sim$snps[sim$snps$region == "SIM001_R1", c("chr", "pos")])
r1_muts <- sim$mutations[sim$mutations$sample == "SIM001_R1", ]
pairs <- select_phasing_snps(r1_muts, r1_snps, max_snps = 16, max_distance = 20000)
head(pairs)
#> chr pos snp_pos distance
#> 1 chr1 10914802 10912428 2374
#> 2 chr1 10914802 10930498 15696
#> 3 chr1 10914802 10931016 16214
#> 4 chr1 10914802 10932885 18083
#> 5 chr1 10914802 10933242 18440
#> 6 chr1 10914802 10933860 19058In a segment with allelic imbalance the B-allele frequency (BAF) of a heterozygous SNP deviates from 0.5. SNPs whose alternative allele is on the more abundant haplotype have BAF > 0.5; this haplotype is called haplotype A, the other haplotype B.
With several tumour regions, a segment can be phased as long as
one region is imbalanced: phase_snps() uses the
region with the strongest allelic imbalance as the reference and applies
its SNP haplotypes to all regions.
cn_long <- sim$cn
cn_long$region <- cn_long$sample
snps_phased <- phase_snps(sim$snps, cn = cn_long)
head(snps_phased[, c("region", "chr", "pos", "baf", "ref_region", "single_region_haplotype", "haplotype")])
#> region chr pos baf ref_region single_region_haplotype haplotype
#> 1 SIM001_R1 chr1 1784 0.7272727 SIM001_R1 A A
#> 2 SIM001_R1 chr1 33879 0.7262570 SIM001_R1 A A
#> 3 SIM001_R1 chr1 111767 0.7327189 SIM001_R1 A A
#> 4 SIM001_R1 chr1 119144 0.7524752 SIM001_R1 A A
#> 5 SIM001_R1 chr1 120435 0.7336683 SIM001_R1 A A
#> 6 SIM001_R1 chr1 121217 0.7211538 SIM001_R1 A AIn region 1, chromosomes 11, 15 and 22 are balanced (BAF ~ 0.5) and cannot be phased from region 1 alone (top), but region 2 carries a gain there, so the multi-region haplotypes phase them (bottom).
Example from real data: single-region (top) and multi-region (bottom) SNP phasing of two regions of the same tumour.
Battenberg reports major and minor copy numbers.
haplotype_copy_number() assigns them to haplotypes A and B
using the phased SNPs, so that copy number and mutations refer to the
same haplotypes in every region. For subclonal segments,
battenberg_average_cn() averages the clone states by their
cellular fractions (resulting in non-integer copy numbers).
cn_r1 <- haplotype_copy_number(sim$cn[sim$cn$sample == "SIM001_R1", ],
snps = snps_phased[snps_phased$region == "SIM001_R1", ])
head(cn_r1[, c("chr", "seg_start", "seg_end", "nMajor", "nMinor", "flip", "cnA", "cnB")])
#> chr seg_start seg_end nMajor nMinor flip cnA cnB
#> 1 chr1 1 102193308 3 1 FALSE 3 1
#> 2 chr1 102193309 248956422 1 0 FALSE 1 0
#> 3 chr2 1 242193529 1 1 FALSE 1 1
#> 4 chr3 1 57739176 2 0 FALSE 2 0
#> 5 chr3 57739177 154447013 4 1 FALSE 4 1
#> 6 chr3 154447014 198295559 2 1 FALSE 2 1Next, the mutation multiplicity (the number of DNA copies carrying a
mutation) is computed from the VAF, the purity p and the
total copy number CN_t (normal copy number
CN_n = 2):
\[n_{mut} = VAF \cdot \frac{1}{p} \left[p\,CN_t + CN_n (1 - p)\right]\]
muts_r1 <- annotate_mutations(sim$mutations[sim$mutations$sample == "SIM001_R1", ],
cn_r1, purity = sim$purity[["SIM001_R1"]])
head(muts_r1[, c("chr", "pos", "vaf", "cnA", "cnB", "tot_cn", "mut_cpn")])
#> chr pos vaf cnA cnB tot_cn mut_cpn
#> 1 chr1 10914802 0.6208531 3 1 4 3.015572
#> 2 chr1 25932362 0.6150235 3 1 4 2.987257
#> 3 chr1 18487166 0.6730769 3 1 4 3.269231
#> 4 chr1 54304392 0.6206897 3 1 4 3.014778
#> 5 chr1 81994754 0.5416667 3 1 4 2.630952
#> 6 chr1 49001200 0.6243386 3 1 4 3.032502For every SNP on the same long reads as a mutation, phaseTime computes the fraction of mutated reads that carry the SNP alternative allele. Using the SNP haplotype, this becomes the fraction of mutated reads on haplotype A (one minus the fraction for SNPs whose alternative allele is on haplotype B). The mean across SNPs decides the haplotype of the mutation: A if >= 0.5, otherwise B.
Phasing mutations using the haplotypes of physically linked SNPs.
phased_r1 <- phase_mutations(sim$linkage[sim$linkage$sample == "SIM001_R1", ], snps_phased, muts_r1)
table(phase = phased_r1$phase, useNA = "ifany")
#> phase
#> A B <NA>
#> 14716 5285 26160In regions of loss of heterozygosity one haplotype has copy number 0,
so any mutation phased to it must be wrong.
phasing_summary() uses this to estimate the phasing
error:
phasing_summary(phased_r1)
#> sample n_mutations n_phased pct_phased n_loh_phased n_incorrect
#> SIM001_R1 SIM001_R1 46161 20001 43.32878 3777 10
#> pct_incorrect pct_incorrect_loh
#> SIM001_R1 0.02166331 0.2647604The simulation also records the true haplotype, which confirms the accuracy:
plot_phased_mutations() shows five views of a tumour
region: copy number only, all mutations, mutations coloured by
haplotype, and haplotype A or B mutations only.
p_types <- c("copy_number", "all", "haplotype", "A", "B")
patchwork::wrap_plots(lapply(p_types, function(t) plot_phased_mutations(phased_r1, cn_r1, type = t)),
ncol = 1) + patchwork::plot_annotation(tag_levels = "A")
Example from real data: the same five views (A-E) for a lung tumour region.
All plots also work for single chromosomes:
On a gained segment, a mutation acquired before the
gain was copied with the segment, so its multiplicity equals the
haplotype copy number (early). A mutation acquired
after the gain is only on some copies
(late). classify_mutation_timing() rounds
multiplicities and compares them with the copy number of the haplotype
the mutation is on. Only clonal (integer) gained haplotypes can be
timed.
phased_r1 <- classify_mutation_timing(phased_r1)
table(phased_r1$timing)
#>
#> early late not_possible
#> 1368 11482 33311(plot_mutation_timing(phased_r1, chr = "chr5", colour_by = "haplotype") + ggtitle("A-C")) /
(plot_mutation_timing(phased_r1, chr = "chr5", haplotypes = "A") + ggtitle("D"))
Example from real data: mutations on one chromosome by haplotype (A-C) and early (orange) / late (purple) mutations on haplotype A (D).
Note how separating the haplotypes matters: without phasing, the haplotype B mutations at multiplicity 1 would be counted as late mutations of the haplotype A gain, making the gain appear earlier than it was.
On segments with several consecutive gains, multiplicities cluster at
integer values below the observed copy number. Each such cluster marks
an ancestral copy number state that existed before a later
gain. assign_cn_states() assigns mutations to the closest
integer state and removes states supported by fewer than 5
mutations.
On chromosome 5 the simulated haplotype A was gained twice (1 -> 2 -> 3), and the mutations form three states:
phased_r1 <- assign_cn_states(phased_r1, min_muts = 5)
table(state = phased_r1$cn_state[phased_r1$chr == "chr5" & phased_r1$phase %in% "A"])
#> state
#> 1 2 3
#> 1189 175 139
plot_cn_states(phased_r1, haplotype = "A", chr = "chr5")
Ancestral copy number states: (A, B) a gain history that results in mutations at multiplicities 1, 2 and 3 on a segment with copy number 3; (C) an example from real data.
The relative numbers of early and late mutations tell us when a gain occurred: few early mutations mean an early gain, many a late one.
Principles of timing copy number gains.
Because more DNA copies accumulate more mutations, each mutation count is divided by the number of copies present when the mutations were acquired. The timing of a gain is then the copy number adjusted number of mutations acquired before the gain, divided by all adjusted mutations. A value of 0 means the gain happened at the start of tumour evolution and 1 means it happened at the time of sampling. The 95% confidence interval is \(T \pm 1.96\sqrt{T(1-T)/N}\).
The worked example below (one mutation before the first gain, two between the gains, three after the second gain) gives timings of 1/3 and 2/3:
gain_timing_from_counts(states = c(3, 2, 1), counts = c(1, 2, 3))
#> gain_index copy_number_before_gain copy_number_after_gain
#> 1 1 1 2
#> 2 2 2 3
#> number_early_muts_cn_adjusted number_all_muts_cn_adjusted gain_timing
#> 1 1 3 0.3333333
#> 2 2 3 0.6666667
#> se loCI hiCI
#> 1 0.2721655 0.000000 0.866768
#> 2 0.2721655 0.133232 1.000000
Estimating the timing of first and subsequent gains: late mutations are divided by the number of copies present when they occurred.
simulate_gain_history() follows a DNA segment over
10,000 cell divisions. Mutations are drawn from a Poisson distribution
whose rate is proportional to the current copy number, and gains happen
after specified numbers of divisions, so the true timing of each gain is
known.
set.seed(11)
scenario_a <- simulate_gain_history(timepoints = c(7000, 8000, 10000), cn_steps = c(1, 3, 10, 10))
plot_simulated_gain_timing(scenario_a)
scenario_b <- simulate_gain_history(timepoints = c(1000, 2000, 5000, 10000), cn_steps = c(1, 2, 3, 5, 5))
simulated_gain_timing(scenario_b)[, c("gain_index", "copy_number_before_gain", "copy_number_after_gain",
"gain_timing", "loCI", "hiCI", "theoretical_gain_timing")]
#> gain_index copy_number_before_gain copy_number_after_gain gain_timing
#> 1 1 1 2 0.1172849
#> 2 2 2 3 0.2177100
#> 3 3 3 5 0.5123882
#> loCI hiCI theoretical_gain_timing
#> 1 0.09313829 0.1414314 0.1
#> 2 0.18673959 0.2486804 0.2
#> 3 0.47487703 0.5498994 0.5
plot_simulated_gain_timing(scenario_b)
Repeating this for 1,000 random scenarios (1 to 5 gains, copy numbers up to 10) at three mutation rates shows that the estimates are accurate (median absolute error < 0.01 at 1.33e-8 per bp) but become very uncertain at low mutation rates:
acc <- simulate_gain_timing_accuracy(n_sim = 1000, seed = 1)
tapply(acc$abs_difference, acc$mutation_rate_bp, median)
#> 1.33e-10 1.33e-09 1.33e-08
#> 0.058032626 0.019593310 0.005530813
tapply(acc$ci_width, acc$mutation_rate_bp, median)
#> 1.33e-10 1.33e-09 1.33e-08
#> 0.54051733 0.20864020 0.06542555
plot_gain_timing_accuracy(acc)
Reference results of the same simulations: example scenarios (A-C, E-G), absolute error (D) and confidence interval widths (H).
time_gains() brings everything together. For every
clonal gained segment and each haplotype it
To time clonal gains, only the mutations of the trunk cluster are used:
timing_r1 <- time_gains(phased_r1, cluster_col = "cluster", tree = sim$tree, sample = "SIM001_R1")
head(timing_r1[, c("segID", "haplotype", "cluster", "gain_index", "copy_number_before_gain",
"copy_number_after_gain", "n_mutations", "gain_timing", "loCI", "hiCI")])
#> segID haplotype cluster gain_index copy_number_before_gain
#> 1 chr1_1_102193308 A 1 1 1
#> 2 chr1_1_102193308 A 1 2 2
#> 3 chr3_1_57739176 A 1 1 1
#> 4 chr3_57739177_154447013 A 1 1 1
#> 5 chr3_57739177_154447013 A 1 2 2
#> 6 chr3_154447014_198295559 A 1 1 1
#> copy_number_after_gain n_mutations gain_timing loCI hiCI
#> 1 2 446 0.6427245 0.5854776 0.6999713
#> 2 3 446 0.7003096 0.6455804 0.7550388
#> 3 2 291 0.3743017 0.3034068 0.4451966
#> 4 2 721 0.3169811 0.2609592 0.3730031
#> 5 4 721 0.4811321 0.4209751 0.5412890
#> 6 2 194 0.4724409 0.3856138 0.5592681plot_gain_timing() produces a genome-wide summary: copy
number; haplotype A mutations (yellow: before the first gain, blue:
after the first gain, grey: not timed); haplotype A gain timing with 95%
confidence intervals; and the same for haplotype B.
Example from real data: gain timing output for a lung tumour region.
Since this is simulated data, the estimates can be compared with the true gain times:
key <- function(d) paste(d$segID, d$haplotype, d$cluster, d$copy_number_before_gain, d$copy_number_after_gain)
timing_r1$truth <- sim$truth$true_timing[match(key(timing_r1), key(sim$truth))]
ggplot(timing_r1[!is.na(timing_r1$truth), ], aes(truth, gain_timing, colour = factor(gain_index))) +
geom_abline(linetype = 2, colour = "grey50") +
geom_errorbar(aes(ymin = loCI, ymax = hiCI), width = 0) + geom_point() +
scale_colour_manual(name = "gain", values = gain_colours(4)[c(2, 4, 6)]) +
coord_equal(xlim = c(0, 1), ylim = c(0, 1)) +
labs(x = "simulated gain timing", y = "phaseTime estimate") + theme_bw()
First gains are recovered accurately. Second gains of the form 2 -> 3 are estimated somewhat too early in this simulation: here the gain duplicates only one of the two copies, so only half of the mutations acquired between the two gains reach the higher multiplicity.
run_phaseTime() runs all steps above for every region at
once, including the multi-region SNP phasing:
res <- run_phaseTime(sim$mutations, sim$cn, sim$purity, sim$snps, sim$linkage,
cluster_col = "cluster", tree = sim$tree)
names(res)
#> [1] "snps" "cn" "mutations" "timing" "phasing"
res$phasing
#> sample n_mutations n_phased pct_phased n_loh_phased n_incorrect
#> SIM001_R1 SIM001_R1 46161 20001 43.32878 3777 10
#> SIM001_R2 SIM001_R2 46494 19900 42.80122 3861 5
#> pct_incorrect pct_incorrect_loh
#> SIM001_R1 0.02166331 0.2647604
#> SIM001_R2 0.01075408 0.1295001The first clonal gains of all regions can be summarised as a heatmap, and gains timed independently in two regions of the same tumour can be compared:
comparison <- compare_region_timing(res$timing)
median(comparison$difference)
#> [1] 0.01909525
plot_region_comparison(comparison)
Example from real data: first clonal gain timing across ten regions of lung tumours (A) and concordance between regions of the same tumour (B, C).
With mutation clusters and a phylogenetic tree, gains can be timed
within the subclonal expansion in which they happened. Use
clonal_only = FALSE to also return gains of subclonal
clusters.
In the simulated tumour, haplotype A of chromosome 14 was gained in the trunk (1 -> 2 copies) after 44% of the trunk’s evolutionary time, and doubled again (2 -> 4 copies) in cluster 2 after 30% of cluster 2’s evolutionary time, while haplotype B was lost. The trunk mutations sit at multiplicity 4 (before the trunk gain) or 2 (after), and the cluster 2 mutations at 2 (before the cluster 2 doubling) or 1 (after):
timing_tree <- time_gains(phased_r1, cluster_col = "cluster", tree = sim$tree, clonal_only = FALSE)
timing_tree[timing_tree$segID == "chr14_1_107043718",
c("haplotype", "cluster", "copy_number_before_gain", "copy_number_after_gain", "gain_timing")]
#> haplotype cluster copy_number_before_gain copy_number_after_gain gain_timing
#> 18 A 1 1 2 0.4638158
#> 19 A 2 2 4 0.3879310
sim$truth[sim$truth$segID == "chr14_1_107043718", ]
#> segID haplotype cluster copy_number_before_gain
#> 24 chr14_1_107043718 A 1 1
#> 25 chr14_1_107043718 A 2 2
#> copy_number_after_gain true_timing
#> 24 2 0.442
#> 25 4 0.300
plot_timing_tree(sim$tree, timing_tree, "chr14_1_107043718", haplotype = "A")
Example from real data: two gains of the same segment, one in the trunk and one in a subclone, placed on the tumour phylogeny.
Gains seen only in region 2 (cluster 4) are placed on the edge leading to cluster 4:
m2 <- res$mutations[res$mutations$sample == "SIM001_R2", ]
t2 <- time_gains(m2, cluster_col = "cluster", tree = sim$tree, clonal_only = FALSE, sample = "SIM001_R2")
t2[t2$cluster == "4", c("segID", "haplotype", "copy_number_before_gain", "copy_number_after_gain", "gain_timing")]
#> segID haplotype copy_number_before_gain
#> 3 chr1_102193309_248956422 A 1
#> 13 chr6_1_107764224 A 1
#> 19 chr11_1_40323325 A 1
#> 23 chr15_49716226_101991189 A 1
#> 30 chr22_1_50818468 A 1
#> copy_number_after_gain gain_timing
#> 3 2 0.06315789
#> 13 2 0.42328042
#> 19 2 0.21621622
#> 23 2 0.12121212
#> 30 2 0.51612903
sim$truth[sim$truth$cluster == "4", ]
#> segID haplotype cluster copy_number_before_gain
#> 3 chr1_102193309_248956422 A 4 1
#> 13 chr6_1_107764224 A 4 1
#> 21 chr11_1_40323325 A 4 1
#> 28 chr15_49716226_101991189 A 4 1
#> 37 chr22_1_50818468 A 4 1
#> copy_number_after_gain true_timing
#> 3 2 0.071
#> 13 2 0.419
#> 21 2 0.265
#> 28 2 0.141
#> 37 2 0.549A typical analysis of Battenberg / CONIPHER output looks like this:
# copy number: average subclonal states, keep a sample column
cn <- do.call(rbind, lapply(regions, function(r) {
bb <- read.table(paste0(r, "_subclones.txt"), header = TRUE)
cbind(sample = r, battenberg_average_cn(bb))
}))
# alternatively use Battenberg's multi-region haplotypes directly:
# flips <- parse_battenberg_msai(read.table("multisample_MSAI.txt", header = TRUE))
# cn_r1 <- haplotype_copy_number(cn[cn$sample == "R1", ], flip = flips[flips$region == "R1", ])
purity <- c(R1 = 0.62, R2 = 0.48)
tree <- as_phaseTime_tree(readRDS("tumour.tree.RDS")) # CONIPHER tree object
res <- run_phaseTime(mutations, cn, purity, snps, linkage,
cluster_col = "cluster", tree = tree, output_dir = "phaseTime_output")adjust = "copies_present") mutation counts are divided by
the number of copies present while the mutations were acquired, as
described above and validated by the simulations. An alternative is to
divide the counts by the multiplicity of each state; use
adjust = "multiplicity" in time_gains() or
gain_timing_from_counts() for this.min_muts = 5 mutations are removed. With low sequencing
depth, noisy multiplicities can create spurious intermediate states;
increase min_muts in that case.detection_limit_cn()
gives the copy number above which a single copy mutation falls below a
given VAF:detection_limit_cn(purity = c(0.3, 0.5, 0.7, 0.9), min_vaf = 0.05)
#> purity max_copy_number
#> 1 0.3 15
#> 2 0.5 18
#> 3 0.7 19
#> 4 0.9 19sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.7.7
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] en_US/en_US/en_US/C/en_US/en_US
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplot2_4.0.0 phaseTime_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] patchwork_1.3.2 vctrs_0.6.5 cli_3.6.5 knitr_1.50
#> [5] rlang_1.1.6 xfun_0.53 generics_0.1.4 S7_0.2.0
#> [9] jsonlite_2.0.0 labeling_0.4.3 glue_1.8.0 htmltools_0.5.8.1
#> [13] sass_0.4.10 scales_1.4.0 rmarkdown_2.30 grid_4.5.1
#> [17] tibble_3.3.0 evaluate_1.0.5 jquerylib_0.1.4 fastmap_1.2.0
#> [21] yaml_2.3.10 lifecycle_1.0.4 compiler_4.5.1 dplyr_1.1.4
#> [25] RColorBrewer_1.1-3 pkgconfig_2.0.3 farver_2.1.2 digest_0.6.37
#> [29] R6_2.6.1 tidyselect_1.2.1 pillar_1.11.1 magrittr_2.0.4
#> [33] bslib_0.9.0 withr_3.0.2 tools_4.5.1 gtable_0.3.6
#> [37] cachem_1.1.0