R/read_gcta_hsq.R

In OchoaLab/genbin: R wrappers for binaries in genetics

#' Read GCTA REML variance components table
#'
#' Parses the output generated by GCTA with option `--reml` (heritability estimation).
#' Returned table has column names standardized for convenience.
#' Validated with GCTA version 1.93.2 beta.
#'
#' @param file The file path to read, possibly missing the extension.
#' @param ext The expected file extension.
#' Set to `NA` to prevent adding any extension to `file`.
#' The default extension "hsq" agrees with the output of [gcta_reml()].
#' @param verbose If `TRUE` (default) function reports the path of the file being loaded (after autocompleting the extensions).
#'
#' @return A list containing these named elements:
#' - `herit`: the estimated heritability.
#' - `se`: the standard error of the heritability estimate.
#' - `data`: the full table of variance component estimates and other miscelaneous info.
#' 
#' @examples 
#' \dontrun{
#' obj <- read_gcta_hsq( file )
#' obj$herit
#' }
#'
#' @seealso
#' [gcta_reml()] wrapper for executing GCTA REML and reading the results.
#' 
#' @export
read_gcta_hsq <- function( file, ext = 'hsq', verbose = TRUE ) {
    if ( missing( file ) )
        stop( '`file` is required!' )

    # add .ext and/or .gz if missing and needed
    file <- genio:::add_ext_read(file, ext)
    
    # announce it if needed
    if (verbose)
        message('Reading: ', file)

    # read table!
    # unfortunately an irregular table
    # this always complains because last few rows only have 2 columns
    # issue avoided by only reading complete rows (will not read trailing data at all, we don't need it anyway)
    data <- readr::read_tsv(
                       file,
                       col_types = 'cdd',
                       n_max = 4 # not including header
                   )
    # find herit estimate
    # row of interest
    index <- data$Source == 'V(G)/Vp'
    # heritability point estimate
    herit <- data$Variance[ index ]
    # standard error from MVN model
    se <- data$SE[ index ]

    # return this
    return(
        list(
            herit = herit,
            se = se,
            data = data # full table
        )
    )
}