R/pwfst.R

In popkin: Estimate Kinship and FST under Arbitrary Population Structure

#' 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)
}