combine_miscalls_and_dropouts.RdThis is the version into which you pass in your own matrices of allele call probabilities W (corresponding to omega) and Ws (omega-star) and vectors of allele-specific dropout rates. Notation follows the paper for the most part. So see it for details.
combine_miscalls_and_dropouts(geno_names, D, W, Ws)the names of the genotypes in the order that CKMRsim has stored them.
vector of allele-specific dropout rates (delta) in the paper. Must be named with the names of the alleles, and should be in the correct order (that which CKMRsim uses).
the matrix omega of allelic call probabilities. Must have rownames and colnames that are the allele names, and have to be in the correct order.
the matrix omega* of allelic call probabilities. Gotta have rownames and colnames.
This function returns matrix C, which is the genotype-call probability matrix. The element in the s-th row and the t-th column gives the probability that a true genotype of index s is called as a genotype of index t. This matrix is returned with rownames and colnames which are the genotype names using the standard CKMRsim alelle separator of " / ".
The old version (general_allele_based_geno_err_model) required that the variables g, gp, h, and gp all be input by the user. This version just uses the names attribute of the geno_freqs variable to take care of that. In other words, you pass it a character vector of genotype names (constructed by CKMRsim using the standard allele separate, " / ", and this function creates the g, gp, h, and hp vectors that are needed.)
Note that this is vectorized over g, gp, h, and hp. Note that g, gp, h, and hp can be character vectors giving the allele names, so long as the vector D has names which are the allele names and W and Ws have rownames and colnames which are the allele names.
# simple example with 2-allele locus.
# first, look at the relevant information you would have about a locus,
# namely the allele frequencies and the genotype frequencies
example_L_biallelic
#> $freqs
#> A B
#> 0.6 0.4
#>
#> $geno_freqs
#> A / A A / B B / B
#> 0.36 0.48 0.16
#>
# Now, we create a matrix of allele-to-allele mis-call rates.
# Let's say that with probability 0.01 you miscall a true A as a B
# and with probability 0.03 you miscall a true B as an A. We set
# these as different just because it provides a better illustration here.
W <- matrix(c(0.99, 0.01, 0.03, 0.97), ncol = 2, byrow = TRUE)
rownames(W) <- c("A", "B")
colnames(W) <- c("A", "B")
# Have a look at that:
W
#> A B
#> A 0.99 0.01
#> B 0.03 0.97
# Ws is a matrix of similar form that specifies the miscall rates given that
# a dropout occurred. We will set it to be the same as W for this example:
Ws <- W
# Now we need a vector of allele-specific dropout rates. Let's say that
# A is very likely to drop out (0.05) and B, not so much (0.02)
D <- c(A = 0.05, B = 0.02)
# See what the resulting genotype-call probability matrix looks like:
combine_miscalls_and_dropouts(names(example_L_biallelic$geno_freqs), D, W, Ws)
#> A / A A / B B / B
#> A / A 0.981090 0.017820 0.001090
#> A / B 0.048921 0.893358 0.057721
#> B / B 0.002064 0.055872 0.942064
# now, let's see what it looks like if it is only the dropouts, no miscalls
W2 <- matrix(0, ncol = 2, nrow = 2)
diag(W2) <- 1
dimnames(W2) <- dimnames(W) # recall, you must have row and column names on these
Ws2 <- W2
combine_miscalls_and_dropouts(names(example_L_biallelic$geno_freqs), D, W2, Ws2)
#> A / A A / B B / B
#> A / A 1.00 0.00 0.00
#> A / B 0.02 0.93 0.05
#> B / B 0.00 0.00 1.00