Nothing
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(fig.width = 7, fig.height = 7, echo = TRUE)
library(junctions)
require(tidyr)
require(dplyr)
require(ggplot2)
require(magrittr)
## ----phased unphased simulation-----------------------------------------------
simulated_pop <- sim_phased_unphased(pop_size = 1000,
freq_ancestor_1 = 0.5,
total_runtime = 100,
size_in_morgan = 1,
markers = sort(runif(n = 1000, 0, 1)),
time_points = seq(0, 100, by = 10))
## ----sim fin chrom------------------------------------------------------------
junctions_fin_chrom <- sim_fin_chrom(pop_size = 1000,
freq_ancestor_1 = 0.5,
total_runtime = 200,
morgan = 1,
R = 100) # the number of crossover sites
junctions_inf_chrom <- sim_inf_chrom(pop_size = 1000,
freq_ancestor_1 = 0.5,
total_runtime = 200,
morgan = 1,
markers = 100)
plot(junctions_inf_chrom$avgJunctions,
type = "l", lwd = 2,
col = "darkgreen",
xlab = "Time since admixture",
ylab = "Number of junctions")
lines(junctions_fin_chrom$avgJunctions, col = "blue", lwd = 2)
lines(junctions_inf_chrom$detectedJunctions, col = "lightgreen", lwd = 2)
legend("topleft", legend = c("Finite chromosome", "Infinite chromosome",
"Infinite chromosome detected"),
col = c("blue", "darkgreen", "lightgreen"), lty = 1, lwd = 2)
## ----backcrossing-------------------------------------------------------------
backcross_result <- sim_backcrossing(population_size = 1000,
freq_ancestor_1 = 0.5,
total_runtime = 20,
size_in_morgan = 1,
number_of_markers = 100)
plot(backcross_result$average_junctions,
type = "l", lwd = 2,
col = "darkgreen",
xlab = "Time since admixture",
ylab = "Number of junctions")
## ----estimate_time------------------------------------------------------------
estimate_time(J = 10, N = 1000, R = 1000, H_0 = 0.5, C = 1)
## ----estimate time markers----------------------------------------------------
estimate_time_one_chrom(J = 10, N = 1000, H_0 = 0.5,
marker_distribution = sort(runif(n = 1000, 0, 1)))
## ----estimate time haploid----------------------------------------------------
ancestry_data <- subset(simulated_pop, simulated_pop$time == 100)
ancestry_matrix <- dplyr::select(ancestry_data, c(individual,
location, anc_chrom_1))
estimate_time_haploid(ancestry_matrix = ancestry_matrix,
N = 1000,
freq_ancestor_1 = 0.5)
## ----estimate time diploid----------------------------------------------------
ancestry_matrix <- dplyr::select(ancestry_data, c(individual,
location,
anc_chrom_1, anc_chrom_2))
ancestry_matrix <- cbind(rep(1, length(ancestry_matrix$individual)),
ancestry_matrix)
t_phased <- estimate_time_diploid(ancestry_information = ancestry_matrix,
analysis_type = "individuals",
phased = TRUE,
pop_size = 1000,
freq_ancestor_1 = 0.5)
t_unphased <- estimate_time_diploid(ancestry_information = ancestry_matrix,
analysis_type = "individuals",
phased = FALSE,
pop_size = 1000,
freq_ancestor_1 = 0.5)
t_phased
t_unphased
## ----likelihood---------------------------------------------------------------
ancestry_matrix <- dplyr::select(ancestry_data, c(individual, location,
anc_chrom_1, anc_chrom_2))
time_points <- 80:120
ll_phased <- log_likelihood_diploid(ancestry_matrix,
pop_size = 1000,
freq_ancestor_1 = 0.5,
t = time_points,
phased = TRUE)
ll_unphased <- log_likelihood_diploid(ancestry_matrix,
pop_size = 1000,
freq_ancestor_1 = 0.5,
t = time_points,
phased = FALSE)
to_plot <- tibble::tibble(time_points, ll_phased, ll_unphased)
to_plot %>%
tidyr::gather(key = "phasing", value = "loglikelihood", -time_points) %>%
ggplot2::ggplot(
ggplot2::aes(x = time_points, y = loglikelihood, col = phasing)) +
ggplot2::geom_line()
## ----likelihood haploid-------------------------------------------------------
ancestry_matrix <- dplyr::select(ancestry_data, c(individual, location,
anc_chrom_1))
ll_haploid <- log_likelihood_haploid(ancestry_matrix,
N = 1000,
freq_ancestor_1 = 0.5,
t = time_points)
plot(ll_haploid ~ time_points, type = "l",
xlab = "Time since admixture",
ylab = "Loglikelihood")
## ----expected number of junctions---------------------------------------------
number_of_junctions(N = 100, R = 100, H_0 = 0.5, C = 1, t = 1000)
## ----expected number of junctions markers-------------------------------------
number_of_junctions_markers(N = 100, H_0 = 0.5, t = 1000,
marker_distribution = sort(runif(100, 0, 1)))
## ----expected number of junctions di------------------------------------------
number_of_junctions_di(N = 100, H_0 = 0.5, t = 1000, di = 1e-5)
## ----expected number of junctions back crossing-------------------------------
number_of_junctions_backcross(H_0 = 0.5, C = 1, t = 10)
## ----K------------------------------------------------------------------------
calc_k(N = 1000, R = 1000, H_0 = 0.5, C = 1)
## ----MAT----------------------------------------------------------------------
calculate_mat(N = 1000, R = 1000, H_0 = 0.5, C = 1)
## ----error--------------------------------------------------------------------
time_error(t = 30, N = 1000, R = 1000, H_0 = 0.5, C = 1)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.