inst/doc/Overview_of_the_junctions_package.R

In junctions: The Breakdown of Genomic Ancestry Blocks in Hybrid Lineages

## ----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 <- junctions::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 <- junctions::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 <- junctions::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 <- junctions::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------------------------------------------------------------
junctions::estimate_time(J = 10, N = 1000, R = 1000, H_0 = 0.5, C = 1)

## ----estimate time markers----------------------------------------------------
junctions::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))
junctions::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 <- junctions::estimate_time_diploid(ancestry_information = ancestry_matrix,
                                             analysis_type = "individuals",
                                             phased = TRUE,
                                             pop_size = 1000,
                                             freq_ancestor_1 = 0.5)

t_unphased <- junctions::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 <- junctions::log_likelihood_diploid(ancestry_matrix,
                                               pop_size = 1000,
                                               freq_ancestor_1 = 0.5,
                                               t = time_points,
                                               phased = TRUE)
ll_unphased <- junctions::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 <- junctions::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---------------------------------------------
junctions::number_of_junctions(N = 100, R = 100, H_0 = 0.5, C = 1, t = 1000)

## ----expected number of junctions markers-------------------------------------
junctions::number_of_junctions_markers(N = 100, H_0 = 0.5, t = 1000,
                                       marker_distribution = sort(runif(100, 0, 1)))

## ----expected number of junctions di------------------------------------------
junctions::number_of_junctions_di(N = 100, H_0 = 0.5, t = 1000, di = 1e-5)

## ----expected number of junctions back crossing-------------------------------
junctions::number_of_junctions_backcross(H_0 = 0.5, C = 1, t = 10)

## ----K------------------------------------------------------------------------
junctions::calc_k(N = 1000, R = 1000, H_0 = 0.5, C = 1)

## ----MAT----------------------------------------------------------------------
junctions::calculate_mat(N = 1000, R = 1000, H_0 = 0.5, C = 1)

## ----error--------------------------------------------------------------------
junctions::time_error(t = 30, N = 1000, R = 1000, H_0 = 0.5, C = 1)