基于R语言构建矩阵并求逆，导入GCTA亲缘关系矩阵

Got it, let's break down how to work with your GCTA-generated Genetic Relationship Matrix (GRM) in R—covering proper import, reconstructing the full symmetric matrix, and safely computing its inverse.

Step 1: Complete the ReadGRMBin Import Function #

Your initial function was cut off, so here's the full, working version to read all GCTA GRM binary files correctly:

ReadGRMBin <- function(prefix, AllN = FALSE, size = 4) {
  # Helper to calculate elements in lower triangle (including diagonal)
  sum_i <- function(i) sum(1:i)
  
  # Define paths to GRM files
  BinFileName <- paste(prefix, ".grm.bin", sep = "")
  NFileName <- paste(prefix, ".grm.N.bin", sep = "")
  IDFileName <- paste(prefix, ".grm.id", sep = "")
  
  # Read sample IDs
  id <- read.table(IDFileName, header = FALSE, stringsAsFactors = FALSE)
  n <- nrow(id)
  
  # Read flattened lower triangle GRM data
  BinFile <- file(BinFileName, "rb")
  grm_vec <- readBin(BinFile, n = sum_i(n), what = numeric(0), size = size)
  close(BinFile)
  
  # Read pairwise SNP count (N) values
  NFile <- file(NFileName, "rb")
  if (AllN) {
    # Read all pairwise N values (matches GRM's lower triangle structure)
    N_vec <- readBin(NFile, n = sum_i(n), what = numeric(0), size = size)
  } else {
    # Only read per-sample diagonal N values
    N_vec <- readBin(NFile, n = n, what = numeric(0), size = size)
  }
  close(NFile)
  
  # Return all imported data as a list
  list(id = id, grm = grm_vec, N = N_vec)
}

Step 2: Reconstruct the Full Symmetric GRM Matrix #

The imported grm_vec is a flattened lower triangle (including the diagonal). We need to convert this into a full symmetric matrix for downstream use:

reconstruct_grm <- function(grm_list) {
  n <- nrow(grm_list$id)
  grm_mat <- matrix(0, nrow = n, ncol = n)
  
  # Fill lower triangle and diagonal with GRM values
  grm_mat[lower.tri(grm_mat, diag = TRUE)] <- grm_list$grm
  
  # Mirror values to the upper triangle to make the matrix symmetric
  grm_mat[upper.tri(grm_mat)] <- t(grm_mat)[upper.tri(grm_mat)]
  
  # Assign sample IDs as row/column names for clarity
  rownames(grm_mat) <- paste(grm_list$id[,1], grm_list$id[,2], sep = "_")
  colnames(grm_mat) <- rownames(grm_mat)
  
  return(grm_mat)
}

# Example usage:
grm_data <- ReadGRMBin(prefix = "your_grm_prefix")  # Replace with your file prefix
full_grm <- reconstruct_grm(grm_data)

Step 3: Compute the GRM Inverse #

GRMs are often near-singular (especially with related samples), so direct inversion might fail. Here are two safe approaches:

Option 1: Direct Inversion (for non-singular matrices) #

If your GRM is well-conditioned, use base R's solve():

grm_inv <- solve(full_grm)

Option 2: Regularized Inversion (safer for near-singular GRMs) #

Add a tiny regularization term to the diagonal to stabilize inversion:

lambda <- 1e-6  # Adjust this if you get singular matrix errors (try 1e-5 if needed)
grm_reg <- full_grm + diag(lambda, nrow(full_grm))
grm_inv_reg <- solve(grm_reg)

Bonus: Sparse Matrix Handling (for large datasets) #

For GRMs with 10k+ samples, use the Matrix package to save memory with sparse operations:

library(Matrix)
# Convert to sparse symmetric matrix
grm_sparse <- Matrix(full_grm, sparse = TRUE, symmetric = TRUE)
# Regularize and invert
grm_sparse_reg <- grm_sparse + Diagonal(nrow(grm_sparse), lambda)
grm_inv_sparse <- solve(grm_sparse_reg)

内容的提问来源于stack exchange，提问作者dcp1234