Move some functions back to pecotmr

DESCRIPTION

@@ -18,6 +18,7 @@ Imports:
     pecotmr,
     readr,
     stats,
+    susieR,
     tidyr,
     vroom
 Suggests:

NAMESPACE

@@ -9,18 +9,24 @@ export(perform_qr_analysis)
 export(qr_screen)
 export(quantile_twas_pipeline)
 export(quantile_twas_weight_pipeline)
-importFrom(bigsnpr,snp_clumping)
-importFrom(bigstatsr,FBM.code256)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(parallel,mclapply)
+importFrom(pecotmr,clean_context_names)
+importFrom(pecotmr,enforce_design_full_rank)
+importFrom(pecotmr,filter_variants_by_ld_reference)
+importFrom(pecotmr,find_data)
+importFrom(pecotmr,harmonize_twas)
+importFrom(pecotmr,ld_clump_by_score)
+importFrom(pecotmr,ld_prune_by_correlation)
+importFrom(pecotmr,parse_variant_id)
+importFrom(pecotmr,twas_analysis)
 importFrom(stats,as.dist)
-importFrom(stats,cor)
 importFrom(stats,cov)
 importFrom(stats,cutree)
 importFrom(stats,hclust)
 importFrom(stats,rnorm)
-importFrom(stats,var)
+importFrom(susieR,univariate_regression)
 importFrom(tidyr,pivot_wider)
 importFrom(tidyr,separate)

R/qQTLR-package.R

@@ -0,0 +1,11 @@
+#' qQTLR: Quantile QTL Analysis
+#'
+#' Tools for quantile QTL analysis, including QUAIL vQTL mapping,
+#' quantile TWAS weight training, and quantile TWAS pipelines.
+#'
+#' @importFrom pecotmr parse_variant_id harmonize_twas twas_analysis find_data
+#'   clean_context_names filter_variants_by_ld_reference
+#'   ld_prune_by_correlation ld_clump_by_score enforce_design_full_rank
+#' @importFrom susieR univariate_regression
+#' @keywords internal
+"_PACKAGE"

R/quail_vqtl.R

@@ -1,95 +1,44 @@
-#' Perform Univariate Regression for Genotype Data
+#' Univariate regression with z/p/q augmentation
 #'
-#' @param X Numeric matrix of genotypes (n x p), where n is the number of samples and p is the number of SNPs.
+#' Thin wrapper around \code{susieR::univariate_regression} that strips rows
+#' with NA in \code{y}, then augments the returned list with z-scores,
+#' p-values, q-values, and names derived from \code{colnames(X)}.
+#'
+#' @param X Numeric matrix of genotypes (n x p).
 #' @param y Numeric vector of phenotypes (length n).
-#' @param Z Optional numeric matrix of covariates (n x k), where k is the number of covariates.
-#' @param center Logical, whether to center the data (default: TRUE).
-#' @param scale Logical, whether to scale the data (default: FALSE).
+#' @param Z Optional numeric matrix of covariates (n x k).
+#' @param center Logical, whether to center (default: TRUE).
+#' @param scale Logical, whether to scale (default: FALSE).
 #' @param return_residuals Logical, whether to return residuals (default: FALSE).
-#' @return A list containing regression results: \code{betahat}, \code{sebetahat}, \code{z_scores}, \code{p_values}, and \code{q_values}.
-#' @examples
-#' X <- matrix(rnorm(1000), 100, 10)
-#' y <- rnorm(100)
-#' results <- univariate_regression(X, y)
-#' @noRd
-calc_stderr <- function(X, residuals) {
-  sqrt(diag(sum(residuals^2) / (nrow(X) - 2) * chol2inv(chol(crossprod(X)))))
-}
-
+#' @return A list with \code{betahat}, \code{sebetahat}, \code{z_scores},
+#'   \code{p_values}, \code{q_values}, and optionally \code{residuals}.
 #' @noRd
 univariate_regression <- function(X, y, Z = NULL, center = TRUE,
                                   scale = FALSE, return_residuals = FALSE) {
   y_na <- which(is.na(y))
   if (length(y_na)) {
     X <- X[-y_na, ]
     y <- y[-y_na]
+    if (!is.null(Z)) Z <- Z[-y_na, , drop = FALSE]
   }
 
-  if (center) {
-    y <- y - mean(y)
-    X <- scale(X, center = TRUE, scale = scale)
-  } else {
-    X <- scale(X, center = FALSE, scale = scale)
-  }
-
-  X[is.nan(X)] <- 0
-
-  if (!is.null(Z)) {
-    if (center) {
-      Z <- scale(Z, center = TRUE, scale = scale)
-    }
-    residual_x <- X # Store X before modifying y
-    y <- .lm.fit(Z, y)$residuals
-  } else {
-    residual_x <- X
-  }
-
-  output <- try(do.call(
-    rbind,
-    lapply(1:ncol(X), function(i) {
-      g <- .lm.fit(cbind(1, X[, i]), y)
-      return(c(coef(g)[2], calc_stderr(cbind(1, X[, i]), g$residuals)[2]))
-    })
-  ), silent = TRUE)
-
-  # Exception handling
-  if (inherits(output, "try-error")) {
-    output <- matrix(0, ncol(X), 2)
-    for (i in 1:ncol(X)) {
-      fit <- summary(lm(y ~ X[, i]))$coef
-      if (nrow(fit) == 2) {
-        output[i, ] <- as.vector(summary(lm(y ~ X[, i]))$coef[2, 1:2])
-      } else {
-        output[i, ] <- c(0, 0)
-      }
-    }
-  }
-
-  # Calculate z-scores (t-statistics)
-  z_scores <- output[, 1] / output[, 2]
-
-  # Calculate p-values from z-scores
-  p_values <- 2 * pnorm(abs(z_scores), lower.tail = FALSE)
-
-  # Calculate q-values using the provided function
-  q_values <- pecotmr:::compute_qvalues(p_values)
-
-  # Use residual_x's column names as the column names for the results
-  rownames(output) <- colnames(residual_x)
-
-  result_list <- list(
-    betahat = setNames(output[, 1], rownames(output)),
-    sebetahat = setNames(output[, 2], rownames(output)),
-    z_scores = setNames(z_scores, rownames(output)),
-    p_values = setNames(p_values, rownames(output)),
-    q_values = setNames(q_values, rownames(output))
+  nm <- colnames(X)
+  result <- susieR::univariate_regression(
+    X, y,
+    Z = Z, center = center, scale = scale,
+    return_residuals = return_residuals
   )
 
-  if (return_residuals && !is.null(Z)) {
-    result_list$residuals <- y
-  }
-
-  return(result_list)
+  z_scores <- result$betahat / result$sebetahat
+  p_values <- 2 * pnorm(abs(z_scores), lower.tail = FALSE)
+  q_values <- compute_qvalues(p_values)
+
+  result$betahat   <- setNames(result$betahat,   nm)
+  result$sebetahat <- setNames(result$sebetahat, nm)
+  result$z_scores  <- setNames(z_scores, nm)
+  result$p_values  <- setNames(p_values, nm)
+  result$q_values  <- setNames(q_values, nm)
+  result
 }
 
 
@@ -117,7 +66,7 @@ run_linear_regression <- function(genotype, phenotype, covariates = NULL, phenot
     scale = FALSE
   )
 
-  snp_info <- pecotmr::parse_variant_id(colnames(genotype))
+  snp_info <- parse_variant_id(colnames(genotype))
 
   data.frame(
     phenotype_id = if (!is.null(phenotype_id)) phenotype_id else NA,

R/quantile_twas.R

@@ -377,7 +377,7 @@ load_quantile_twas_weights <- function(weight_db_files, tau_values = seq(0.01, 0
           if (!is.null(context_data$twas_weight) && !is.null(context_data$pseudo_R2)) {
             list(
               gene = gene,
-              context = pecotmr::clean_context_names(context, gene),
+              context = clean_context_names(context, gene),
               data = context_data
             )
           } else {

R/quantile_twas_pipeline.R

@@ -165,7 +165,7 @@ quantile_twas_pipeline <- function(twas_weights_data,
   }
 
   # Step 1: Harmonize TWAS data
-  twas_data_qced_result <- pecotmr::harmonize_twas(twas_weights_data, ld_meta_file_path, gwas_meta_file,
+  twas_data_qced_result <- harmonize_twas(twas_weights_data, ld_meta_file_path, gwas_meta_file,
                                                     ld_reference_sample_size = ld_reference_sample_size,
                                                     column_file_path = column_file_path, comment_string = comment_string)
   twas_results_db <- lapply(names(twas_weights_data), function(weight_db) {
@@ -202,7 +202,7 @@ quantile_twas_pipeline <- function(twas_weights_data,
         if (length(twas_variants) == 0) {
           return(list(twas_rs_df = data.frame()))
         }
-        twas_rs <- pecotmr::twas_analysis(
+        twas_rs <- twas_analysis(
           twas_data_qced[[weight_db]][["weights_qced"]][[context]][[study]][["weights"]],
           twas_data_qced[[weight_db]][["gwas_qced"]][[study]],
           twas_data_qced[[weight_db]][["LD"]], twas_variants
@@ -212,8 +212,8 @@ quantile_twas_pipeline <- function(twas_weights_data,
         }
         twas_rs_df <- data.frame(
           gwas_study = study, method = sub("_[^_]+$", "", names(twas_rs)),
-          twas_z = pecotmr::find_data(twas_rs, c(2, "z")),
-          twas_pval = pecotmr::find_data(twas_rs, c(2, "pval")),
+          twas_z = find_data(twas_rs, c(2, "z")),
+          twas_pval = find_data(twas_rs, c(2, "pval")),
           context = context, molecular_id = weight_db
         )
         return(list(twas_rs_df = twas_rs_df))