# tests/testthat/test-substitution.R In nodeSub: Simulate DNA Alignments Using Node Substitutions

```context("mutations")

test_that("mutations use", {

levels <- c("a", "c", "g", "t")
lbf <- length(levels)

# default is c(0.25, 0.25, 0.25, 0.25)
bf <- rep(1 / lbf, lbf)
Q <- rep(1, lbf * (lbf - 1) / 2) # default is JC69      # nolint

# only extract the 6 important rates.
if (is.matrix(Q)) Q <- Q[lower.tri(Q)] # nolint

eig_q <- phangorn::edQt(Q, bf) # eigen values
m <- length(levels) # always 4 (bases)

node_transition_matrix <- make_transition_matrix(mut_double = 0)
eigen_obj <- eigen(node_transition_matrix, FALSE)
eigen_obj\$inv <- solve.default(eigen_obj\$vec)

parentseq <- rep("a", 100)

found <- c()
for (repl in 1:100) {
rate <- 0.3
node_time <- 1

P <- nodeSub::slow_matrix(eig = eig_q,  # nolint
node_time,
rate = rate)

offspr1 <- parentseq
offspr2 <- parentseq

for (j in 1:m) {
ind <- parentseq == levels[j]
offspr1[ind] <- sample(levels, sum(ind), replace = TRUE, prob = P[, j])
offspr2[ind] <- sample(levels, sum(ind), replace = TRUE, prob = P[, j])
}
offspr1
offspr2
# total substitutions
sub1 <- sum(offspr1 != parentseq)
sub2 <- sum(offspr2 != parentseq)

# and now with the nodesub method:

p_matrix <- slow_matrix(eig = eigen_obj,
node_time,
rate = rate)
p_matrix_bl <- slow_matrix(eig = eig_q,
node_time,
rate = rate)

offsprings <- get_mutated_sequences(parentseq, p_matrix)

sub3 <- sum(offsprings[[1]] != parentseq)
sub4 <- sum(offsprings[[2]] != parentseq)
found <- rbind(found, cbind("regular", c(sub1, sub2)))
found <- rbind(found, cbind("nodesub", c(sub3, sub4)))
}

colnames(found) <- c("method", "substitutions")
found <- tibble::as_tibble(found)
found\$substitutions <- as.numeric(found\$substitutions)

vx <- found %>%
dplyr::group_by(method) %>%
dplyr::summarise("mean_sub" = mean(substitutions))

testthat::expect_equal(vx\$mean_sub[1], vx\$mean_sub[2] / 2,
tolerance = 0.5)
})

test_that("matrices", {

levels <- c("a", "c", "g", "t")
lbf <- length(levels)

# default is c(0.25, 0.25, 0.25, 0.25)
bf <- rep(1 / lbf, lbf)
Q <- rep(1, lbf * (lbf - 1) / 2) # nolint

# only extract the 6 important rates.
if (is.matrix(Q)) Q <- Q[lower.tri(Q)]  # nolint

eig_q <- phangorn::edQt(Q, bf) # eigen values

branch_length <- 1
rate <- 0.1

P1 <- nodeSub::slow_matrix(eig = eig_q,  # nolint
branch_length,
rate = rate)

P2 <- get_p_matrix(branch_length, eig_q, rate) # nolint

testthat::expect_equal(P1, P2)
})
```

## Try the nodeSub package in your browser

Any scripts or data that you put into this service are public.

nodeSub documentation built on Nov. 14, 2023, 5:10 p.m.