Possible inversion of LOCO logic in build_k_matrix function

[gzygaozhiyuan]

Hi Dr. Endelman,
I’ve been benchmarking the kinship-calculating section of the script on GitHub and noticed a small but critical inversion in the way the LOCO (Leave-One-Chromosome-Out) matrices are built.
Current code (set.K):

f1 <- function(marks, data) {
      M <- scale(data@geno[, marks], center = T, scale = F)
      K <- tcrossprod(M)
      K <- K/mean(diag(K))

markers <- split(data@map$Marker, data@map$Chrom)   # per-chr SNPs
K <- lapply(markers, f1, data = data)               # f1 uses those SNPs

This produces “chromosome-specific” kinship matrices, i.e. each matrix is built from the SNPs on the same chromosome that will later be tested. In the standard LOCO implementation, we do the opposite: when we scan chromosome i we estimate kinship from all SNPs except those on chromosome i, so that the random-effect term does not absorb the very signal we are trying to detect.

allSNPidx <- seq_len(ncol(data@geno))
markers   <- split(data@map$Marker, data@map$Chrom)
K <- lapply(names(markers), function(chr) {
  out   <- which(data@map$Chrom == chr)  
  M     <- scale(data@geno[, -out], center = T, scale = F)
  tcrossprod(M) / mean(diag(tcrossprod(M)))
})

That would make the resulting list K[[i]] correspond to the kinship matrix used when scanning chromosome i, consistent with the usual LOCO definition.

Thanks very much for maintaining the code—it's been incredibly useful.

Best regards,
Zhiyuan
China

[jendelman]

Hello Zhiyuan, I am glad when people check my code for mistakes!
As you pointed out, that code snippet from set.K produces a list of K matrices using the markers on each chromosome. But the LOCO part is implemented in the GWASpoly function:

LOCO <- length(data@K) > 1
if (LOCO) {
  K <- vector("list",n.chr)
  for (i in 1:n.chr) {
    K[[i]] <- makeLOCO(data@K,exclude=i)
  }

And the makeLOCO function is defined separately:

makeLOCO <- function(K,exclude) {
  #K is list
  #exclude is vector of indices
  n.chr <- length(K)
  tmp <- K[[1]]*0
  keep <- setdiff(1:n.chr,exclude)
  for (i in keep) {
    tmp <- tmp + K[[i]]
  }
  return(tmp/length(keep))
}

Your solution is much simpler and would suffice for the genome scan. But a kinship matrix is also needed in fit.QTL, and when there are multiple QTL on different chromosomes, the approach I implemented allows for rapid construction of kinship matrices for each situation.