From 7d5c7f0f132d4082db7141df4fd65dda0f4db6a4 Mon Sep 17 00:00:00 2001 From: Jean Morrison Date: Fri, 20 Feb 2026 11:03:26 -0500 Subject: [PATCH 1/9] add comparisons option to R_ldsc --- DESCRIPTION | 4 +- R/estimate_R.R | 156 ++++++++++++++++++++++++++++--------------------- 2 files changed, 90 insertions(+), 70 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9606cea..6620dc8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: GFA Title: Genetic Factor Analysis (GFA) -Version: 1.0.0.0449 -Date: 2025-11-21 +Version: 1.0.0.0450 +Date: 2026-02-20 Authors@R: person(given = "Jean", family = "Morrison", diff --git a/R/estimate_R.R b/R/estimate_R.R index b0470a8..509a9a5 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -57,7 +57,8 @@ R_ldsc <- function(Z_hat, make_well_conditioned = TRUE, cond_num = 1e5-1, blocks = NULL, - ncores = 1){ + ncores = 1, + comparisons = NULL){ M <- ncol(Z_hat) J <- nrow(Z_hat) stopifnot(class(ldscores) == "numeric") @@ -69,20 +70,20 @@ R_ldsc <- function(Z_hat, N <- matrix(rep(N, each = J), nrow = J) } - # h2 <- lapply(1:M, function(m){ - # ix <- which(!is.na(Z_hat[,m])) - # snp_ldsc(ld_score = ldscores[ix], ld_size = ld_size, chi2 = Z_hat[ix,m]^2, sample_size = N[ix,m], blocks = NULL) - # }) - res <- expand.grid(trait1 = 1:M, trait2 = 1:M) %>% - filter(trait1 <= trait2) + if(is.null(comparisons)){ + res <- expand.grid(trait1 = 1:M, trait2 = 1:M) %>% + filter(trait1 <= trait2) + return_matrix <- TRUE + }else{ + if(make_well_conditioned){ + stop("Cannot use make_well_conditioned and comparisons arguement together.") + } + res <- comparisons + names(res) <- c("trait1", "trait2") + return_matrix <- FALSE + } vals <- map2(res$trait1, res$trait2, function(i, j){ - #cat(i, j, "\n") - # if(i == j){ - # ix <- which(!is.na(Z_hat[,i])) - # h2 <- snp_ldsc(ld_score = ldscores[ix], ld_size = ld_size, chi2 = Z_hat[ix,i]^2, sample_size = N[ix,i], blocks = blocks) - # return(c(h2, 1)) - # } ix <- which(!is.na(Z_hat[,i]) & !is.na(Z_hat[,j])) rg <- ldsc_rg(ld_score = ldscores[ix], ld_size = ld_size, z1 = Z_hat[ix,i], z2 = Z_hat[ix,j], @@ -96,25 +97,32 @@ R_ldsc <- function(Z_hat, val_int <- map(vals, 1) %>% unlist() - res$value <- val_int - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) + res$intercept <- val_int - Se <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "value") - Se <- as.matrix(Se[,-1]) - colnames(Se) <- rownames(Se) <- NULL + if(return_matrix){ + res_copy <- filter(res, trait1 != trait2) %>% + rename(n1c = trait2, n2c = trait1) %>% + rename(trait1 = n1c, trait2 = n2c) - if(make_well_conditioned){ - #Se <- condition(Se, cond_num) - Se <- Matrix::nearPD(Se, corr = FALSE, keepDiag = TRUE, + Se <- bind_rows(res, res_copy) %>% + reshape2::dcast(trait1 ~ trait2, value.var = "intercept") + Se <- as.matrix(Se[,-1]) + colnames(Se) <- rownames(Se) <- NULL + + + if(make_well_conditioned){ + #Se <- condition(Se, cond_num) + Se <- Matrix::nearPD(Se, corr = FALSE, keepDiag = TRUE, posd.tol = 1/cond_num)$mat - Se <- as.matrix(Se) + Se <- as.matrix(Se) + } + ret <- list("Se" = Se) + }else{ + ret <- res } if(!return_gencov & is.null(blocks)){ - return(list("Se" = Se)) + return(ret) } if(!is.null(blocks)){ @@ -128,65 +136,77 @@ R_ldsc <- function(Z_hat, val_gencor <- map(vals, 7) %>% unlist() } - ret <- list("Se" = Se) - if(return_gencov){ ## genetic covariance matrix - res$value <- val_gencov - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) + res$gencov <- val_gencov - Sg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "value") - ret$Sg <- as.matrix(Sg[,-1]) - colnames(ret$Sg) <- rownames(ret$Sg) <- NULL + if(return_matrix){ + res_copy <- filter(res, trait1 != trait2) %>% + rename(n1c = trait2, n2c = trait1) %>% + rename(trait1 = n1c, trait2 = n2c) + Sg <- bind_rows(res, res_copy) %>% + reshape2::dcast(trait1 ~ trait2, value.var = "gencov") + ret$Sg <- as.matrix(Sg[,-1]) + colnames(ret$Sg) <- rownames(ret$Sg) <- NULL + } ## genetic correlation matrix - res$value <- val_gencor - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - Rg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "value") - ret$Rg <- as.matrix(Rg[,-1]) - colnames(ret$Rg) <- rownames(ret$Rg) <- NULL - } - if(!is.null(blocks)){ - res$value <- val_resid_ve^2 - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) + res$gencor <- val_gencor - Ve <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "value") - ret$Ve <- as.matrix(Ve[,-1]) - colnames(ret$Ve) <- rownames(ret$Ve) <- NULL - if(return_gencov){ - ## genetic covariance matrix - res$value <- val_gencov_ve^2 + if(return_matirx){ res_copy <- filter(res, trait1 != trait2) %>% rename(n1c = trait2, n2c = trait1) %>% rename(trait1 = n1c, trait2 = n2c) - Vg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "value") - ret$Vg <- as.matrix(Vg[,-1]) - colnames(ret$Vg) <- rownames(ret$Vg) <- NULL + Rg <- bind_rows(res, res_copy) %>% + reshape2::dcast(trait1 ~ trait2, value.var = "gencor") + ret$Rg <- as.matrix(Rg[,-1]) + colnames(ret$Rg) <- rownames(ret$Rg) <- NULL + } + } + if(!is.null(blocks)){ + res$intercept_var <- val_resid_ve^2 - ## genetic correlation matrix - res$value <- val_gencor_ve^2 + if(return_matrix){ res_copy <- filter(res, trait1 != trait2) %>% rename(n1c = trait2, n2c = trait1) %>% rename(trait1 = n1c, trait2 = n2c) - VRg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "value") - ret$VRg <- as.matrix(VRg[,-1]) - colnames(ret$VRg) <- rownames(ret$VRg) <- NULL + Ve <- bind_rows(res, res_copy) %>% + reshape2::dcast(trait1 ~ trait2, value.var = "intercept_var") + ret$Ve <- as.matrix(Ve[,-1]) + colnames(ret$Ve) <- rownames(ret$Ve) <- NULL + } + if(return_gencov){ + ## genetic covariance matrix + res$gencov_var <- val_gencov_ve^2 + + if(return_matrix){ + res_copy <- filter(res, trait1 != trait2) %>% + rename(n1c = trait2, n2c = trait1) %>% + rename(trait1 = n1c, trait2 = n2c) + + Vg <- bind_rows(res, res_copy) %>% + reshape2::dcast(trait1 ~ trait2, value.var = "gencov_var") + ret$Vg <- as.matrix(Vg[,-1]) + colnames(ret$Vg) <- rownames(ret$Vg) <- NULL + } + + ## genetic correlation matrix + res$gencor_var <- val_gencor_ve^2 + + if(return_matrix){ + res_copy <- filter(res, trait1 != trait2) %>% + rename(n1c = trait2, n2c = trait1) %>% + rename(trait1 = n1c, trait2 = n2c) + + VRg <- bind_rows(res, res_copy) %>% + reshape2::dcast(trait1 ~ trait2, value.var = "gencor_var") + ret$VRg <- as.matrix(VRg[,-1]) + colnames(ret$VRg) <- rownames(ret$VRg) <- NULL + } } } return(ret) From e3a6508f7165d2a04e6e02ba2e72a6805d1d6d4d Mon Sep 17 00:00:00 2001 From: Jean Morrison Date: Fri, 20 Feb 2026 11:34:47 -0500 Subject: [PATCH 2/9] add symmetric matrix function --- DESCRIPTION | 2 +- R/estimate_R.R | 87 +++++++++++++------------------------------------- 2 files changed, 23 insertions(+), 66 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6620dc8..2db2cd0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: GFA Title: Genetic Factor Analysis (GFA) -Version: 1.0.0.0450 +Version: 1.0.0.0451 Date: 2026-02-20 Authors@R: person(given = "Jean", diff --git a/R/estimate_R.R b/R/estimate_R.R index 509a9a5..966871c 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -95,23 +95,12 @@ R_ldsc <- function(Z_hat, rg }) - val_int <- map(vals, 1) %>% unlist() - - res$intercept <- val_int + res$intercept <- map(vals, 1) %>% unlist() if(return_matrix){ - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - Se <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "intercept") - Se <- as.matrix(Se[,-1]) - colnames(Se) <- rownames(Se) <- NULL - + Se <- make_symm_matrix(res, "trait1", "trait2", "intercept") if(make_well_conditioned){ - #Se <- condition(Se, cond_num) Se <- Matrix::nearPD(Se, corr = FALSE, keepDiag = TRUE, posd.tol = 1/cond_num)$mat Se <- as.matrix(Se) @@ -139,80 +128,48 @@ R_ldsc <- function(Z_hat, if(return_gencov){ ## genetic covariance matrix res$gencov <- val_gencov - - if(return_matrix){ - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - Sg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "gencov") - ret$Sg <- as.matrix(Sg[,-1]) - colnames(ret$Sg) <- rownames(ret$Sg) <- NULL - } - - ## genetic correlation matrix res$gencor <- val_gencor - - if(return_matirx){ - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - Rg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "gencor") - ret$Rg <- as.matrix(Rg[,-1]) - colnames(ret$Rg) <- rownames(ret$Rg) <- NULL + if(return_matrix){ + Sg <-make_symm_matrix(res, "trait1", "trait2", "gencov") + Rg <- make_symm_matrix(res, "trait1", "trait2", "gencor") } } if(!is.null(blocks)){ res$intercept_var <- val_resid_ve^2 - if(return_matrix){ - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - Ve <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "intercept_var") - ret$Ve <- as.matrix(Ve[,-1]) - colnames(ret$Ve) <- rownames(ret$Ve) <- NULL + Ve <- make_symm_matrix(res, "trait1", "trait2", "intercept_var") } if(return_gencov){ ## genetic covariance matrix res$gencov_var <- val_gencov_ve^2 - - if(return_matrix){ - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - Vg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "gencov_var") - ret$Vg <- as.matrix(Vg[,-1]) - colnames(ret$Vg) <- rownames(ret$Vg) <- NULL - } - ## genetic correlation matrix res$gencor_var <- val_gencor_ve^2 if(return_matrix){ - res_copy <- filter(res, trait1 != trait2) %>% - rename(n1c = trait2, n2c = trait1) %>% - rename(trait1 = n1c, trait2 = n2c) - - VRg <- bind_rows(res, res_copy) %>% - reshape2::dcast(trait1 ~ trait2, value.var = "gencor_var") - ret$VRg <- as.matrix(VRg[,-1]) - colnames(ret$VRg) <- rownames(ret$VRg) <- NULL + Vg <- make_symm_matrix(res, "trait1", "trait2", "gencov_var") + VRg <- make_symm_matrix(res, "trait1", "trait2", "gencor_var") } } } + if(!return_matrix){ + ret <- res + } return(ret) } +make_symm_matrix <- function(x, row_name, col_name, value_name){ + form <- as.formula(paste0(row_name, "~", col_name)) + x <- dplyr::select(x, !!row_name, !!col_name, !!value_name) + x_copy <- dplyr::filter(x, .data[[row_name]] != .data[[col_name]]) + names(x_copy) <- c(col_name, row_name, value_name) + x <- bind_rows(x, x_copy) %>% reshape2::dcast(form, value.var = value_name) + x <- as.matrix(x[, -1]) + rownames(x) <- colnames(x) <- NULL + return(x) +} + #'@title Calculate matrix of error correlations using p-value threshold method. #'@param B_hat Matrix of effect size estimates #'@param S_hat Matrix of standard errors From fa58a10b002aebd395eb425ecba9ea5dc48f3ba5 Mon Sep 17 00:00:00 2001 From: Jean Morrison Date: Fri, 20 Feb 2026 11:38:56 -0500 Subject: [PATCH 3/9] add some input checking --- DESCRIPTION | 2 +- R/estimate_R.R | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2db2cd0..872c050 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: GFA Title: Genetic Factor Analysis (GFA) -Version: 1.0.0.0451 +Version: 1.0.0.0452 Date: 2026-02-20 Authors@R: person(given = "Jean", diff --git a/R/estimate_R.R b/R/estimate_R.R index 966871c..41fef26 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -78,6 +78,12 @@ R_ldsc <- function(Z_hat, if(make_well_conditioned){ stop("Cannot use make_well_conditioned and comparisons arguement together.") } + if(!ncol(comparisons) == 1){ + stop("comparisons should have two columns") + } + if(! (all(comparisons[,1] %in% 1:M) & all(comparisons[,2] %in% 1:M))){ + stop(paste0("comparisons contains values out of range 1 to ", M)) + } res <- comparisons names(res) <- c("trait1", "trait2") return_matrix <- FALSE From fc839f1cd400c6203a311ba114d743898bc826cf Mon Sep 17 00:00:00 2001 From: Jean Morrison Date: Fri, 27 Feb 2026 10:03:56 -0500 Subject: [PATCH 4/9] Change comparisons column check from 1 to 2 --- R/estimate_R.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/estimate_R.R b/R/estimate_R.R index 41fef26..9089f8b 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -78,7 +78,7 @@ R_ldsc <- function(Z_hat, if(make_well_conditioned){ stop("Cannot use make_well_conditioned and comparisons arguement together.") } - if(!ncol(comparisons) == 1){ + if(!ncol(comparisons) == 2){ stop("comparisons should have two columns") } if(! (all(comparisons[,1] %in% 1:M) & all(comparisons[,2] %in% 1:M))){ From 1e97273611ac105af7ddc25d257a4b2b3dc83b45 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Fri, 24 Apr 2026 12:22:48 -0400 Subject: [PATCH 5/9] changed align_beta in gwas_format.R to a data.table version --- R/gwas_format.R | 111 ++++++++++++++++++++++++++++-------------------- man/R_ldsc.Rd | 3 +- 2 files changed, 66 insertions(+), 48 deletions(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 080530d..951410c 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -118,11 +118,14 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, X <- remove_ambiguous(X, upper = TRUE) cat("Removed ", n-nrow(X), " variants with ambiguous strand.\n") + # make X into a data.table for I can do the new align_beta on it. ideally this would be for everything + data.table::setDT(X) + cat("Flipping strand and effect allele so A1 is always A\n") - X <- align_beta(X, "beta_hat", "allele_freq", TRUE) + align_beta(X) - - X <- X %>% select(chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, allele_freq) + # data table syntax + X <- X[, .(chrom, pos, snp, A1, A2, beta_hat, se, p_value, sample_size, allele_freq)] if(!missing(output_file)){ cat("Writing out ", nrow(X), " variants to file.\n") @@ -141,50 +144,6 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, # return(dat) # } -#Flip signs and strands so that allele 1 is always A -align_beta <- function(X, beta_hat_name, af_name, upper=TRUE){ - flp = c("A" = "T", "G" = "C", "T" = "A", - "C" = "G", "a" = "t", "t" = "a", - "c" = "g", "g" = "c") - if(upper){ - X <- X %>% mutate( flip_strand = A1 == "T" | A2 == "T") - }else{ - X <- X %>% mutate( flip_strand = A1 == "t" | A2 == "t") - } - if(missing(af_name)){ - X <- mutate(X, af = NA) - af_name <- "tempaf" - af_missing <- TRUE - }else{ - af_missing <- FALSE - } - - X <- X %>% mutate(A1flp = case_when(flip_strand ~ flp[A1], - TRUE ~ A1), - A2flp = case_when(flip_strand ~ flp[A2], - TRUE ~ A2), - # afflp = case_when(flip_strand ~ 1-get(af_name), - # TRUE ~ get(af_name)), - tempbh = case_when(A1flp == "A" | A1flp == "a" ~ get(beta_hat_name), - TRUE ~ -1*get(beta_hat_name)), - tempaf = case_when(A1flp == "A" | A1flp == "a" ~ get(af_name), - TRUE ~ 1-get(af_name))) %>% - select(-A1, -A2) %>% - select(-all_of(c(af_name, beta_hat_name))) %>% - mutate(A1 = case_when(A1flp == "A" | A1flp=="a" ~ A1flp, - TRUE ~ A2flp), - A2 = case_when(A1flp == "A" | A1flp=="a" ~ A2flp, - TRUE ~ A1flp)) %>% - select(-A1flp, -A2flp, -flip_strand) - - ix <- which(names(X)== "tempbh") - names(X)[ix] <- beta_hat_name - ix <- which(names(X)== "tempaf") - names(X)[ix] <- af_name - if(af_missing) X <- select(X, -tempaf) - return(X) -} - remove_ambiguous <- function(X, upper=TRUE){ if(upper){ X <- X %>% dplyr::filter(!(A1 == "G" & A2 == "C") & @@ -200,3 +159,61 @@ remove_ambiguous <- function(X, upper=TRUE){ return(X) } +# flip signs and strands so that allele 1 is always A +# now modifies X in-place w/ data table for speed and memory savings +# believe beta_hat and af are always the names assigned in gwas_format. but kept for consistency w/ old func +align_beta <- function(X, upper = TRUE, + beta_col = "beta_hat", + af_col = "af") { + + stopifnot(data.table::is.data.table(X)) + stopifnot(all(c("A1", "A2") %in% names(X))) + stopifnot(beta_col %in% names(X)) + + flp <- c("A"="T","G"="C","T"="A","C"="G", + "a"="t","t"="a","c"="g","g"="c") + + af_present <- af_col %in% names(X) + if (!af_present) { + X[, (af_col) := NA_real_] # create af col temporarily so code can be uniform + } + + # flip strands if we have Ts to get As + X[, ..__alignbeta_flip__ := + if (upper) (A1 == "T" | A2 == "T") else (A1 == "t" | A2 == "t")] + X[, `:=`( + ..__alignbeta_A1flp__ = data.table::fifelse(..__alignbeta_flip__, flp[A1], A1), + ..__alignbeta_A2flp__ = data.table::fifelse(..__alignbeta_flip__, flp[A2], A2) + )] + + # flag that is true if the A1 is A (we want) + X[, ..__alignbeta_condA__ := ..__alignbeta_A1flp__ %chin% c("A","a")] + + # swap A1 and A2 if A1 is not A + X[, `:=`( + A1 = data.table::fifelse(..__alignbeta_condA__, ..__alignbeta_A1flp__, ..__alignbeta_A2flp__), + A2 = data.table::fifelse(..__alignbeta_condA__, ..__alignbeta_A2flp__, ..__alignbeta_A1flp__) + )] + + # flip beta hats if the A1 is not A + # sd means subset of data (select just beta column) + X[, (beta_col) := { + b <- .SD[[1L]] + if (!is.numeric(b)) b <- as.numeric(b) + data.table::fifelse(..__alignbeta_condA__, b, -b) + }, .SDcols = beta_col] + + # flip af if the A1 is not A + X[, (af_col) := { + p <- .SD[[1L]] + if (!is.numeric(p)) p <- as.numeric(p) + data.table::fifelse(..__alignbeta_condA__, p, 1 - p) + }, .SDcols = af_col] + + # delete columns we don't need anymore + X[, c("..__alignbeta_flip__", "..__alignbeta_A1flp__", "..__alignbeta_A2flp__", "..__alignbeta_condA__") := NULL] + if (!af_present) X[, (af_col) := NULL] + + # since these are in-place mods, we can call func w/o assignment to new var. but nobody wants to see the whole table + return(invisible(X)) +} diff --git a/man/R_ldsc.Rd b/man/R_ldsc.Rd index 6e1317f..17b3051 100644 --- a/man/R_ldsc.Rd +++ b/man/R_ldsc.Rd @@ -13,7 +13,8 @@ R_ldsc( make_well_conditioned = TRUE, cond_num = 1e+05 - 1, blocks = NULL, - ncores = 1 + ncores = 1, + comparisons = NULL ) } \arguments{ From a6e28f6a8b3726c166b0821e298e598b7926577b Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Sun, 3 May 2026 11:04:14 -0400 Subject: [PATCH 6/9] added import for data.table --- DESCRIPTION | 3 ++- NAMESPACE | 1 + R/gwas_format.R | 19 ++++++++++--------- R/package.R | 1 + 4 files changed, 14 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 872c050..1689eba 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -21,7 +21,8 @@ Imports: intervals, purrr, lpSolve, - ggplot2 + ggplot2, + data.table Suggests: rmarkdown, knitr diff --git a/NAMESPACE b/NAMESPACE index 157e060..485f9f0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,6 +17,7 @@ export(min_norm) export(plot_factors) export(plot_factors_bars) export(snp_ldsc) +import(data.table) import(dplyr) import(flashier) import(ggplot2) diff --git a/R/gwas_format.R b/R/gwas_format.R index 951410c..21c04bb 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -119,7 +119,7 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, cat("Removed ", n-nrow(X), " variants with ambiguous strand.\n") # make X into a data.table for I can do the new align_beta on it. ideally this would be for everything - data.table::setDT(X) + setDT(X) cat("Flipping strand and effect allele so A1 is always A\n") align_beta(X) @@ -129,7 +129,8 @@ gwas_format <- function(X, snp, beta_hat, se, A1, A2, if(!missing(output_file)){ cat("Writing out ", nrow(X), " variants to file.\n") - write_tsv(X, path = output_file) + # changed from path= to file= + write_tsv(X, file = output_file) return(NULL) } cat("Returning ", nrow(X), " variants.\n") @@ -166,7 +167,7 @@ align_beta <- function(X, upper = TRUE, beta_col = "beta_hat", af_col = "af") { - stopifnot(data.table::is.data.table(X)) + stopifnot(is.data.table(X)) stopifnot(all(c("A1", "A2") %in% names(X))) stopifnot(beta_col %in% names(X)) @@ -182,8 +183,8 @@ align_beta <- function(X, upper = TRUE, X[, ..__alignbeta_flip__ := if (upper) (A1 == "T" | A2 == "T") else (A1 == "t" | A2 == "t")] X[, `:=`( - ..__alignbeta_A1flp__ = data.table::fifelse(..__alignbeta_flip__, flp[A1], A1), - ..__alignbeta_A2flp__ = data.table::fifelse(..__alignbeta_flip__, flp[A2], A2) + ..__alignbeta_A1flp__ = fifelse(..__alignbeta_flip__, flp[A1], A1), + ..__alignbeta_A2flp__ = fifelse(..__alignbeta_flip__, flp[A2], A2) )] # flag that is true if the A1 is A (we want) @@ -191,8 +192,8 @@ align_beta <- function(X, upper = TRUE, # swap A1 and A2 if A1 is not A X[, `:=`( - A1 = data.table::fifelse(..__alignbeta_condA__, ..__alignbeta_A1flp__, ..__alignbeta_A2flp__), - A2 = data.table::fifelse(..__alignbeta_condA__, ..__alignbeta_A2flp__, ..__alignbeta_A1flp__) + A1 = fifelse(..__alignbeta_condA__, ..__alignbeta_A1flp__, ..__alignbeta_A2flp__), + A2 = fifelse(..__alignbeta_condA__, ..__alignbeta_A2flp__, ..__alignbeta_A1flp__) )] # flip beta hats if the A1 is not A @@ -200,14 +201,14 @@ align_beta <- function(X, upper = TRUE, X[, (beta_col) := { b <- .SD[[1L]] if (!is.numeric(b)) b <- as.numeric(b) - data.table::fifelse(..__alignbeta_condA__, b, -b) + fifelse(..__alignbeta_condA__, b, -b) }, .SDcols = beta_col] # flip af if the A1 is not A X[, (af_col) := { p <- .SD[[1L]] if (!is.numeric(p)) p <- as.numeric(p) - data.table::fifelse(..__alignbeta_condA__, p, 1 - p) + fifelse(..__alignbeta_condA__, p, 1 - p) }, .SDcols = af_col] # delete columns we don't need anymore diff --git a/R/package.R b/R/package.R index f016673..d121965 100644 --- a/R/package.R +++ b/R/package.R @@ -7,5 +7,6 @@ #' @import flashier #' @import dplyr readr reshape2 lpSolve #' @import purrr ggplot2 +#' @import data.table #' @name GFA NULL From 7743a840370a056b5845a0433d8f56c488df517a Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Sun, 3 May 2026 14:30:13 -0400 Subject: [PATCH 7/9] improved align_beta readability --- R/gwas_format.R | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/R/gwas_format.R b/R/gwas_format.R index 21c04bb..138243a 100644 --- a/R/gwas_format.R +++ b/R/gwas_format.R @@ -180,20 +180,20 @@ align_beta <- function(X, upper = TRUE, } # flip strands if we have Ts to get As - X[, ..__alignbeta_flip__ := + X[, flip_strands_flag := if (upper) (A1 == "T" | A2 == "T") else (A1 == "t" | A2 == "t")] X[, `:=`( - ..__alignbeta_A1flp__ = fifelse(..__alignbeta_flip__, flp[A1], A1), - ..__alignbeta_A2flp__ = fifelse(..__alignbeta_flip__, flp[A2], A2) + flipped_A1 = fifelse(flip_strands_flag, flp[A1], A1), + flipped_A2 = fifelse(flip_strands_flag, flp[A2], A2) )] # flag that is true if the A1 is A (we want) - X[, ..__alignbeta_condA__ := ..__alignbeta_A1flp__ %chin% c("A","a")] + X[, A1_A_flag := flipped_A1 %chin% c("A","a")] # swap A1 and A2 if A1 is not A X[, `:=`( - A1 = fifelse(..__alignbeta_condA__, ..__alignbeta_A1flp__, ..__alignbeta_A2flp__), - A2 = fifelse(..__alignbeta_condA__, ..__alignbeta_A2flp__, ..__alignbeta_A1flp__) + A1 = fifelse(A1_A_flag, flipped_A1, flipped_A2), + A2 = fifelse(A1_A_flag, flipped_A2, flipped_A1) )] # flip beta hats if the A1 is not A @@ -201,18 +201,18 @@ align_beta <- function(X, upper = TRUE, X[, (beta_col) := { b <- .SD[[1L]] if (!is.numeric(b)) b <- as.numeric(b) - fifelse(..__alignbeta_condA__, b, -b) + fifelse(A1_A_flag, b, -b) }, .SDcols = beta_col] # flip af if the A1 is not A X[, (af_col) := { p <- .SD[[1L]] if (!is.numeric(p)) p <- as.numeric(p) - fifelse(..__alignbeta_condA__, p, 1 - p) + fifelse(A1_A_flag, p, 1 - p) }, .SDcols = af_col] # delete columns we don't need anymore - X[, c("..__alignbeta_flip__", "..__alignbeta_A1flp__", "..__alignbeta_A2flp__", "..__alignbeta_condA__") := NULL] + X[, c("flip_strands_flag", "flipped_A1", "flipped_A2", "A1_A_flag") := NULL] if (!af_present) X[, (af_col) := NULL] # since these are in-place mods, we can call func w/o assignment to new var. but nobody wants to see the whole table From 5946abd5ce70e0eac7b27e1615cc099e682bf016 Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Tue, 12 May 2026 11:10:35 -0400 Subject: [PATCH 8/9] make R_ldsc handle trait names in comparisons --- R/estimate_R.R | 96 ++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 73 insertions(+), 23 deletions(-) diff --git a/R/estimate_R.R b/R/estimate_R.R index 9089f8b..a9976ab 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -64,25 +64,42 @@ R_ldsc <- function(Z_hat, stopifnot(class(ldscores) == "numeric") stopifnot(length(ldscores) == J) if("matrix" %in% class(N)){ - stopifnot(nrow(N) == J & ncol(N) == M) + stopifnot(nrow(N) == J & identical(colnames(N),colnames(Z_hat))) }else if(class(N) == "numeric"){ - stopifnot(length(N) == M) + stopifnot(identical(names(N),colnames(Z_hat))) N <- matrix(rep(N, each = J), nrow = J) } if(is.null(comparisons)){ - res <- expand.grid(trait1 = 1:M, trait2 = 1:M) %>% - filter(trait1 <= trait2) + #res <- expand.grid(trait1 = 1:M, trait2 = 1:M) %>% + # filter(trait1 <= trait2) + + # data.table better for ram, but not necessary to do all that here + # generate M (M+1) / 2 unique pairs instead of expand M^2 and filter to M (M+1)/2 + + # trait1_idx: 1, 1, 2, 1, 2, 3, 1, 2, 3, 4 + # trait2_idx: 1, 2, 2, 3, 3, 3, 4, 4, 4, 4 + + res <- data.frame( + trait1_idx = sequence(seq_len(M)), + trait2_idx = rep.int(seq_len(M), times = seq_len(M)) + ) + # use trait names instead of indices + trait_names <- colnames(Z_hat) + res$trait1 <- trait_names[res$trait1_idx] + res$trait2 <- trait_names[res$trait2_idx] + return_matrix <- TRUE + }else{ if(make_well_conditioned){ - stop("Cannot use make_well_conditioned and comparisons arguement together.") + stop("Cannot use make_well_conditioned and comparisons argument together.") } if(!ncol(comparisons) == 2){ - stop("comparisons should have two columns") + stop("Comparisons should have two columns") } - if(! (all(comparisons[,1] %in% 1:M) & all(comparisons[,2] %in% 1:M))){ - stop(paste0("comparisons contains values out of range 1 to ", M)) + if(! (all(comparisons[,1] %in% colnames(Z_hat)) & all(comparisons[,2] %in% colnames(Z_hat)))){ + stop(paste0("Comparisons contains values out of provided trait names in Z_hat")) } res <- comparisons names(res) <- c("trait1", "trait2") @@ -104,7 +121,7 @@ R_ldsc <- function(Z_hat, res$intercept <- map(vals, 1) %>% unlist() if(return_matrix){ - Se <- make_symm_matrix(res, "trait1", "trait2", "intercept") + Se <- make_symm_matrix(res, "trait1", "trait2", "intercept", traits_ordered=colnames(Z_hat)) if(make_well_conditioned){ Se <- Matrix::nearPD(Se, corr = FALSE, keepDiag = TRUE, @@ -136,15 +153,15 @@ R_ldsc <- function(Z_hat, res$gencov <- val_gencov res$gencor <- val_gencor if(return_matrix){ - Sg <-make_symm_matrix(res, "trait1", "trait2", "gencov") - Rg <- make_symm_matrix(res, "trait1", "trait2", "gencor") + Sg <-make_symm_matrix(res, "trait1", "trait2", "gencov", traits_ordered=colnames(Z_hat)) + Rg <- make_symm_matrix(res, "trait1", "trait2", "gencor", traits_ordered=colnames(Z_hat)) } } if(!is.null(blocks)){ res$intercept_var <- val_resid_ve^2 if(return_matrix){ - Ve <- make_symm_matrix(res, "trait1", "trait2", "intercept_var") + Ve <- make_symm_matrix(res, "trait1", "trait2", "intercept_var", traits_ordered=colnames(Z_hat)) } if(return_gencov){ ## genetic covariance matrix @@ -153,8 +170,8 @@ R_ldsc <- function(Z_hat, res$gencor_var <- val_gencor_ve^2 if(return_matrix){ - Vg <- make_symm_matrix(res, "trait1", "trait2", "gencov_var") - VRg <- make_symm_matrix(res, "trait1", "trait2", "gencor_var") + Vg <- make_symm_matrix(res, "trait1", "trait2", "gencov_var", traits_ordered=colnames(Z_hat)) + VRg <- make_symm_matrix(res, "trait1", "trait2", "gencor_var", traits_ordered=colnames(Z_hat)) } } } @@ -165,15 +182,48 @@ R_ldsc <- function(Z_hat, } -make_symm_matrix <- function(x, row_name, col_name, value_name){ - form <- as.formula(paste0(row_name, "~", col_name)) - x <- dplyr::select(x, !!row_name, !!col_name, !!value_name) - x_copy <- dplyr::filter(x, .data[[row_name]] != .data[[col_name]]) - names(x_copy) <- c(col_name, row_name, value_name) - x <- bind_rows(x, x_copy) %>% reshape2::dcast(form, value.var = value_name) - x <- as.matrix(x[, -1]) - rownames(x) <- colnames(x) <- NULL - return(x) +#make_symm_matrix <- function(x, row_name, col_name, value_name){ +# form <- as.formula(paste0(row_name, "~", col_name)) +# x <- dplyr::select(x, !!row_name, !!col_name, !!value_name) +# x_copy <- dplyr::filter(x, .data[[row_name]] != .data[[col_name]]) +# names(x_copy) <- c(col_name, row_name, value_name) +# x <- bind_rows(x, x_copy) %>% reshape2::dcast(form, value.var = value_name) +# x <- as.matrix(x[, -1]) +# rownames(x) <- colnames(x) <- NULL +# return(x) +#} + +# simpler version of make_symm_matrix that reduces copies of data & preserves trait names +# if you leave traits null, we will infer it and give back alphabetically. or can give traits to get the exact order you wanted +make_symm_matrix <- function(x, row_name, col_name, value_name, + traits_ordered = NULL, + keep_trait_names = TRUE) { + r <- as.character(x[[row_name]]) + c <- as.character(x[[col_name]]) + v <- x[[value_name]] + + if (is.null(traits_ordered)) { + traits_ordered <- sort(unique(c(as.character(r), as.character(c)))) + } + + ri <- match(r, traits_ordered) + ci <- match(c, traits_ordered) + + if (anyNA(ri) || anyNA(ci)) { + stop("Some row/column labels in `x` were not found in `traits`.") + } + + out <- matrix( + NA_real_, + nrow = length(traits_ordered), + ncol = length(traits_ordered), + dimnames = if (keep_trait_names) list(traits_ordered, traits_ordered) else NULL + ) + + out[cbind(ri, ci)] <- v + out[cbind(ci, ri)] <- v + + return(out) } #'@title Calculate matrix of error correlations using p-value threshold method. From 50dbfb5e222b6cfaf87a8810d0f88507392b212d Mon Sep 17 00:00:00 2001 From: Charlotte Edwards Date: Mon, 18 May 2026 10:46:33 -0400 Subject: [PATCH 9/9] improved readability of changes to estimate_R.R --- R/estimate_R.R | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/R/estimate_R.R b/R/estimate_R.R index a9976ab..dd0e07c 100644 --- a/R/estimate_R.R +++ b/R/estimate_R.R @@ -82,7 +82,7 @@ R_ldsc <- function(Z_hat, res <- data.frame( trait1_idx = sequence(seq_len(M)), - trait2_idx = rep.int(seq_len(M), times = seq_len(M)) + trait2_idx = rep(seq_len(M), times = seq_len(M)) ) # use trait names instead of indices trait_names <- colnames(Z_hat) @@ -182,19 +182,8 @@ R_ldsc <- function(Z_hat, } -#make_symm_matrix <- function(x, row_name, col_name, value_name){ -# form <- as.formula(paste0(row_name, "~", col_name)) -# x <- dplyr::select(x, !!row_name, !!col_name, !!value_name) -# x_copy <- dplyr::filter(x, .data[[row_name]] != .data[[col_name]]) -# names(x_copy) <- c(col_name, row_name, value_name) -# x <- bind_rows(x, x_copy) %>% reshape2::dcast(form, value.var = value_name) -# x <- as.matrix(x[, -1]) -# rownames(x) <- colnames(x) <- NULL -# return(x) -#} - # simpler version of make_symm_matrix that reduces copies of data & preserves trait names -# if you leave traits null, we will infer it and give back alphabetically. or can give traits to get the exact order you wanted +# if you leave traits_ordered null, we will give order back alphabetically make_symm_matrix <- function(x, row_name, col_name, value_name, traits_ordered = NULL, keep_trait_names = TRUE) { @@ -206,11 +195,15 @@ make_symm_matrix <- function(x, row_name, col_name, value_name, traits_ordered <- sort(unique(c(as.character(r), as.character(c)))) } + # strip NA entries from traits_ordered + # x will have no NAs, but the colnames of Z_hat might because of the construction of Z_work in 3_R_ldsc_strip.R + traits_ordered <- traits_ordered[!is.na(traits_ordered)] + ri <- match(r, traits_ordered) ci <- match(c, traits_ordered) if (anyNA(ri) || anyNA(ci)) { - stop("Some row/column labels in `x` were not found in `traits`.") + stop("Some row/column labels in `x` were not found in `traits_ordered`.") } out <- matrix(