intro_2.R
In segtest: Tests for Segregation Distortion in Polyploids

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
set.seed(1)

## ----setup--------------------------------------------------------------------
library(segtest)

## -----------------------------------------------------------------------------
drbounds(ploidy = 8) ## DR bounds
gamfreq(g = 4, ploidy = 8, type = "polysomic") ## no DR
gamfreq(g = 4, ploidy = 8, alpha = c(0.1, 0.01), type = "polysomic") ## Some DR

## -----------------------------------------------------------------------------
seg[seg$ploidy == 8 & seg$g == 4 & seg$mode %in% c("disomic", "both"), "p"]

## -----------------------------------------------------------------------------
n_pp_mix(g = 4, ploidy = 8)

## -----------------------------------------------------------------------------
gamfreq(g = 4, ploidy = 8, gamma = c(1, 0, 0), type = "mix")
gamfreq(g = 4, ploidy = 8, gamma = c(0, 1, 0), type = "mix")
gamfreq(g = 4, ploidy = 8, gamma = c(0, 0, 1), type = "mix")

## -----------------------------------------------------------------------------
gamfreq(g = 4, ploidy = 8, gamma = c(1, 1, 1)/3, type = "mix") 

## -----------------------------------------------------------------------------
n_pp_mix(g = 1, ploidy = 8)
gamfreq(g = 1, ploidy = 8, gamma = 1, type = "mix")
n_pp_mix(g = 7, ploidy = 8)
gamfreq(g = 7, ploidy = 8, gamma = 1, type = "mix")

## -----------------------------------------------------------------------------
beta_bounds(ploidy = 8)
gamfreq(g = 1, ploidy = 8, gamma = 1, beta = 0.03, type = "mix")
gamfreq(g = 7, ploidy = 8, gamma = 1, beta = 0.03, type = "mix")

## ----fig.alt="Bar plot of genotype frequencies."------------------------------
gf <- gf_freq(
  p1_g = 2, 
  p1_ploidy = 6,
  p1_gamma = c(0.7, 0.3), 
  p1_type = "mix",
  p2_g = 4,
  p2_ploidy = 6, 
  p2_gamma = c(0.5, 0.5), 
  p2_type = "mix")
plot(gf, type = "h", xlab = "Genotype", ylab = "Frequency")

## -----------------------------------------------------------------------------
x <- c(rmultinom(n = 1, size = 10, prob = gf))
x

## -----------------------------------------------------------------------------
gl <- simgl(nvec = x)
gl

## -----------------------------------------------------------------------------
## With known genotypes
sout1 <- seg_lrt(x = x, p1_ploidy = 6, p2_ploidy = 6, p1 = 2, p2 = 4)
sout1$p_value
## With genotype likelihoods
sout2 <- seg_lrt(x = gl, p1_ploidy = 6, p2_ploidy = 6, p1 = 2, p2 = 4)
sout2$p_value

## -----------------------------------------------------------------------------
gl10 <- gl / log(10)
seg_lrt(x = gl10, p1_ploidy = 6, p2_ploidy = 6, p1 = 2, p2 = 4)$p_value