#' Estimate the individual-level pairwise FST matrix
#'
#' This function construct the individual-level pairwise FST matrix implied by the input kinship matrix.
#' If the input is the true kinship matrix, the return value corresponds to the true pairwise FST matrix.
#' On the other hand, if the input is the estimated kinship returned by [popkin()], the same code results in the return value being the pairwise FST estimates described in our paper.
#' In all cases the diagonal of the pairwise FST matrix is zero by definition.
#'
#' @param kinship The `n`-by-`n` kinship matrix
#'
#' @return The `n`-by-`n` pairwise FST matrix
#'
#' @examples
#' # Construct toy data
#' X <- matrix(c(0,1,2,1,0,1,1,0,2), nrow=3, byrow=TRUE) # genotype matrix
#' subpops <- c(1,1,2) # subpopulation assignments for individuals
#'
#' # NOTE: for BED-formatted input, use BEDMatrix!
#' # "file" is path to BED file (excluding .bed extension)
#' ## library(BEDMatrix)
#' ## X <- BEDMatrix(file) # load genotype matrix object
#'
#' # estimate the kinship matrix from the genotypes "X"!
#' kinship <- popkin(X, subpops) # calculate kinship from X and optional subpop labels
#'
#' # lastly, compute pairwise FST matrix from the kinship matrix
#' pwF <- pwfst(kinship)
#'
#' @export
pwfst <- function(kinship) {
# die if this is missing
if (missing(kinship))
stop('`kinship` matrix is required!')
# run additional validations
validate_kinship(kinship)
# the below code works best with kinship scaled like a coancestry matrix
kinship <- inbr_diag(kinship)
# extract inbreeding coefficients
inbrs <- diag(kinship)
# construct other things for estimation
n <- length(inbrs)
inbrMat <- matrix(inbrs, n, n) # repeats inbreeding coefficients along columns
# return the following desired pairwise Fst matrix!
( ( inbrMat + t(inbrMat) ) / 2 - kinship ) / (1 - kinship)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.