In thoree/mut2: Reversible mutation matrices and pairwise likelihood ratios
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(pedmut, quietly = T)
library(pedsuite, quietly = T)
library(mut2, quietly = T)
library(xtable, quietly = T)
library(forrel, quietly = T)
Example 1: Tutorial examples
Consider a marker with $n = 4$ alleles with allele frequencies
$p = (0.1, 0.3, 0.4, 0.2)$. Assume an equal mutation model with mutation rate $0.003$, which in this coincides with the expected mutation rate, i.e.,
\begin{align}
Q=
\begin{pmatrix}
0.99700 & 0.00100 & 0.00100 & 0.00100 \
0.00100 & 0.99700 & 0.00100 & 0.00100 \
0.00100 & 0.00100 & 0.99700 & 0.00100 \
0.00100 & 0.00100 & 0.00100 & 0.99700 \
\end{pmatrix}
.
\end{align}
The proposed mutation matrix and the allele frequency vector,
$(Q,p)$, does not correspond to a balanced mutation model. For instance
[
p_2 q_{21} = 0.3\cdot 0.001 = 0.0003 \neq p_1q_{12} = 0.1\cdot 0.001 = 0.0001.
]
In this case the proposed matrix $Q$ is symmetric and
therefore default balancing simplifies to
\begin{align}
p_{ij}=q_{ij}\min {1, \frac{p_j}{p_i} }, \;i \neq j;
p_{ii}=1-\sum_{i\neq j} p_{ij},
\end{align}
Balancing is achieved by
Q = mutationMatrix("equal", rate = 0.003, alleles = 1:4)
p = c(0.1,0.3,0.4,0.2)
names(p) = 1:4
P1 = makeReversible(Q, afreq = p)
P1
The pair $(P_1,p)$ is balanced with expected mutation rate 0.002, which is lower than the expected for Q, 0.003. In this case the condition
$q_{ji} < p_i,\, j\neq i$ holds and we can transform using the alternative transformation (method = PM)
\begin{align}
p_{ij}=\frac{p_iq_{ij}+p_j q_{ji}}{2p_i},\;i \neq j,\;
p_{ii}=1-\sum_{i\neq j} p_{ij},
\end{align}
which preserves the expected mutation rate:
P2 = makeReversible(Q, afreq = p, method = "PM")
P2
Alternatively, the first balanced mutation matrix can be adjusted to have the required
expected mutation rate
P1.adj = adjustReversible(Q, P1, afreq = p)
P1.adj
There is also third option for balancing, the Barker option
P3 = makeReversible(Q, afreq = p, method = "BA")
P3
Checking mutation matrices
Create matrices based on database `NorwegianFrequencies'.
db = NorwegianFrequencies
nAlleles = unlist(lapply(db, length))
mN = names(db)
M = lapply(db, function(x){
mat = mutationMatrix("proportional",
alleles = names(x),
rate = 0.003,
rate2 = 0.000001,
range = 0.1,
afreq = x)
mat
})
Check matrices
gamma = lapply(M, function(x) expectedMutationRate(x, afreq = NULL))
numerator = lapply(M, function(x) {p = attr(x, "afreq");1-sum(p^2)})
K = unlist(gamma)/unlist(numerator)
res = lapply(M, function(x) reasonableMutationMatrix(x, p = NULL, verbose = F))
names(res) = mN
R = lapply(M, function(x) makeReversible(x, method = "BA"))
res2 = lapply(R, function(x) reasonableMutationMatrix(x, p = NULL, verbose = F))
names(res2) = mN
sum(unlist(res2))
Example 3: Comparing transformations
We first compare PM to MH.
db = NorwegianFrequencies
nAlleles = unlist(lapply(db, length))
mN = names(db)
mN = mN[nAlleles < 11]
db= db[mN]
M = lapply(db, function(x){
mat = mutationMatrix("stepwise",
alleles = names(x),
rate = 0.001,
rate2 = 1e-05,
range = 0.1,
afreq = x)
mat
})
R = lapply(M, function(x) makeReversible(x, method = "PM"))
resPM = compareLRR(M, R)
R = lapply(M, function(x) makeReversible(x, method = "MH"))
resMH = compareLRR(M, R)
res1 = cbind(t(resPM), t(resMH))
colnames(res1) = c("PMmin", "PMmuM", "PMsigma", "PMmax",
"MHmin", "MHmuM", "MHsigma", "MHmax")
tab1 = xtable(res1, digits = c(2,3,4,2,2,3,4,2,1), caption = "The")
We then compare PM to MH on a log scale.
db = NorwegianFrequencies
nAlleles = unlist(lapply(db, length))
mN = names(db)
mN = mN[nAlleles < 11]
db= db[mN]
M = lapply(db, function(x){
mat = mutationMatrix("stepwise",
alleles = names(x),
rate = 0.001,
rate2 = 1e-05,
range = 0.1,
afreq = x)
mat
})
R = lapply(M, function(x) makeReversible(x, method = "PM"))
resPM = compareLRR(M, R, ln = TRUE)
R = lapply(M, function(x) makeReversible(x, method = "MH"))
resMH = compareLRR(M, R, ln = TRUE)
res2 = cbind(t(resPM), t(resMH))
colnames(res2) = c("PMmin", "PMmuM", "PMsigma", "PMmax",
"MHmin", "MHmuM", "MHsigma", "MHmax")
tab2 = xtable(res2, digits = c(2,3,4,2,2,3,4,2,1), caption = "The")
We finally compare BA to MH, both adjusted, on a log scale.
db = NorwegianFrequencies
nAlleles = unlist(lapply(db, length))
mN = names(db)
mN = mN[nAlleles < 11]
n = length(mN)
db = db[mN]
M = lapply(db, function(x){
mat = mutationMatrix("stepwise",
alleles = names(x),
rate = 0.001,
rate2 = 1e-05,
range = 0.1,
afreq = x)
mat
})
R = lapply(M, function(x) makeReversible(x, method = "BA"))
RA = R
for(i in 1:n)
RA[[i]] = adjustReversible(M[[i]], R[[i]])
resBA = compareLRR(M, RA, ln = TRUE)
R = lapply(M, function(x) makeReversible(x, method = "MH"))
for(i in 1:n)
RA[[i]] = adjustReversible(M[[i]], R[[i]])
resMH = compareLRR(M, RA, ln = TRUE)
res3 = cbind(t(resBA), t(resMH))
colnames(res3) = c("BAmin", "BAmuM", "BAsigma", "BAmax",
"MHmin", "MHmuM", "MHsigma", "MHmax")
tab3 = xtable(res3, digits = c(2,3,4,2,2,3,4,2,1), caption = "The")
We calculate the root mean square errors:
RMSQ.PM = sqrt(sum(res1[,"PMmuM"]^2) + sum(res1[,"PMsigma"]^2))
RMSQ.BA = sqrt(sum(res3[,"BAmuM"]^2) + sum(res3[,"BAsigma"]^2))
RMSQ.MH = sqrt(sum(res3[,"MHmuM"]^2) + sum(res3[,"MHsigma"]^2))
data.frame(RMSQ.MH, RMSQ.PM, RMSQ.BA)
Example 3: Extreme cases
Trivial mutation matrix
Assume zero mutation rate. Then balancing leaves the matrix unchanged as it should:
Q = mutationMatrix("trivial", alleles = 1:4)
p = c(0.1,0.3,0.4,0.2)
names(p) = 1:4
P= makeReversible(Q, afreq = p)
all(diag(P) == 1)
P =makeReversible(Q, afreq = p, method = "PM")
all(diag(P) == 1)
Always mutation
Consider
Q = mutationMatrix("custom", matrix = matrix(c(0,1,1,0), ncol = 2),
alleles = 1:2)
p = c(0.1,0.9)
P = makeReversible(Q, afreq = p)
P
The two other methods, PM and BA do not work in this case.
Equally likely mutations
In this case only "MH" and "BA" work
p = c(0.1,0.9)
Q = mutationMatrix("custom",
matrix = matrix(rep(0.5,4), ncol = 2),
alleles = 1:2, afreq = p)
attr(Q, "rate") = expectedMutationRate(Q, p)
P1 = makeReversible(Q, afreq = p)
P3 = makeReversible(Q, afreq = p, method = "BA")
cat("Expected mutation rates:", "\n")
c("R.Q" = attr(Q, "rate"), "R.MH" = attr(P1, "rate"),
"R.BA" = attr(P3, "rate") )
One step model
Q = mutationMatrix("onestep", rate = 0.003, alleles = 1:4)
p = c(0.1,0.3,0.4,0.2)
names(p) = 1:4
P1 = makeReversible(Q, afreq = p)
P2 = makeReversible(Q, afreq = p, method = "BA")
SNPS
The input of the below function specifies
the allele frequency $p = (p_1, p_2)$ and a mutation matrix
\begin{align}
Q =
\begin{pmatrix}
1-q_{12} & q_{12} \
q_{21} & 1-q_{21}
\end{pmatrix}
\end{align}
LR.Q, the LR comparing "AF father of CH" to "AF and CH are unrelated", is
calculated for all genotype combinations.
The LR calculation is next done with Q replaced by a reversible
mutation matrix with the same expected mutation rate giving LR.M.
The ratio LRR = LR.M/LR.Q is found and also its expected value, ELRR.
LRR = function(p1 = 0.5, q12 = 0.05, q21 = 0.05){
p = c("1" = p1, "2" = 1 - p1)
q = matrix(c(1 - q12, q21, q12, 1 - q21), ncol = 2)
Q = mutationMatrix("custom", matrix = q, afreq = p, alleles = 1:2)
r = expectedMutationRate(Q, p)
# Reversible mutation model A with same mutation rate:
if (max(r/(2*p[1]), r/(2*p[2])) > 1)
stop("Reversible, mutation preserved matrix does not exist")
A = matrix(c(1-r/(2*p[1]), r/(2*p[2]), r/(2*p[1]), 1-r/(2*p[2])), ncol = 2)
x1 = nuclearPed(father = "AF", mother = "MO", child = "CH")
m = marker(x1, afreq = p)
x1 = setMarkers(x1, m)
numerator = oneMarkerDistribution(x1, partialmarker = 1, ids = c("AF", "MO"),
verbose = F)
x1 = setMutationModel(x1, "custom", matrix = A)
L1 = oneMarkerDistribution(x1, partialmarker = 1, ids = c("AF", "CH"),
verbose = F)
x2 = setMutationModel(x1, "custom", matrix = Q)
L2 = oneMarkerDistribution(x2, partialmarker = 1, ids = c("AF", "CH"),
verbose = F)
ELRR = sum(L1^2/L2)
list(ELRR = ELRR,
tab = list(LRM = L1/numerator,
LRQ = L2/numerator, LRR = L1/L2))
}
The default corresponds to a stationary model and hence we get an expected LRR equal to 1:
LRR()$ELRR
Next follows an example where the input matrix is far from stationary:
LRR(p1 = 0.01, q12 = 0.001, q21 = 0.001)
thoree/mut2 documentation built on May 16, 2023, 7:56 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(pedmut, quietly = T) library(pedsuite, quietly = T) library(mut2, quietly = T) library(xtable, quietly = T) library(forrel, quietly = T)
Example 1: Tutorial examples
Consider a marker with $n = 4$ alleles with allele frequencies $p = (0.1, 0.3, 0.4, 0.2)$. Assume an equal mutation model with mutation rate $0.003$, which in this coincides with the expected mutation rate, i.e., \begin{align} Q= \begin{pmatrix} 0.99700 & 0.00100 & 0.00100 & 0.00100 \ 0.00100 & 0.99700 & 0.00100 & 0.00100 \ 0.00100 & 0.00100 & 0.99700 & 0.00100 \ 0.00100 & 0.00100 & 0.00100 & 0.99700 \ \end{pmatrix} . \end{align} The proposed mutation matrix and the allele frequency vector, $(Q,p)$, does not correspond to a balanced mutation model. For instance [ p_2 q_{21} = 0.3\cdot 0.001 = 0.0003 \neq p_1q_{12} = 0.1\cdot 0.001 = 0.0001. ] In this case the proposed matrix $Q$ is symmetric and therefore default balancing simplifies to \begin{align} p_{ij}=q_{ij}\min {1, \frac{p_j}{p_i} }, \;i \neq j; p_{ii}=1-\sum_{i\neq j} p_{ij}, \end{align} Balancing is achieved by
Q = mutationMatrix("equal", rate = 0.003, alleles = 1:4) p = c(0.1,0.3,0.4,0.2) names(p) = 1:4 P1 = makeReversible(Q, afreq = p) P1
The pair $(P_1,p)$ is balanced with expected mutation rate 0.002, which is lower than the expected for Q, 0.003. In this case the condition
$q_{ji} < p_i,\, j\neq i$ holds and we can transform using the alternative transformation (method = PM)
\begin{align}
p_{ij}=\frac{p_iq_{ij}+p_j q_{ji}}{2p_i},\;i \neq j,\;
p_{ii}=1-\sum_{i\neq j} p_{ij},
\end{align}
which preserves the expected mutation rate:
P2 = makeReversible(Q, afreq = p, method = "PM") P2
Alternatively, the first balanced mutation matrix can be adjusted to have the required expected mutation rate
P1.adj = adjustReversible(Q, P1, afreq = p) P1.adj
There is also third option for balancing, the Barker option
P3 = makeReversible(Q, afreq = p, method = "BA") P3
Checking mutation matrices
Create matrices based on database `NorwegianFrequencies'.
db = NorwegianFrequencies nAlleles = unlist(lapply(db, length)) mN = names(db) M = lapply(db, function(x){ mat = mutationMatrix("proportional", alleles = names(x), rate = 0.003, rate2 = 0.000001, range = 0.1, afreq = x) mat })
Check matrices
gamma = lapply(M, function(x) expectedMutationRate(x, afreq = NULL)) numerator = lapply(M, function(x) {p = attr(x, "afreq");1-sum(p^2)}) K = unlist(gamma)/unlist(numerator) res = lapply(M, function(x) reasonableMutationMatrix(x, p = NULL, verbose = F)) names(res) = mN R = lapply(M, function(x) makeReversible(x, method = "BA")) res2 = lapply(R, function(x) reasonableMutationMatrix(x, p = NULL, verbose = F)) names(res2) = mN sum(unlist(res2))
Example 3: Comparing transformations
We first compare PM to MH.
db = NorwegianFrequencies nAlleles = unlist(lapply(db, length)) mN = names(db) mN = mN[nAlleles < 11] db= db[mN] M = lapply(db, function(x){ mat = mutationMatrix("stepwise", alleles = names(x), rate = 0.001, rate2 = 1e-05, range = 0.1, afreq = x) mat }) R = lapply(M, function(x) makeReversible(x, method = "PM")) resPM = compareLRR(M, R) R = lapply(M, function(x) makeReversible(x, method = "MH")) resMH = compareLRR(M, R) res1 = cbind(t(resPM), t(resMH)) colnames(res1) = c("PMmin", "PMmuM", "PMsigma", "PMmax", "MHmin", "MHmuM", "MHsigma", "MHmax") tab1 = xtable(res1, digits = c(2,3,4,2,2,3,4,2,1), caption = "The")
We then compare PM to MH on a log scale.
db = NorwegianFrequencies nAlleles = unlist(lapply(db, length)) mN = names(db) mN = mN[nAlleles < 11] db= db[mN] M = lapply(db, function(x){ mat = mutationMatrix("stepwise", alleles = names(x), rate = 0.001, rate2 = 1e-05, range = 0.1, afreq = x) mat }) R = lapply(M, function(x) makeReversible(x, method = "PM")) resPM = compareLRR(M, R, ln = TRUE) R = lapply(M, function(x) makeReversible(x, method = "MH")) resMH = compareLRR(M, R, ln = TRUE) res2 = cbind(t(resPM), t(resMH)) colnames(res2) = c("PMmin", "PMmuM", "PMsigma", "PMmax", "MHmin", "MHmuM", "MHsigma", "MHmax") tab2 = xtable(res2, digits = c(2,3,4,2,2,3,4,2,1), caption = "The")
We finally compare BA to MH, both adjusted, on a log scale.
db = NorwegianFrequencies nAlleles = unlist(lapply(db, length)) mN = names(db) mN = mN[nAlleles < 11] n = length(mN) db = db[mN] M = lapply(db, function(x){ mat = mutationMatrix("stepwise", alleles = names(x), rate = 0.001, rate2 = 1e-05, range = 0.1, afreq = x) mat }) R = lapply(M, function(x) makeReversible(x, method = "BA")) RA = R for(i in 1:n) RA[[i]] = adjustReversible(M[[i]], R[[i]]) resBA = compareLRR(M, RA, ln = TRUE) R = lapply(M, function(x) makeReversible(x, method = "MH")) for(i in 1:n) RA[[i]] = adjustReversible(M[[i]], R[[i]]) resMH = compareLRR(M, RA, ln = TRUE) res3 = cbind(t(resBA), t(resMH)) colnames(res3) = c("BAmin", "BAmuM", "BAsigma", "BAmax", "MHmin", "MHmuM", "MHsigma", "MHmax") tab3 = xtable(res3, digits = c(2,3,4,2,2,3,4,2,1), caption = "The")
We calculate the root mean square errors:
RMSQ.PM = sqrt(sum(res1[,"PMmuM"]^2) + sum(res1[,"PMsigma"]^2)) RMSQ.BA = sqrt(sum(res3[,"BAmuM"]^2) + sum(res3[,"BAsigma"]^2)) RMSQ.MH = sqrt(sum(res3[,"MHmuM"]^2) + sum(res3[,"MHsigma"]^2)) data.frame(RMSQ.MH, RMSQ.PM, RMSQ.BA)
Example 3: Extreme cases
Trivial mutation matrix
Assume zero mutation rate. Then balancing leaves the matrix unchanged as it should:
Q = mutationMatrix("trivial", alleles = 1:4) p = c(0.1,0.3,0.4,0.2) names(p) = 1:4 P= makeReversible(Q, afreq = p) all(diag(P) == 1) P =makeReversible(Q, afreq = p, method = "PM") all(diag(P) == 1)
Always mutation
Consider
Q = mutationMatrix("custom", matrix = matrix(c(0,1,1,0), ncol = 2), alleles = 1:2) p = c(0.1,0.9) P = makeReversible(Q, afreq = p) P
The two other methods, PM and BA do not work in this case.
Equally likely mutations
In this case only "MH" and "BA" work
p = c(0.1,0.9) Q = mutationMatrix("custom", matrix = matrix(rep(0.5,4), ncol = 2), alleles = 1:2, afreq = p) attr(Q, "rate") = expectedMutationRate(Q, p) P1 = makeReversible(Q, afreq = p) P3 = makeReversible(Q, afreq = p, method = "BA") cat("Expected mutation rates:", "\n") c("R.Q" = attr(Q, "rate"), "R.MH" = attr(P1, "rate"), "R.BA" = attr(P3, "rate") )
One step model
Q = mutationMatrix("onestep", rate = 0.003, alleles = 1:4) p = c(0.1,0.3,0.4,0.2) names(p) = 1:4 P1 = makeReversible(Q, afreq = p) P2 = makeReversible(Q, afreq = p, method = "BA")
SNPS
The input of the below function specifies the allele frequency $p = (p_1, p_2)$ and a mutation matrix
\begin{align} Q = \begin{pmatrix} 1-q_{12} & q_{12} \ q_{21} & 1-q_{21} \end{pmatrix} \end{align}
LR.Q, the LR comparing "AF father of CH" to "AF and CH are unrelated", is calculated for all genotype combinations. The LR calculation is next done with Q replaced by a reversible mutation matrix with the same expected mutation rate giving LR.M. The ratio LRR = LR.M/LR.Q is found and also its expected value, ELRR.
LRR = function(p1 = 0.5, q12 = 0.05, q21 = 0.05){ p = c("1" = p1, "2" = 1 - p1) q = matrix(c(1 - q12, q21, q12, 1 - q21), ncol = 2) Q = mutationMatrix("custom", matrix = q, afreq = p, alleles = 1:2) r = expectedMutationRate(Q, p) # Reversible mutation model A with same mutation rate: if (max(r/(2*p[1]), r/(2*p[2])) > 1) stop("Reversible, mutation preserved matrix does not exist") A = matrix(c(1-r/(2*p[1]), r/(2*p[2]), r/(2*p[1]), 1-r/(2*p[2])), ncol = 2) x1 = nuclearPed(father = "AF", mother = "MO", child = "CH") m = marker(x1, afreq = p) x1 = setMarkers(x1, m) numerator = oneMarkerDistribution(x1, partialmarker = 1, ids = c("AF", "MO"), verbose = F) x1 = setMutationModel(x1, "custom", matrix = A) L1 = oneMarkerDistribution(x1, partialmarker = 1, ids = c("AF", "CH"), verbose = F) x2 = setMutationModel(x1, "custom", matrix = Q) L2 = oneMarkerDistribution(x2, partialmarker = 1, ids = c("AF", "CH"), verbose = F) ELRR = sum(L1^2/L2) list(ELRR = ELRR, tab = list(LRM = L1/numerator, LRQ = L2/numerator, LRR = L1/L2)) }
The default corresponds to a stationary model and hence we get an expected LRR equal to 1:
LRR()$ELRR
Next follows an example where the input matrix is far from stationary:
LRR(p1 = 0.01, q12 = 0.001, q21 = 0.001)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.