Read in a binary GRM from GCTA as square, symmetric matrix into R

read_grm_bin.R

## Read in a binary GRM from GCTA as a matrix
## see https://cnsgenomics.com/software/gcta/#MakingaGRM

library(Matrix)

read_grm_bin <- function(prefix, size=4) {

        bin_filename <- paste(prefix, 'grm.bin', sep='.')
        n_filename <- paste(prefix, 'grm.N.bin', sep='.')
        id_filename <- paste(prefix, 'grm.id', sep='.')

        id <- read.table(id_filename, col.names=c('FID', 'IID'))

        n <- nrow(id)

        bin_file <- file(bin_filename, 'rb')
        # binary GRM is stored as the upper-triangle of a dense, symmetric, packed matrix 
        GRM <- new("dspMatrix", uplo='U', x=readBin(bin_file, n=n*(n+1)/2, what=numeric(0), size=size), Dim=c(n, n))
        close(bin_file)
        
        n_file <- file(n_filename, 'rb')
        N <- new("dspMatrix", uplo='U', x=readBin(n_file, n=n*(n+1)/2, what=numeric(0), size=size), Dim=c(n, n))
        close(n_file)

        return(list(id=id, n=n, GRM=GRM, N=N))

}

## Example
## $ gcta64 --bfile test --make-grm --out test
## $ R
## > test <- read_grm_bin('test')
## > test$GRM[1:5,1:5]
## 5 x 5 Matrix of class "dspMatrix"
##             [,1]          [,2]          [,3]        [,4]        [,5]
## [1,]  0.97306669 -0.0426406004 -0.0108031128  0.04988819  0.01807305
## [2,] -0.04264060  1.0836025476 -0.0007884779 -0.05592454  0.04774712
## [3,] -0.01080311 -0.0007884779  1.1463223696  0.11909716  0.12911978
## [4,]  0.04988819 -0.0559245422  0.1190971583  1.06905055 -0.01259598
## [5,]  0.01807305  0.0477471203  0.1291197836 -0.01259598  0.93640965
## > head(test$id)
##   FID IID
## 1   1  11
## 2   2  21
## 3   3  31
## 4   4  41
## 5   5  51
## 6   6  61
## > test$N[1:5,1:5]
## 5 x 5 Matrix of class "dspMatrix"
##      [,1] [,2] [,3] [,4] [,5]
## [1,] 1000  994  992 1000 1000
## [2,]  994  994  986  994  994
## [3,]  992  986  992  992  992
## [4,] 1000  994  992 1000 1000
## [5,] 1000  994  992 1000 1000