graph_to_pwts = function(graph, leaves) {
# input igraph object
# output: numedges*(numpops-1) matrix 'pwts' which indicates all paths from pop to outpop; admixture edges are mapped onto parent edges and weighted
root = get_root(graph)
pwts = matrix(0, length(E(graph)), length(leaves))
colnames(pwts) = leaves
rownames(pwts) = attr(E(graph), 'vnames')
admixnodes = which(degree(graph, mode='in') == 2)
admixedges = unlist(incident_edges(graph, admixnodes, mode='in'))
allpaths = all_simple_paths(graph, root, leaves, mode='out')
pathcounts = table(names(sapply(allpaths, function(x) tail(x, 1))))
for(i in seq_len(length(allpaths))) {
target = names(tail(allpaths[[i]],1))
ln = length(allpaths[[i]])
pth2 = allpaths[[i]][c(1, 1+rep(seq_len(ln-2), each=2), ln)]
rowind = as.vector(E(graph)[get.edge.ids(graph, pth2)])
pwts[rowind,target] = pwts[rowind,target] + 1/pathcounts[target]
}
if(!is.null(admixedges)) pwts = pwts[-admixedges,]
#pwts[,-1] - pwts[,1]
pwts
}
tree_to_pwts = function(graph, leaves) {
# should give same output as graph_to_pwts for trees, but faster
# returns numedges*(numpops-1) matrix
# in tree, each edge maps to one node
pwts = matrix(0, length(E(graph))+1, length(leaves))
colnames(pwts) = leaves
rownames(pwts) = names(V(graph))
enames = attr(E(graph), 'vnames') %>% str_replace('.+\\|', '')
paths = shortest_paths(graph, 'R', leaves)$vpath
indmat = paths %>% map(as.numeric) %>% map(tail, -1) %>% imap(~cbind(.x, .y)) %>% do.call(rbind, .)
pwts[indmat] = 1
pwts = pwts[,-1] - pwts[,1]
pwts[enames,]
}
expand_path = function(path) {
# igraph path (sequence of vertices)
# duplicates inner vertices, so that this works with igraph functions that process vertex sequences pairwise
ln = length(path)
path[c(1, 1+rep(seq_len(ln-2), each=2), ln)]
}
graph_to_weightind = function(graph) {
# input igraph object
# output:
# map (leaf, edge) -> paths
# map path -> weights
# ultimately: indices for weights into paths, indices for paths into pwts
# room for improvement here...
leaves = get_leaves(graph)
root = get_root(graph)
admixnodes = which(degree(graph, mode='in') == 2)
admixedges = unlist(incident_edges(graph, admixnodes, mode='in'))
normedges = setdiff(1:length(E(graph)), admixedges)
paths = all_simple_paths(graph, root, leaves, mode='out')
ends = sapply(paths, tail, 1)
edge_per_path = paths %>% map(expand_path) %>% map(~get.edge.ids(graph, .))
weight_per_path = edge_per_path %>% map(~(which(admixedges %in% .)))
path_edge_table = do.call(rbind, lapply(seq_len(length(weight_per_path)),
function(i) tibble(path=i, edge=c(edge_per_path[[i]])))) %>%
mutate(edge2 = match(edge, normedges)) %>% filter(!is.na(edge2)) %>%
mutate(leaf = as.vector(ends[path]), leaf2 = match(leaf, leaves)) %>%
left_join(enframe(c(table(ends))) %>% transmute(leaf=as.numeric(name), numpaths=value), by='leaf') %>%
group_by(leaf2, edge2) %>% mutate(cnt = n(), keep = cnt < numpaths) %>% filter(keep) %>% as.matrix
path_admixedge_table = do.call(rbind, lapply(seq_len(length(weight_per_path)),
function(i) tibble(path=i, admixedge=c(weight_per_path[[i]])))) %>%
as.matrix
list(path_edge_table, path_admixedge_table, length(paths))
}
fill_pwts = function(pwts, weights, path_edge_table, path_admixedge_table, numpaths = NULL) {
# puts weights onto pwts, using index matrix and vectors
if(length(weights)==0) return(pwts)
wts2 = rep(weights, each=2)*c(1,-1) + (0:1)
path_weights = path_admixedge_table %>% as_tibble %>% mutate(w = wts2[admixedge]) %>%
group_by(path) %>% summarize(w = prod(w))
pwts_weights = path_edge_table %>% as_tibble %>% left_join(path_weights, by='path') %>%
mutate(w = ifelse(is.na(w), 1, w)) %>% group_by(edge2, leaf2) %>% summarize(w = sum(w)) %>% ungroup %>% as.matrix
pwts[pwts_weights[,1:2]] = pwts_weights[,3]
pwts
}
optimweightsfun = function(weights, args) {
# likelihood function used in optimizing admixture weights
# weights is vector of admixture weights to be optmized; only values for first incoming edge; 2nd is 1 - first
pwts = args[[1]]
ppinv = args[[2]]
f3_est = args[[3]]
path_edge_table = args[[4]] # indices into pwts with weight positions
path_admixedge_table = args[[5]]
#index3 = args[[6]]
cmb = args[[7]]
qpsolve = args[[8]]
lower = args[[9]]
upper = args[[10]]
fudge = args[[11]]
baseind = args[[12]]
constrained = args[[13]]
pwts = fill_pwts(pwts, weights, path_edge_table, path_admixedge_table)
pwts = pwts[,-baseind] - pwts[,baseind]
ppwts_2d = t(pwts[,cmb[1,]]*pwts[,cmb[2,]])
branch_lengths = opt_edge_lengths(ppwts_2d, ppinv, f3_est, qpsolve, lower, upper, fudge = fudge, constrained = constrained)
f3_fit = ppwts_2d %*% branch_lengths
qpgraph_score(f3_fit, f3_est, ppinv)
}
opt_edge_lengths = function(ppwts_2d, ppinv, f3_est, qpsolve, lower, upper, fudge = 1e-4, constrained = TRUE) {
# finds optimal edge lengths
# pwts2d: npair x nedge design matrix with paths to outpop
# ppinv: inverse of npair x npair covariance matrix of jackknife f3 stats
# f3_est: estimated f3 stats
pppp = t(ppwts_2d) %*% ppinv
cc = pppp %*% ppwts_2d
nc = ncol(cc)
diag(cc) = diag(cc) + fudge*mean(diag(cc))
sc = sqrt(diag(cc))
q1 = (pppp %*% f3_est)[,1]/sc
cc = cc/tcrossprod(sc)
if(constrained) qpsolve(cc, q1, cbind(diag(nc), -diag(nc)), c(lower*sc, -upper*sc))/sc
else solve(cc, q1)/sc
}
qpgraph_score = function(f3_fit, f3_est, ppinv = diag(length(f3_fit))) {
res = f3_est - f3_fit
lik = t(res) %*% ppinv %*% res
lik[1,1]
}
treemix_score = function(f3_fit, f3_est, ppinv) {
res = f3_est - f3_fit
se = sqrt(diag(solve(ppinv)))
sum(res^2/(2*se^2) + log(se * sqrt(2*pi)))
}
#' Compute the fit of an admixture graph
#'
#' Computes the fit of a given admixturegraph from f2-statistics. Drift edge weights and admixture edges weights are optimized until the (negative) likelihood score is minimized. The likelihood score is based on the squared difference between estimated and fitted f3-statistics.
#' @export
#' @param data Input data in one of three forms:
#' \enumerate{
#' \item A 3d array of blocked f2 statistics, output of \code{\link{f2_from_precomp}} or \code{\link{extract_f2}} (fastest option)
#' \item A directory which contains pre-computed f2-statistics
#' \item The prefix of genotype files (slowest option)
#' }
#' @param graph An admixture graph represented as a matrix of edges, an \code{\link{igraph}} object, or the path to a *qpGraph* graph file. Edges can be constrained by providing a matrix or data frame of edges with columns titled `lower` and `upper` with lower and upper bounds, respectively. By default, admixture edges are constrained to be between zero and one (with paired edges summing to one), and drift edges have a lower bound at zero.
#' @param lambdascale Scales f2-statistics. This has no effect on the fit, but is used in the original *qpGraph* program to display branch weights on a scale that corresponds to FST distances.
#' @param boot If `FALSE` (the default), each block will be left out at a time and the covariance matrix of
#' f3 statistics will be computed using block-jackknife. Otherwise bootstrap resampling is performed `n` times,
#' where `n` is either equal to `boot` if it is an integer, or equal to the number of blocks if `boot` is `TRUE`.
#' The covariance matrix of f3 statistics will be computed using bootstrap resampling.
#' @param diag Regularization term added to the diagonal elements of the covariance matrix of fitted branch lengths (after scaling by the matrix trace). Default is 0.0001.
#' @param diag_f3 Regularization term added to the diagonal elements of the covariance matrix of estimated f3 statistics (after scaling by the matrix trace). In the original *qpGraph* program, this is fixed at 0.00001.
#' @param lsqmode Least-squares mode. If `TRUE`, the likelihood score will be computed using a diagonal matrix with `1/(sum(diag(f3_var)) * diag_f3)`, in place of the inverse f3-statistic covariance matrix.
#'
#' `lsqmode = 2` will use the identity matrix instead, which is equivalent to computing the score as the sum of squared residuals (`sum((f3_est-f3_fit)^2)`).
#'
#' Both of these options do not take the covariance of f3-statistics into account. This can lead to bias, but is more stable in cases where the inverse f3-statistics covariance matrix can not be estimated precisely (for example because the number of populations is large). An alternative to `lsqmode = TRUE` that doesn't completely ignore the covariance of f3-statistics is to increase `diag_f3`.
#' @param numstart Number of random initializations of starting weights. Defaults to 10. Increasing this number will make the optimization slower, but reduce the risk of not finding the optimal weights. Check the `opt` output to see how much the optimization depends on the starting weights.
#' @param seed Random seed for generating starting weights.
#' @param cpp Use C++ functions. Setting this to `FALSE` will be slower but can help with debugging.
#' @param return_fstats Return estimated and fitted f2- and f4-statistics, as well as the worst f4-statistic residual Z-score. Defaults to `FALSE` because this can be slow.
#' @param f3precomp Optional precomputed f3-statistics. This should be the output of \code{\link{qpgraph_precompute_f3}} and can be provided instead of `data`. This can speed things up if many graphs are evaluated using the same set of f3-statistics.
#' @param f3basepop Optional f3-statistics base population. Inference will be based on f3-statistics of the form `f3(f3basepop; i, j)` for all population pairs `(i, j)`. Defaults to the outgroup population if the graph has one. This option is ignored if `f3precomp` is provided. Changing `f3basepop` should make very little difference.
#' @param constrained Constrain estimated drift edge weights to be non-negative, and admixture edge weights to be between zero and one.
#' @param allsnps Compute f3 from different SNPs for each population triplet (if data is missing for some SNPs and populations). This only has an effect when `data` is the prefix of genotype files.
#' @param ppinv Optional inverse f3-statistics covariance matrix
#' @param f2_blocks_test An optional 3d array of f2-statistics used for computing an out-of-sample score. This should contain only SNP blocks which are not part of `f2_blocks`. This allows to estimate the fit of a graph without overfitting and will not be used during the optimization step
#' @param verbose Print progress updates
#' @return `qpgraph` returns a list with data describing the model fit:
#' \itemize{
#' \item `edges`: A data frame where each row is an edge in the graph. For regular edges,
#' the column `weight` is the estimated edge length, and for admixture edges, it is the estimated admixture weight.
#' \item `score`: The likelihood score of the fitted graph. Lower values correspond to better fits.
#' The score is calculated as the inner product of the residuals (difference between estimated and
#' fitted f3 statistics), weighted by the inverse of the f3 covariance matrix. See \code{\link{qpgraph_score}}
#' \item `f2`: Estimated and fitted f2 statistics (if `return_fstats = TRUE`). p-values and z-scores test the significance of the difference.
#' \item `f3`: Estimated and fitted f3 statistics. p-values and z-scores test the significance of the difference.
#' \item `f4`: Estimated and fitted f4 statistics (if `return_fstats = TRUE`). p-values and z-scores test the significance of the difference.
#' \item `opt`: A data frame with details of the weight-fitting step, including the randomly sampled starting weights. The column `value` contains the score for each set of starting weights. Columns starting with `x` denote initial weights, and columns starting with `y` denote fitted weights.
#' \item `worst_residual`: The highest absolute z-score of f4-statistics residuals (fitted - estimated f4); (returned if `return_fstats = TRUE`)
#' }
#' @references Patterson, N. et al. (2012) \emph{Ancient admixture in human history.} Genetics
#' @seealso \code{\link{qpgraph_wrapper}} for a wrapper functions which calls the original *qpGraph* program.
#' @examples
#' out = qpgraph(example_f2_blocks, example_graph)
#' plot_graph(out$edges)
qpgraph = function(data, graph, lambdascale = 1, boot = FALSE, diag = 1e-4, diag_f3 = 1e-5,
lsqmode = FALSE, numstart = 10, seed = NULL, cpp = TRUE, return_fstats = FALSE,
return_pvalue = FALSE, f3precomp = NULL, f3basepop = NULL,
constrained = TRUE, allsnps = FALSE, ppinv = NULL, f2_blocks_test = NULL, verbose = FALSE) {
#----------------- process graph -----------------
if('matrix' %in% class(graph)) {
edges = as.data.frame(graph, stringsAsFactors = FALSE)
} else if('character' %in% class(graph)) {
edges = parse_qpgraph_graphfile(graph)
} else if('igraph' %in% class(graph)) {
edges = igraph::as_edgelist(graph) %>% as.data.frame(stringsAsFactors = FALSE)
} else if('data.frame' %in% class(graph)) {
edges = graph
} else stop(paste0('Cannot parse graph of class ', class(graph),'!'))
if(lambdascale == -1) lambdascale = 1
if(!lambdascale > 0) stop("'lambdascale' has to be > 0!")
if(!isFALSE(return_fstats) && is.null(data)) stop("Can't compute f4 without f2 data!")
if(cpp) {
optimweightsfun = cpp_optimweightsfun
opt_edge_lengths = cpp_opt_edge_lengths
fill_pwts = cpp_fill_pwts
get_pairindex = cpp_get_pairindex
}
if(class(graph)[1] != 'igraph') graph = graph_from_edgelist(as.matrix(edges[,1:2]))
nedges = length(E(graph))
admixnodes = which(degree(graph, mode='in') == 2)
nadmix = length(admixnodes)
admixedgesfull = sapply(seq_len(nadmix), function(i) incident_edges(graph, admixnodes, mode='in')[[i]][1:2])
normedges = setdiff(1:nedges, admixedgesfull)
pops = get_leafnames(graph)
npop = length(pops)
cmb = combn(0:(npop-1), 2)+(1:0)
if(!is.null(data) && !is.null(f3precomp)) stop("'f2_blocks' and 'f3precomp' can't both be provided!")
if(!is.null(f3precomp)) {
precomp = f3precomp
f3pops = attr(precomp$f3_est, 'pops')
pairmatch = get_pairindex(match(pops, f3pops))
precomp$f3_est = precomp$f3_est[pairmatch]
precomp$ppinv = precomp$ppinv[pairmatch, pairmatch]
precomp$f3out %<>% slice(pairmatch)
baseind = which(pops == f3pops[1])
} else if(is_geno_prefix(data)) {
if(!isFALSE(return_fstats)) stop("The 'return_fstats' option currently doesn't work when the first argument is the prefix of a genotype file!")
precomp = qpgraph_precompute_f3_geno(data, pops, f3basepop, allsnps = allsnps, outgroupmode = TRUE,
apply_corr = 2, poly_only = FALSE, adjust_pseudohaploid = FALSE)
#precomp = qpgraph_precompute_f3_geno(data, pops, f3basepop, ...)
baseind = if(is.null(f3basepop)) 1 else which(pops == f3basepop)
} else {
if(is.data.frame(data) || is.matrix(data)) {
# sets f3 covariance matrix to identity matrix
if(is.data.frame(data)) {
data %<>% select(1:3) %>% set_colnames(c('pop1', 'pop2', 'f2')) %>% filter(pop1 < pop2) %>%
bind_rows(rename(., pop1=pop2,pop2=pop1)) %>% bind_rows(tibble(pop1=pops,pop2=pops,f2=0)) %>%
arrange(pop1,pop2)
f2mat = data %>% pivot_wider(names_from=pop2, values_from=f2) %>% column_to_rownames('pop1') %>%
as.matrix
} else {
f2mat = data
data = f2mat %>% as_tibble(rownames = 'pop1') %>%
pivot_longer(-pop1, names_to = 'pop2', values_to = 'f2')
}
f2mat = f2mat[pops,pops]
precomp = list()
f3mat = (t(t(-f2mat + f2mat[,1])+f2mat[,1])/2)[-1,-1]
precomp$f3_est = c(f3mat[!upper.tri(f3mat)])
precomp$ppinv = diag(choose(length(pops), 2))
precomp$f3out = data %>% transmute(pop1,pop2,est=f2,se=1) %>% filter(pop1 < pop2)
} else {
f2_blocks = get_f2(data, pops, afprod = FALSE, verbose = verbose)
precomp = qpgraph_precompute_f3(f2_blocks, pops, f3basepop = f3basepop, lambdascale = lambdascale, boot = boot,
seed = seed, diag_f3 = diag_f3, lsqmode = lsqmode)
}
baseind = if(is.null(f3basepop)) 1 else which(pops == f3basepop)
}
stopifnot(all(!is.na(precomp$ppinv)))
if(!is.null(ppinv)) {
if(!is.null(f3precomp)) stop("'f3precomp' and 'ppinv' can't both be provided!")
f3pops = attr(ppinv, 'pops')
pairmatch = get_pairindex(match(pops, f3pops))
precomp$ppinv = ppinv[pairmatch, pairmatch]
}
f3_est = precomp$f3_est
#f3_est = pmax(0, f3_est)
ppinv = precomp$ppinv
weight = low = high = rep(NA, nedges)
pwts = graph_to_pwts(graph, pops)
opt = NULL
mim = .Machine$integer.max
if('lower' %in% names(edges)) {
elower = replace_na(edges$lower[normedges], 0)
eupper = replace_na(edges$upper[normedges], mim)
} else {
elower = rep(0, length(normedges))
eupper = rep(mim, length(normedges))
}
if(nadmix > 0) {
if(!is.null(seed)) set.seed(seed)
if('lower' %in% names(edges)) {
alower = replace_na(pmax(edges$lower[admixedgesfull[1,]], 1-edges$upper[admixedgesfull[2,]], na.rm=T), 0)
aupper = replace_na(pmin(edges$upper[admixedgesfull[1,]], 1-edges$lower[admixedgesfull[2,]], na.rm=T), 1)
aupper = pmin(1, aupper) + 1e-9
} else if(constrained) {
alower = rep(0, nadmix)
aupper = rep(1, nadmix)
} else {
alower = rep(-Inf, nadmix)
aupper = rep(Inf, nadmix)
}
parmat = matrix(runif(numstart*nadmix), numstart)
if(verbose) alert_info(paste0('Testing ', nrow(parmat), ' combinations of admixture weight starting values...\n'))
weightind = graph_to_weightind(graph)
arglist = list(pwts, ppinv, f3_est, weightind[[1]], weightind[[2]], weightind[[3]],
cmb, qpsolve, elower, eupper, diag, baseind, constrained)
oo = multistart(parmat, optimweightsfun, args = arglist, method = 'L-BFGS-B',
lower = alower, upper = aupper, control=list(maxit = 1e4, fnscale = 1e-6),
verbose = verbose)
best = oo %>% top_n(1, -value)
opt = data.frame(parmat, oo, stringsAsFactors = F)
admnames = names(V(graph))[admixnodes]
colnames(opt)[1:(nadmix*2)] = paste0(rep(c('i.', 'e.'), each = nadmix), rep(admnames, 2))
hilo = apply(as.matrix(oo[,1:nadmix]), 2, function(x) quantile(x, c(0, 1), na.rm = TRUE))
wts = as.matrix(best[,1:nadmix])[1,]
ta = c(t(admixedgesfull))
weight[ta] = c(wts, 1-wts)
low[ta] = c(pmin(hilo[1,], hilo[2,]), pmin(1-hilo[1,], 1-hilo[2,]))
high[ta] = c(pmax(hilo[1,], hilo[2,]), pmax(1-hilo[1,], 1-hilo[2,]))
pwts = fill_pwts(pwts, wts, weightind[[1]], weightind[[2]], weightind[[3]])
}
pwts = pwts[,-baseind] - pwts[,baseind]
ppwts_2d = t(pwts[,cmb[1,]]*pwts[,cmb[2,]])
branch_lengths = opt_edge_lengths(ppwts_2d, ppinv, f3_est, qpsolve, elower, eupper, fudge = diag, constrained = constrained)
f3_fit = ppwts_2d %*% branch_lengths
score = qpgraph_score(f3_fit, f3_est, ppinv)
if(!is.null(f2_blocks_test)) {
precomp_test = qpgraph_precompute_f3(f2_blocks_test, pops, lambdascale = lambdascale, boot = boot,
seed = seed, diag_f3 = diag_f3, lsqmode = lsqmode)
score_test = qpgraph_score(f3_fit, precomp_test$f3_est, ppinv)
}
# if(constrained) weight[normedges] = pmax(0, branch_lengths)
# else weight[normedges] = branch_lengths
weight[normedges] = branch_lengths
edges %<>% select(1:2) %>% set_colnames(c('from', 'to')) %>% as_tibble %>%
mutate(type = ifelse(1:n() %in% normedges, 'edge', 'admix'), weight = weight, low = low, high = high)
f2 = precomp$f2out
f3 = precomp$f3out %>% mutate(fit = c(f3_fit), diff = est - fit, z = diff/se, p = ztop(z))
out = namedList(edges, score, f3, opt, ppinv)
if(!is.null(f2_blocks_test)) {
out[['score_test']] = score_test
out[['f3_test']] = precomp_test$f3out
}
if(isTRUE(return_fstats) || return_fstats == 'f4') {
if(verbose) alert_info(paste0('Computing f4\n'))
out$f4 = fitf4(f2_blocks[pops, pops, ], f2, f3)
out$f2 = out$f4 %>% filter(pop1 == pop3, pop2 == pop4) %>% select(-pop3, -pop4)
out$worst_residual = max(abs(out$f4$z))
} else if(return_fstats == 'f2') {
out$worst_residual = fitf2(f2_blocks[pops, pops, ], f2, f3)$z %>% abs %>% max
} else if(return_fstats == 'f3') {
out$worst_residual = max(abs(f3$z))
}
if(return_pvalue) out$p.value = qpgraph_pvalue(f2_blocks, graph)
out
}
#' Compute f3-statistics from f2-statistics.
#'
#' Takes a 3d array of f2 block jackknife estimates and computes f3-statistics between the
#' first population \eqn{p1} and all population pairs \eqn{i, j}: \eqn{f3(p1; p_i, p_j)}
#' @export
#' @param data Input data in one of three forms:
#' \enumerate{
#' \item A 3d array of blocked f2 statistics, output of \code{\link{f2_from_precomp}} or \code{\link{extract_f2}} (fastest option)
#' \item A directory which contains pre-computed f2-statistics
#' \item The prefix of genotype files (slowest option)
#' }
#' @param pops Populations for which to compute f3-statistics
#' @param f3basepop f3-statistics base population. If `NULL` (the default),
#' the first population in `pops` will be used as the basis.
#' @param lambdascale Scales f2-statistics. This has no effect on the fit, but is used in the original *qpGraph* program to display branch weights on a scale that corresponds to FST distances.
#' @param boot If `FALSE` (the default), block-jackknife resampling will be used to compute standard errors.
#' Otherwise, block-bootstrap resampling will be used to compute standard errors. If `boot` is an integer, that number
#' will specify the number of bootstrap resamplings. If `boot = TRUE`, the number of bootstrap resamplings will be
#' equal to the number of SNP blocks.
#' @param seed Random seed used if `boot` is `TRUE`.
#' @param diag_f3 Regularization term added to the diagonal elements of the covariance matrix of estimated f3 statistics (after scaling by the matrix trace). In the original *qpGraph* program, this is fixed at 0.00001.
#' @param lsqmode Least-squares mode. If `TRUE`, the likelihood score will be computed using a diagonal matrix with `1/(sum(diag(f3_var)) * diag_f3)`, in place of the inverse f3-statistic covariance matrix. `lsqmode = 2` will use the identity matrix instead, which is equivalent to computing the score as the sum of squared residuals (`sum((f3_est-f3_fit)^2)`). Both of these options do not take the covariance of f3-statistics into account. This can lead to bias, but is more stable in cases where the inverse f3-statistics covariance matrix can not be estimated precisely (for example because the number of populations is large). An alternative to using `lsqmode = TRUE` which doesn't completely ignore the covariance of f3-statistics is to increase `diag_f3`.
#' @return A list with four items
#' \enumerate{
#' \item `f3_est` a matrix with f3-statistics for all population pairs with the output
#' \item `ppinv` a matrix with the inverse of the f3-statistic covariance matrix
#' \item `f2out` a data frame with f2 estimates
#' \item `f3out` a data frame with f3 estimates
#' }
#' @examples
#' pops = get_leafnames(example_igraph)
#' qpgraph_precompute_f3(example_f2_blocks, pops)$f3_est
#' \dontrun{
#' qpgraph_precompute_f3(f2_dir, pops)
#' }
qpgraph_precompute_f3 = function(data, pops, f3basepop = NULL, lambdascale = 1, boot = FALSE,
seed = NULL, diag_f3 = 1e-5, lsqmode = FALSE) {
# returns list of f3_est and ppinv for subset of populations.
# f3_est and ppinv are required for qpgraph_slim; f2out and f3out are extra output
# f2_blocks may contain more populations than the ones used in qpgraph
# f2_blocks input here should be subset which is used by qpgraph function
#----------------- read f-stats -----------------
if(!is.null(f3basepop)) pops = c(f3basepop, setdiff(pops, f3basepop))
if(!is.null(seed)) set.seed(seed)
samplefun = ifelse(boot, function(x) est_to_boo(x, boot), est_to_loo)
matstatfun = ifelse(boot, boot_mat_stats, jack_mat_stats)
arrstatfun = ifelse(boot, boot_arr_stats, jack_arr_stats)
f2_blocks = (get_f2(data, pops) * lambdascale) %>% samplefun
#f2_blocks = array(pmax(0, f2_blocks), dim(f2_blocks), dimnames(f2_blocks))
block_lengths = parse_number(dimnames(f2_blocks)[[3]])
npop = length(pops)
npair = choose(npop, 2)
cmb = combn(0:(npop-1), 2)+(1:0)
f2 = arrstatfun(f2_blocks, block_lengths)
f2out = tibble(pop1 = combn(pops, 2)[1,],
pop2 = combn(pops, 2)[2,],
est = f2[[1]][lower.tri(f2[[1]])],
se = sqrt(f2[[2]][lower.tri(f2[[2]])]))
f3_blocks = (f2_blocks[,rep(1, npop),,drop=F] + f2_blocks[rep(1, npop),,,drop=F] - f2_blocks)/2
#f3_blocks = array(pmax(0, f3_blocks), dim(f3_blocks), dimnames(f3_blocks))
f3_blocks_2d = arr3d_to_pairmat(f3_blocks[-1,-1,,drop=F])
f3dat = matstatfun(f3_blocks_2d, block_lengths)
#f3dat = jack_mat_stats(f3_blocks_2d, block_lengths)
f3_est = f3dat$est
f3_var = f3dat$var
f3out = tibble(pop1 = pops[1],
pop2 = pops[cmb[1,]+1],
pop3 = pops[cmb[2,]+1],
est = f3_est, se = sqrt(diag(f3_var)))
add_diag = sum(diag(f3_var)) * diag_f3
diag(f3_var) = diag(f3_var) + add_diag
# in qpGraph diag_f3 is 1e-5; has large effect on magnitude of likelihood score
if(lsqmode == 2) ppinv = diag(nrow(f3_var))
else if(lsqmode) ppinv = diag(nrow(f3_var)) / add_diag
else ppinv = solve(f3_var)
f3_est %<>% structure(pops = pops)
ppinv %<>% structure(pops = pops)
namedList(f3_est, ppinv, f2out, f3out, f3_blocks_2d)
}
qpgraph_precompute_f3_geno = function(data, pops, f3basepop = NULL, lambdascale = 1, boot = FALSE,
seed = NULL, diag_f3 = 1e-5, lsqmode = FALSE, outgroupmode = TRUE, adjust_pseudohaploid = FALSE,
poly_only = FALSE, apply_corr = TRUE, allsnps = FALSE, verbose = TRUE) {
if(!is.null(f3basepop)) pops = c(f3basepop, setdiff(pops, f3basepop))
f3basepop = pops[1]
popcombs = expand_grid(pop1 = pops[1], pop2 = pops[-1], pop3 = pops[-1]) %>%
filter(pop2 <= pop3)
f3blockdat = f3blockdat_from_geno(data, popcombs, outgroupmode = outgroupmode, apply_corr = apply_corr,
poly_only = poly_only, adjust_pseudohaploid = adjust_pseudohaploid,
allsnps = allsnps, verbose = verbose)
f3blocks = f3blockdat %>% select(-pop1) %>% rename(pop1 = pop2, pop2 = pop3, f2 = est) %>%
f2dat_to_f2blocks(fill_diag = FALSE)
#f3blocks = f3blocks[pops[-1],pops[-1],-apply(f3blocks, 3, function(x) any(is.na(x)))]
f3blocks = f3blocks[pops[-1],pops[-1],]
if(!is.null(seed)) set.seed(seed)
samplefun = ifelse(boot, function(x) est_to_boo(x, boot), est_to_loo)
matstatfun = ifelse(boot, boot_mat_stats, jack_mat_stats)
arrstatfun = ifelse(boot, boot_arr_stats, jack_arr_stats)
f3blocks = f3blocks %>% samplefun
block_lengths = parse_number(dimnames(f3blocks)[[3]])
npop = length(pops)
npair = choose(npop, 2)
cmb = combn(0:(npop-1), 2)+(1:0)
f3 = arrstatfun(f3blocks, block_lengths)
f3_blocks_2d = arr3d_to_pairmat(f3blocks)
f3dat = matstatfun(f3_blocks_2d, block_lengths)
f3_est = f3dat$est
f3_var = f3dat$var
f3out = tibble(pop1 = pops[1],
pop2 = pops[cmb[1,]+1],
pop3 = pops[cmb[2,]+1],
est = f3_est, se = sqrt(diag(f3_var)))
# est = f3[[1]][lower.tri(f3[[1]],T)],
# se = sqrt(f3[[2]][lower.tri(f3[[2]],T)]))
add_diag = sum(diag(f3_var)) * diag_f3
diag(f3_var) = diag(f3_var) + add_diag
# in qpGraph diag_f3 is 1e-5; has large effect on magnitude of likelihood score
if(lsqmode == 2) ppinv = diag(nrow(f3_var))
else if(lsqmode) ppinv = diag(nrow(f3_var)) / add_diag
else ppinv = solve(f3_var)
f3_est %<>% structure(pops = pops)
ppinv %<>% structure(pops = pops)
namedList(f3_est, ppinv, f3out, f3_blocks_2d)
}
get_pairindex = function(perm) {
# returns index vector that matches population pairs in
# f3_est and ppinv (which were computed using pops) to graph populations
cmb = combn(0:(length(perm)-1), 2)+(1:0)
popind = setdiff(perm, 1)
orig_order = apply(cmb+1, 2, paste0, collapse='')
new_order = apply(matrix(popind[c(cmb)], 2), 2, function(x) paste0(sort(x), collapse=''))
match(new_order, orig_order)
}
qpgraph_anorexic = function(f3precomp, graph, diag = 1e-4,
numstart = 10, seed = NULL, verbose = FALSE, cpp = TRUE) {
# only works for trees at the moment, because weightind order is coupled to pwts order
admixnodes = which(degree(graph, mode='in') == 2)
stopifnot(length(admixnodes) == 0)
if(cpp) {
opt_edge_lengths = cpp_opt_edge_lengths
}
graphpops = get_leafnames(graph)
f3_est = f3precomp$f3_est
ppinv = f3precomp$ppinv
f3pops = attr(f3_est, 'pops')
pwts = tree_to_pwts(graph, graphpops)
pwts = pwts[, match(f3pops[-1], setdiff(graphpops, f3pops[1]))]
cmb = combn(0:(length(f3pops)-1), 2)+(1:0)
ppwts_2d = t(pwts[,cmb[1,]] * pwts[,cmb[2,]])
branch_lengths = opt_edge_lengths(ppwts_2d, ppinv, f3_est, qpsolve,
lower = rep(0, nrow(pwts)), upper = rep(.Machine$integer.max, nrow(pwts)),
fudge = diag, constrained = TRUE)
f3_fit = ppwts_2d %*% branch_lengths
score = qpgraph_score(f3_fit, f3_est, ppinv)
edges = as_tibble(as_edgelist(graph), .name_repair = ~c('from', 'to')) %>%
mutate(type = 'edge', weight = c(branch_lengths))
namedList(edges, score, opt = NULL)
}
# not used
f3out_to_fittedf2out = function(f2out, f3out) {
# computes fitted f2 statistics data frame from f2 and f3 data frames
# will not include f2(outgroup, X)
f2out %>%
right_join(f3out %>% filter(pop2 != pop3) %>% transmute(pop1=pop2, pop2=pop3, f3 = fit), by = c('pop1', 'pop2')) %>%
left_join(f3out %>% filter(pop2 == pop3) %>% transmute(pop1=pop2, f2_1 = fit), by = c('pop1')) %>%
left_join(f3out %>% filter(pop2 == pop3) %>% transmute(pop2=pop3, f2_2 = fit), by = c('pop2')) %>%
transmute(pop1, pop2, est, se, fit = (f2_1 + f2_2 - f3*2), diff = est - fit, z = diff/se, p = ztop(z))
}
fitf2 = function(f2_blocks, f2, f3) {
# returns a tibble with estimated and fitted f2-statistics
cmb = combn(0:(dim(f2_blocks)[1]-1), 2)+(1:0)
f2_out = f3 %>% filter(pop2 == pop3) %$% fit
f2_fit = f3 %>% mutate(f21 = f2_out[cmb[1,]], f22 = f2_out[cmb[2,]], f2fit = (f21 + f22 - fit*2))
f2_fit2 = f2 %>%
left_join(f2_fit, by = c('pop1'='pop2', 'pop2'='pop3')) %>%
filter(!is.na(f2fit)) %>%
transmute(pop1, pop2, fit = f2fit) %>%
bind_rows(f2_fit %>% filter(pop2 == pop3) %>% select(pop1, pop2, fit))
f2(f2_blocks, verbose = FALSE) %>%
inner_join(f2_fit2, by = c('pop1', 'pop2')) %>%
mutate(diff = est - fit, z = diff/se, p = ztop(z))
}
fitf4 = function(f2_blocks, f2, f3) {
# returns a tibble with estimated and fitted f4-statistics
cmb = combn(0:(dim(f2_blocks)[1]-1), 2)+(1:0)
f2_out = f3 %>% filter(pop2 == pop3) %$% fit
f2_fit = f3 %>% mutate(f21 = f2_out[cmb[1,]], f22 = f2_out[cmb[2,]], f2fit = (f21 + f22 - fit*2))
f2_fit2 = f2 %>%
left_join(f2_fit, by = c('pop1'='pop2', 'pop2'='pop3')) %>%
filter(!is.na(f2fit)) %>%
select(pop1, pop2, f2fit) %>%
bind_rows(f2_fit %>% filter(pop2 == pop3) %>% transmute(pop1, pop2, f2fit = fit)) %>%
bind_rows(rename(., pop1 = pop2, pop2 = pop1)) %>%
bind_rows(tibble(pop1 = unique(.$pop1), pop2 = pop1, f2fit = 0))
x = f4(f2_blocks, unique_only = FALSE, verbose = FALSE) %>% select(-z, -p)
x %>%
left_join(f2_fit2 %>% rename(c1 = f2fit), by = c('pop1' = 'pop1', 'pop4' = 'pop2')) %>%
left_join(f2_fit2 %>% rename(c2 = f2fit), by = c('pop2' = 'pop1', 'pop3' = 'pop2')) %>%
left_join(f2_fit2 %>% rename(c3 = f2fit), by = c('pop1' = 'pop1', 'pop3' = 'pop2')) %>%
left_join(f2_fit2 %>% rename(c4 = f2fit), by = c('pop2' = 'pop1', 'pop4' = 'pop2')) %>%
mutate(fit = (c1 + c2 - c3 - c4)/2, diff = est - fit, z = diff/se, p = ztop(z)) %>%
select(-c1:-c4)
}
#' Compare the fit of two qpgraph models
#'
#' Takes two data frames with model fits computed on two graphs for on the same populations and tests whether the scores of one graph are significantly better than the scores of the other
#' @param fits1 The fits of the first graph
#' @param fits2 The fits of the second graph
#' @param boot should match the `boot` parameter in `qpgraph_resample_snps` (`FALSE` by default)
#' @examples
#' \dontrun{
#' nblocks = dim(example_f2_blocks)[3]
#' train = sample(1:nblocks, round(nblocks/2))
#' fits1 = qpgraph_resample_snps(example_f2_blocks[,,train], graph = graph1,
#' f2_blocks_test = example_f2_blocks[,,-train])
#' fits2 = qpgraph_resample_snps(example_f2_blocks[,,train], graph = graph2,
#' f2_blocks_test = example_f2_blocks[,,-train])
#' compare_fits2(fit1, fit2)
#' }
compare_fits2 = function(fits1, fits2, boot = FALSE) {
matstatfun = ifelse(boot, boot_mat_stats, jack_mat_stats)
stats = matstatfun(t(fits1$score_test - fits2$score_test), rep(1, length(fits1$score_test)))
diff = stats$est
se = sqrt(stats$var)
z = diff/se
p = ztop(z)
namedList(diff, se, z, p, scores1 = fits1$score_test, scores2 = fits2$score_test, boot)
}
#' Compare the fit of two qpgraph models
#'
#' Takes the bootstrap score distribution of two fits on the same populations and tests whether the scores of one graph are significantly higher or lower than the scores of the other graph.
#' @export
#' @param scores1 Scores for the first graph. Each score should be for same model evaluated on different bootstrap samples of the SNP blocks. See \code{\link{qpgraph_resample_multi}}
#' @param scores2 Scores for the second graph, evaluated on the same bootstrap samples as the first graph
#' @return A list with statistics comparing the two models
#' \itemize{
#' \item `p_emp`: The two-sided bootstrap p-value testing whether the scores of one model are higher or lower than the scores of the other model. It is two times the fraction of bootstrap replicates in which model 1 has a lower score than model 2 (or vice-vera, whichever is lower). This is simply a bootstrap test for testing whether some quantity (the the score difference of two models in this case) is significantly different from zero.
#' \item `p_emp_nocorr`: `p_emp` is truncated to be no less than the reciprocal of the number of bootstrap iterations (the length of the score vectors). `p_emp_nocorr` is not truncated and can be equal to 0.
#' \item `ci_low`: The 2.5% quantile of distribution of score differences
#' \item `ci_high`: The 97.5% quantile of distribution of score differences
#' \item The other items in this list are less important and can be ignored. In particular, `p` is not as reliable as `p_emp` because it assumes that the score differences follow a normal distribution.
#' }
#' @seealso \code{\link{qpgraph_resample_multi}}
#' @examples
#' \dontrun{
#' fits = qpgraph_resample_multi(f2_blocks, list(graph1, graph2), nboot = 100)
#' compare_fits(fits[[1]]$score, fits[[2]]$score)
#' }
compare_fits = function(scores1, scores2) {
if(length(scores1) != length(scores2)) stop('Length of scores1 and scores2 differ!')
scorediff = scores1 - scores2
ci_low = unname(quantile(scorediff, 0.025, na.rm = T))
ci_high = unname(quantile(scorediff, 0.975, na.rm = T))
scorediff = na.omit(scores1 - scores2)
stats = boot_mat_stats(t(scorediff), rep(1, length(scorediff)))
diff = stats$est
se = sqrt(stats$var)
z = diff/se
p = ztop(z)[1]
frac1 = mean(scorediff < 0)
frac2 = mean(scorediff > 0)
p_emp_nocorr = min(frac1, frac2)*2
p_emp = max(p_emp_nocorr, 1/length(scorediff))
namedList(diff, se, z, p, p_emp, p_emp_nocorr, ci_low, ci_high)
}
#' Evaluate a qpgraph model many times
#'
#' This function is used in combination with \code{\link{compare_fits}} in order to test whether one graph has a significantly better fit than another. using a list of bootstrap resampled f2-statistics and corresponding test-set f2-statistics
#' @export
#' @param f2_blocks A list of bootstrap resampled f2-statistics
#' @param graph An admixture graph
#' @param f2_blocks_test A list of f2-statistics from all blocks which were not used in the corresponding f2-array in `f2_blocks`
#' @param verbose Print progress updates
#' @param ... Arguments passed to \code{\link{qpgraph}}
#' @return A data frame with \code{\link{qpgraph}} results for each iteration of bootstrap resampled f2-statistics
#' @seealso \code{\link{compare_fits}} \code{\link{boo_list}}
qpgraph_resample_snps2 = function(f2_blocks, graph, f2_blocks_test, verbose = TRUE, ...) {
#progressr::with_progress({
#pb = progressr::progressor(steps = length(f2_blocks))
ell = list(...)
fun = function(f2dat, f2dat_test, g) function() {
#pb$tick()
#pb()
safely(qpgraph)(data = f2dat, graph = g, f2_blocks_test = f2dat_test, verbose = FALSE, ...)
}
#pb = progress::progress_bar$new(total = length(f2_blocks))
tibble(id = seq_len(length(f2_blocks)), graph = list(graph), f2_blocks, f2_blocks_test) %>%
mutate(fun2 = pmap(list(f2_blocks, f2_blocks_test, graph), fun)) %>%
mutate(out = furrr::future_invoke_map(fun2, .progress = verbose, .options = furrr::furrr_options(seed = TRUE)),
result = map(out, 'result', .null = tibble()), error = map(out, 'error')) %>%
select(-out, -fun2) %>% unnest_wider(result)
#})
}
# qpgraph_resample_snps2 = function(f2_blocks, graph, f2_blocks_test, verbose = TRUE, ...) {
#
# fun = function(f2dat, f2dat_test, g) function() safely(qpgraph)(data = f2dat, graph = g, f2_blocks_test = f2dat_test, verbose = FALSE, ...)
#
# tibble(id = seq_len(length(f2_blocks)), graph = list(graph), f2_blocks, f2_blocks_test) %>%
# mutate(fun2 = pmap(list(f2_blocks, f2_blocks_test, graph), fun)) %>%
# mutate(out = furrr::future_invoke_map(fun2, .progress = verbose),
# result = map(out, 'result', .null = tibble()), error = map(out, 'error')) %>%
# select(-out, -fun2) %>% unnest_wider(result)
# }
#' Evaluate a qpgraph models many times
#'
#' This function is used in combination with \code{\link{compare_fits}} in order to test whether one graph has a significantly better fit than another. It creates bootstrap resampled SNP block training and test sets, and uses them to evaluate multiple graphs.
#' @export
#' @param f2_blocks 3d array of f2-statistics
#' @param graphlist A list of admixture graphs
#' @param nboot Number of bootstrap iterations
#' @param verbose Print progress updates
#' @param ... Arguments passed to \code{\link{qpgraph}}
#' @return A list of same length as `graphlist` with data frames with \code{\link{qpgraph}} results for each iteration of bootstrap resampled f2-statistics
#' @examples
#' \dontrun{
#' fits = qpgraph_resample_multi(f2_blocks, list(graph1, graph2), nboot = 100)
#' compare_fits(fits[[1]]$score_test, fits[[2]]$score_test)
#' }
#' @seealso \code{\link{compare_fits}}
qpgraph_resample_multi = function(f2_blocks, graphlist, nboot, verbose = TRUE, ...) {
argnam = names(formals(qpgraph))
if(!all(...names() %in% argnam)) {
notused = setdiff(...names(), argnam)
stop(paste0("The following arguments are not recognized: '", paste0(notused, collapse = "', '"), "'"))
}
boo = boo_list(f2_blocks, nboot = nboot)
#f3pre = map(graphlist, ~qpgraph_precompute_f3(f2_blocks, get_leafnames(.))$ppinv)
#map2(graphlist, f3pre, function(.x, .y, ...) qpgraph_resample_snps2(
# boo$boo, .x, boo$test, ppinv = .y, verbose = verbose, ...), ...)
map(graphlist, function(.x, ...) qpgraph_resample_snps2(
boo$boo, .x, boo$test, verbose = verbose, ...), ...)
}
qpgraph_pvalue = function(f2_blocks, graph, nboot = 100, ...) {
boo = boo_list(f2_blocks, nboot = nboot)
score_fit = qpgraph_resample_snps2(boo$boo, graph, boo$test, numstart = 1, ...)$score_test
pops = get_leafnames(graph)
precomp_train = map(boo$boo, ~qpgraph_precompute_f3(., pops))
precomp_test = map(boo$test, ~qpgraph_precompute_f3(., pops))
score_emp = map2_dbl(precomp_train, precomp_test, ~qpgraph_score(.x$f3_est, .y$f3_est, .x$ppinv))
out = tibble(rep = 1:nboot, score_fit, score_emp) %>% mutate(diff = score_fit - score_emp)
t.test(out$diff, alternative = 'greater')$p.value
}
#' Compare the fit of two qpgraph models
#'
#' Takes two data frames with model fits computed on two graphs for on the same populations and tests whether the scores of one graph are significantly better than the scores of the other.
#' @param fit1 The fit of the first graph
#' @param fit2 The fit of the second graph
#' @param f2_blocks f2 blocks used for fitting `fit1` and `fit2`. Used in combination with `f2_blocks_test` to compute f-statistics covariance matrix.
#' @param f2_blocks_test f2 blocks which were not used for fitting `fit1` and `fit2`
#' @param boot If `FALSE` (the default), block-jackknife resampling will be used to compute standard errors.
#' Otherwise, block-bootstrap resampling will be used to compute standard errors. If `boot` is an integer, that number
#' will specify the number of bootstrap resamplings. If `boot = TRUE`, the number of bootstrap resamplings will be
#' equal to the number of SNP blocks. If bootstrap resampling is enabled, empirical p-values (`p_emp`) and 95 confidence intervals (`ci_low` and `ci_high`) will be reported.
#' @param seed Random seed used if `boot` is `TRUE`. Does not need to match a seed used in fitting the models
#' @examples
#' \dontrun{
#' nblocks = dim(example_f2_blocks)[3]
#' train = sample(1:nblocks, round(nblocks/2))
#' fit1 = qpgraph(example_f2_blocks[,,train], graph1)
#' fit2 = qpgraph(example_f2_blocks[,,train], graph2)
#' compare_fits4(fit1, fit2, example_f2_blocks[,,train], example_f2_blocks[,,-train])
#' }
compare_fits4 = function(fit1, fit2, f2_blocks, f2_blocks_test, boot = FALSE, seed = NULL) {
stopifnot(all.equal(sort(attr(fit1$ppinv, 'pops')), sort(attr(fit2$ppinv, 'pops'))))
matstatfun = ifelse(boot, boot_mat_stats, jack_mat_stats)
pops = attr(fit1$ppinv, 'pops')
ppinv = qpgraph_precompute_f3(abind::abind(f2_blocks, f2_blocks_test), pops, boot = boot, seed = seed)$ppinv
f3_test = qpgraph_precompute_f3(f2_blocks_test, pops, boot = boot, seed = seed)$f3_blocks_2d
f3_fit = fit1$f3 %>%
left_join(fit2$f3 %>% bind_rows(rename(., pop2=pop3, pop3=pop2) %>% filter(pop2 != pop3)),
by = c('pop1', 'pop2', 'pop3'))
scores1 = map_dbl(1:dim(f2_blocks_test)[3], ~qpgraph_score(f3_fit$fit.x, f3_test[,.], ppinv))
scores2 = map_dbl(1:dim(f2_blocks_test)[3], ~qpgraph_score(f3_fit$fit.y, f3_test[,.], ppinv))
scorediff = na.omit(scores1 - scores2)
stats = matstatfun(t(scorediff), rep(1, length(scorediff)))
diff = stats$est
se = sqrt(stats$var)
z = diff/se
p = ztop(z)
frac = mean(scorediff < 0)
p_emp = ci_low = ci_high = NA
if(boot) {
p_emp = min(frac, 1-frac)*2
ci_low = unname(quantile(scorediff, 0.025, na.rm = T))
ci_high = unname(quantile(scorediff, 0.975, na.rm = T))
}
namedList(diff, se, z, p, p_emp, ci_low, ci_high, scores1, scores2)
}
fit_to_worstf3z = function(fit) {
fit$f3 %>% left_join(fit$f3_test, by = c('pop1', 'pop2', 'pop3')) %>%
mutate(z = (est.y-fit)/sqrt(se.x^2+se.y^2)) %>% slice_max(abs(z)) %>% pluck('z', 1)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.