Non-numeric argument to mathematical function in XtX_pro

[guarelin]

I am trying to called mesusie_core and getting the following error:

Error in sqrt(XtX.diag[[x]]) :
  non-numeric argument to mathematical function
Calls: XtX_pro -> lapply -> FUN -> diag
Execution halted

I traced through the source code and first had to troubleshoot the issue that the underlying code needs a column called "Se" not "SE." Then, I was able to recreate the error from the package.

Here is the code I'm running:

library(dplyr)
library(MESuSiE)

args <- commandArgs(TRUE)
paths <- strsplit(args, " ")
sumstats_paths <- paths[grep("sumstats", paths)]
LD_mtx_paths <- paths[grep("LD_mtx", paths)]
ancestries <- setdiff(paths, union(sumstats_paths, LD_mtx_paths))

sumstats_frames <- lapply(sumstats_paths, function(path) {
  df <- read.csv(path)
  df$SNP <- gsub(":", ".", df$SNP)
  df$Z <- df$BETA / df$SE
  df$Se <- df$SE
  return(df)
})
names(sumstats_frames) <- ancestries

LD_mtx_frames <- lapply(LD_mtx_paths, function(path) {
  df <- read.csv(path, row.names=1)
  rownames(df) <- gsub(":", ".", rownames(df))
  return(df)
})

names(LD_mtx_frames) <- ancestries

gwas_snps <- lapply(sumstats_frames, function(df) {
  df$SNP
})
gwas_snps <- Reduce(intersect, gwas_snps)

ld_snps <- lapply(LD_mtx_frames, function(df) {
  rownames(df)
})
ld_snps <- Reduce(intersect, ld_snps)

snp_list <-intersect(gwas_snps, ld_snps)
print(length(snp_list))

sumstats_frames <- lapply(sumstats_frames, function(df) {
  df[df$SNP %in% snp_list,]
})

LD_mtx_frames <- lapply(LD_mtx_frames, function(df) {
  df[snp_list, snp_list]
})

###Process XtX
XtX_pro <- function(R, XtX.diag, Name_list){
  return(lapply(Name_list,function(x){
    print(XtX.diag[[x]])
    return(diag(sqrt(XtX.diag[[x]])) %*% R[[x]] %*% diag(sqrt(XtX.diag[[x]])))
  }))
}

var_y <- 1

XtX_diag <- function(Summary_Stat, Name_list){
    return( lapply(Name_list, function(x){
      R2 <- (Summary_Stat[[x]]$Z^2)/(Summary_Stat[[x]]$Z^2+Summary_Stat[[x]]$N-2)
      sigma2 <- var_y*(1-R2)*(Summary_Stat[[x]]$N-1)/(Summary_Stat[[x]]$N-2)
      return(sigma2/(Summary_Stat[[x]]$Se)^2)
    }))
  }

Name_list <- names(sumstats_frames)
R <- LD_mtx_frames
Summary_Stat <- sumstats_frames

XtX.diag <- XtX_diag(Summary_Stat, Name_list)
print(XtX.diag)
print(names(XtX.diag))

XtX_list <- XtX_pro(R, XtX.diag, Name_list)
print(XtX_list)

The issue occurs because in the XtX_pro function, XtX.diag has no names. When the XtX_diag function is called, the names are lost. I'm not sure if this is possibly an R versioning issue, but please let me know if there's something I'm doing wrong!