diff --git a/DESCRIPTION b/DESCRIPTION
index c4b14df..6cdb205 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: GFA
 Title: Genetic Factor Analysis (GFA)
 Version: 1.0.0.0456
-Date: 2026-02-18
+Date: 2026-07-08
 Authors@R: 
     person(given = "Jean",
            family = "Morrison",
@@ -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/estimate_R.R b/R/estimate_R.R
index b0470a8..dd0e07c 100644
--- a/R/estimate_R.R
+++ b/R/estimate_R.R
@@ -57,32 +57,56 @@ 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")
   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)
   }
 
-  # 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)
+
+    # 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(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 argument together.")
+    }
+    if(!ncol(comparisons) == 2){
+      stop("Comparisons should have two columns")
+    }
+    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")
+    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],
@@ -94,27 +118,23 @@ R_ldsc <- function(Z_hat,
     rg
   })
 
-  val_int <- map(vals, 1) %>% unlist()
+  res$intercept <- map(vals, 1) %>% unlist()
 
-  res$value <- val_int
-  res_copy <- filter(res, trait1 != trait2) %>%
-    rename(n1c = trait2, n2c = trait1) %>%
-    rename(trait1 = n1c, trait2 = n2c)
+  if(return_matrix){
+    Se <- make_symm_matrix(res, "trait1", "trait2", "intercept", traits_ordered=colnames(Z_hat))
 
-  Se <- bind_rows(res, res_copy)  %>%
-    reshape2::dcast(trait1 ~ trait2, value.var = "value")
-  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,
+    if(make_well_conditioned){
+      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,71 +148,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)
-
-    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
-
-
-    ## 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
+    res$gencov <- val_gencov
+    res$gencor <- val_gencor
+    if(return_matrix){
+      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$value <- val_resid_ve^2
-    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 = "value")
-    ret$Ve <- as.matrix(Ve[,-1])
-    colnames(ret$Ve) <- rownames(ret$Ve) <- NULL
+    res$intercept_var <- val_resid_ve^2
+    if(return_matrix){
+      Ve <- make_symm_matrix(res, "trait1", "trait2", "intercept_var", traits_ordered=colnames(Z_hat))
+    }
     if(return_gencov){
       ## genetic covariance matrix
-      res$value <- val_gencov_ve^2
-      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
-
-
+      res$gencov_var <- val_gencov_ve^2
       ## genetic correlation matrix
-      res$value <- val_gencor_ve^2
-      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
+      res$gencor_var <- val_gencor_ve^2
+
+      if(return_matrix){
+        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))
+      }
     }
   }
+  if(!return_matrix){
+    ret <- res
+  }
   return(ret)
 }
 
 
+# simpler version of make_symm_matrix that reduces copies of data & preserves trait names
+# 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) {
+  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))))
+  }
+
+  # 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_ordered`.")
+  }
+
+  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.
 #'@param B_hat Matrix of effect size estimates
 #'@param S_hat Matrix of standard errors
diff --git a/R/gwas_format.R b/R/gwas_format.R
index 080530d..138243a 100644
--- a/R/gwas_format.R
+++ b/R/gwas_format.R
@@ -118,15 +118,19 @@ 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
+  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")
-    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")
@@ -141,50 +145,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 +160,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(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[, flip_strands_flag :=
+       if (upper) (A1 == "T" | A2 == "T") else (A1 == "t" | A2 == "t")]
+  X[, `:=`(
+    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[, A1_A_flag := flipped_A1 %chin% c("A","a")]
+
+  # swap A1 and A2 if A1 is not A
+  X[, `:=`(
+    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
+  # sd means subset of data (select just beta column)
+  X[, (beta_col) := {
+    b <- .SD[[1L]]
+    if (!is.numeric(b)) b <- as.numeric(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(A1_A_flag, p, 1 - p)
+  }, .SDcols = af_col]
+
+  # delete columns we don't need anymore
+  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
+  return(invisible(X))
+}
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
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{