Phased and unphased data
Thijs Janzen
2022-02-24
Phased and unphased data
In this vignette, we will explore how to analyze phased or unphased
data, using data simulated by the junctions’ package. We will go over
how the input data needs to be structured, and how to interpret the
results.
Simulating data
To have some example data, we will first simulate some artificial
data. We can do so as follows:
simulated_data <- sim_phased_unphased(pop_size = 1000,
freq_ancestor_1 = 0.5,
total_runtime = 100,
size_in_morgan = 1,
time_points = 100,
markers = 1000)
simulated_data
## # A tibble: 10,000 × 5
## time individual location anc_chrom_1 anc_chrom_2
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 100 0 0.00134 0 0
## 2 100 0 0.00217 0 0
## 3 100 0 0.00590 0 0
## 4 100 0 0.00694 1 0
## 5 100 0 0.00741 1 0
## 6 100 0 0.00799 1 0
## 7 100 0 0.00991 1 0
## 8 100 0 0.0101 1 0
## 9 100 0 0.0129 1 0
## 10 100 0 0.0131 1 0
## # … with 9,990 more rows
This returns a tibble with the following columns: time, individual,
location, anc_chrom_1 and anc_chrom_2. This contains ancestry on both
chromosomes for 10 sampled individuals, where 0 indicates an allele from
ancestor 0, and a 1 indicates an allele from ancestor 1. Locations are
given in Morgan.
Inferring the time since admixture
To infer the time since admixture, we need to provide the junctions
package with either our (in this case phased) ancestry data. Secondly,
we need to provide a vector with the location in Morgan. This vector is
converted into recombination distances between markers (very similar to
how Ancestry HMM treats markers as well). Simulation output already
contains this information and we can use this to infer the time since
admixture:
focal_data <- subset(simulated_data, simulated_data$time == 100)
admixture_time <- estimate_time_diploid(ancestry_information =
cbind(1,
1,
focal_data$location,
focal_data$anc_chrom_1,
focal_data$anc_chrom_2),
pop_size = 1000,
freq_ancestor_1 = 0.5,
phased = FALSE)
admixture_time
## # A tibble: 1 × 3
## individual time loglikelihood
## <dbl> <dbl> <dbl>
## 1 1 97.6 -2838.
Which returns two answers: the estimate (“minimum”) and the
-loglikelihood (“objective”)
Because the raw data is already phased (this comes easily with
simulated data), we can also use the phased framework. Here, we need to
provide a matrix with two columns, where each column reflects ancestry
in a specific chromosome. The size of the matrix (e.g. the number of
rows) should correspond to the length of the locations vector, where
these are locations in Morgan. Using again the simulated data, we
obtain:
morgan_locations <- focal_data$location
phased_data <- cbind(focal_data$anc_chrom_1, focal_data$anc_chrom_2)
admixture_time_phased <- estimate_time_diploid(ancestry_information =
cbind(1,
1,
morgan_locations,
phased_data),
phased = TRUE,
pop_size = 1000,
freq_ancestor_1 = 0.5)
admixture_time_phased
## # A tibble: 1 × 3
## individual time loglikelihood
## <dbl> <dbl> <dbl>
## 1 1 97.6 -3212.
Which yields a very similar time estimate.
Population size
Because we are using simulated data, we know beforehand the size of
the population. In reality, this is hardly ever the case. However, we
can explore the impact of population size on the time estimate, by
varying this systematically:
found <- c()
for (N in c(100, 1000, 10000, 100000, 1e6, 1e7)) {
admixture_time_phased <- estimate_time_diploid(ancestry_information =
cbind(1,
1,
morgan_locations,
phased_data),
phased = TRUE,
pop_size = N,
freq_ancestor_1 = 0.5)
found <- rbind(found, c(N, admixture_time_phased$time[[1]],
admixture_time_phased$loglikelihood[[1]]))
}
found
## [,1] [,2] [,3]
## [1,] 1e+02 93.4375 -3322.073
## [2,] 1e+03 97.5625 -3212.240
## [3,] 1e+04 95.6875 -3213.363
## [4,] 1e+05 95.6875 -3214.066
## [5,] 1e+06 95.6875 -3214.150
## [6,] 1e+07 95.6875 -3214.158
plot(found[, 2] ~ found[, 1], log = "x",
xlab = "Population Size",
ylab = "Admixture Time")
plot((-1 * found[, 3]) ~ found[, 1], log = "x",
xlab = "Population Size",
ylab = "Log Likelihood")
Because the loglikelihood surface is typically very flat for higher
population sizes, the population size estimate is typically not correct.
If we want to explore the effect of population size, we can also
directly calculate the likelihood (without optimization) and explore the
impact.
found <- c()
for (N in 10 ^ (seq(1, 6, length.out = 100))) {
ll <- junctions::log_likelihood_diploid(local_anc_matrix =
cbind(1,
morgan_locations,
phased_data),
phased = TRUE,
pop_size = N,
freq_ancestor_1 = 0.5,
t = 100)
found <- rbind(found, c(N, ll))
}
plot(found, xlab = "Population Size", ylab = "Log Likelihood", log = "x")
Which again shows that the likelihood increases strongly with
population size, but that past about 1000 individuals, there is no
significant difference anymore between population sizes. This is a
typical result, as the impact of 1/(2N) diminishes strongly with
increasing N and becomes much smaller than the impact of recombination
over time (see also Janzen et al. 2018).
Simulating data with error
To reflect phasing error, we can also simulate data with imposed
phasing error, e.g. with a fixed probability, the assignment of
chromosomes is swapped - for instance, if the true ancestry on
chromosome 1 is 0 for marker j, and is 1 on chromosome 2, with
probability m, the recording is swapped and appears as ancestry of 1 on
chromosome 1 and ancestry of 0 on chromosome 2. To simulate such data,
we can use the junctions package:
simulated_data <- sim_phased_unphased(pop_size = 1000,
freq_ancestor_1 = 0.5,
total_runtime = 100,
size_in_morgan = 1,
time_points = 100,
markers = 1000,
error_rate = 0.01)
simulated_data$true_data
## # A tibble: 10,000 × 5
## time individual location anc_chrom_1 anc_chrom_2
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 100 0 0.00105 0 0
## 2 100 0 0.00130 0 0
## 3 100 0 0.00378 1 0
## 4 100 0 0.00410 1 0
## 5 100 0 0.00736 1 0
## 6 100 0 0.00820 1 0
## 7 100 0 0.00847 1 0
## 8 100 0 0.00900 1 0
## 9 100 0 0.00961 1 0
## 10 100 0 0.0105 1 0
## # … with 9,990 more rows
simulated_data$phased_data
## # A tibble: 10,000 × 5
## time individual location anc_chrom_1 anc_chrom_2
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 100 0 0.00105 0 0
## 2 100 0 0.00130 0 0
## 3 100 0 0.00378 1 0
## 4 100 0 0.00410 1 0
## 5 100 0 0.00736 1 0
## 6 100 0 0.00820 1 0
## 7 100 0 0.00847 1 0
## 8 100 0 0.00900 1 0
## 9 100 0 0.00961 1 0
## 10 100 0 0.0105 1 0
## # … with 9,990 more rows
Now, the function returns the true data, and the phased data. We can
use either to infer the time since admixture (100 generations) and see
what error is induced:
focal_true_data <- subset(simulated_data$true_data,
simulated_data$true_data$individual == 0)
true_data <- cbind(focal_true_data$anc_chrom_1,
focal_true_data$anc_chrom_2)
true_loc <- focal_true_data$location
admixture_time_true <- estimate_time_diploid(ancestry_information =
cbind(1,
1,
true_loc,
true_data),
phased = TRUE,
pop_size = 1000,
freq_ancestor_1 = 0.5)
focal_phased_data <- subset(simulated_data$phased_data,
simulated_data$phased_data$individual == 0)
phased_data <- cbind(focal_phased_data$anc_chrom_1,
focal_phased_data$anc_chrom_2)
phased_loc <- focal_phased_data$location
admixture_time_error <- estimate_time_diploid(ancestry_information =
cbind(1,
1,
phased_loc,
phased_data),
phased = TRUE,
pop_size = 1000,
freq_ancestor_1 = 0.5)
admixture_time_true
## # A tibble: 1 × 3
## individual time loglikelihood
## <dbl> <dbl> <dbl>
## 1 1 107. -346.
## # A tibble: 1 × 3
## individual time loglikelihood
## <dbl> <dbl> <dbl>
## 1 1 150. -496.
Thus, we find that when using data with phasing error, a considerably
higher age is estimated, which is due to the introduction of “fake”
junctions as a result of incorrect phasing.