Merge pull request #19 from alexmccreight/main

general improvements, documentation, bug fixes

Author: alexmccreight

DESCRIPTION

@@ -16,12 +16,17 @@ Authors@R: c(
 License: MIT + file LICENSE
 Imports:
     matrixStats,
-    Rcpp (>= 1.0.3)
+    Rcpp (>= 1.0.3),
+    susieR,
+    tibble
 Suggests:
     testthat,
     knitr,
-    rmarkdown
+    rmarkdown,
+    BEDMatrix,
+    data.table,
+    mvtnorm
 LinkingTo: Rcpp
 NeedsCompilation: yes
 VignetteBuilder: knitr
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3

NAMESPACE

@@ -1,11 +1,13 @@
 # Generated by roxygen2: do not edit by hand
 
+export(calculate_sumstat)
 export(compute_s2g)
 export(gene_data)
 export(gene_data_null1)
 export(gene_data_null2)
 export(gene_data_null3)
 export(generate_eqtl_data)
+export(get_correlation)
 export(get_lower_chol)
 export(get_random_A)
 export(parse_num_causal_snps)
@@ -18,14 +20,17 @@ export(sim_multi_traits)
 export(sim_sumstats)
 export(simulate_causal_config)
 export(simulate_cis_expression)
-export(simulate_infinitesimal_effects)
-export(simulate_oligogenic_effects)
 export(simulate_polygenic_trait)
-export(simulate_sparse_effects)
 export(simulate_trans_expression)
 export(simulate_trans_mixture_celltype)
+importFrom(BEDMatrix,BEDMatrix)
 importFrom(MASS,mvrnorm)
 importFrom(Matrix,chol)
+importFrom(Matrix,chol2inv)
+importFrom(data.table,fread)
+importFrom(mvtnorm,rmvnorm)
 importFrom(purrr,map2_dfc)
 importFrom(stats,pnorm)
 importFrom(stats,rnorm)
+importFrom(susieR,univariate_regression)
+importFrom(tibble,tibble)

R/ld_utils.R

@@ -1,30 +1,16 @@
-
-
-# - Calculate LD matrix
-get_ld <- function(X, intercepte = FALSE){
-    X = t(X)
-    # Center each variable
-    if (!intercepte){
-        X = X - rowMeans(X)
-    }
-    # Standardize each variable
-    X = X / sqrt(rowSums(X^2))
-    # Calculate correlations
-    cr = tcrossprod(X)
-    return(cr)
-}
-
-# - Get positive definite LD matrix
-# Function to check if a matrix is positive definite
+# Check if a matrix is positive definite
+#' @keywords internal
 is.positive.definite <- function(LD) {
     ev <- eigen(LD, symmetric = TRUE)$values
     return(all(ev > 0))
 }
 
+# Get positive definite LD matrix
+#' @keywords internal
 get_ld_pd <- function(LD, lambda = 1e-3){
 
     if (!isSymmetric(LD))
-        stop("Provided LD matrix should be sysmatric!")
+        stop("Provided LD matrix should be symmetric!")
 
     while (!is.positive.definite(LD)) {
         LD <- LD + lambda * diag(nrow(LD))
@@ -33,29 +19,32 @@ get_ld_pd <- function(LD, lambda = 1e-3){
     return(LD)
 }
 
-
-
-# - Calculate pseduo inverse of LD matrix
+# Calculate pseudo inverse of LD matrix
+#' @importFrom Matrix chol2inv
+#' @keywords internal
 get_ld_inverse <- function(LD, method = "svd", 
                            tol.exact = sqrt(.Machine$double.eps),
                            tol.pct = 1e-3){
 
     if (!isSymmetric(LD))
-        stop("Provided LD matrix should be sysmatric!")
-    
-    if (method == "svd")
+        stop("Provided LD matrix should be symmetric!")
+
+    if (method == "svd") {
         LD_inv <- get_ld_inv_svd(LD, tol = tol.exact)
-    
-    if (method == "eigen")
+    } else if (method == "eigen") {
         LD_inv <- get_ld_inv_eigen(LD, tol = tol.pct)
-    
-    if (method == "fastchol")
-        LD_inv <- Matrix::chol2inv(LD)
-    
+    } else if (method == "fastchol") {
+        LD_inv <- chol2inv(LD)
+    } else {
+        stop("Invalid method. Choose 'svd', 'eigen', or 'fastchol'.")
+    }
+
     return(LD_inv)
-    
+
 }
 
+# Calculate LD inverse using eigendecomposition
+#' @keywords internal
 get_ld_inv_eigen <- function(LD, tol = 1e-3){
 
     eigen_LD <- eigen(LD)
@@ -67,16 +56,17 @@ get_ld_inv_eigen <- function(LD, tol = 1e-3){
     return(LD_inv)
 }
 
+# Calculate LD inverse using SVD
+#' @keywords internal
 get_ld_inv_svd <- function(LD, tol = sqrt(.Machine$double.eps)){
 
     svd_LD <- svd(LD)
     Positive <- svd_LD$d > max(tol * svd_LD$d[1L], 0)
     if (all(Positive))
-        LD_inv = svd_LD$v %*% (1/svd_LD$d * t(svd_LD$u))
-    else 
-        LD_inv = svd_LD$v[,Positive,drop = FALSE] %*%
+        LD_inv <- svd_LD$v %*% (1/svd_LD$d * t(svd_LD$u))
+    else
+        LD_inv <- svd_LD$v[,Positive,drop = FALSE] %*%
         ((1/svd_LD$d[Positive]) * t(svd_LD$u[,Positive,drop = FALSE]))
 
     return(LD_inv)
 }
-

R/simulate_eQTL.R

@@ -14,7 +14,7 @@
 #'   \item{beta}{A vector of effect sizes with nonzero entries for the sparse SNPs.}
 #'   \item{sentinel_index}{Index of the sentinel SNP.}
 #'   \item{other_sparse_indices}{Indices of the additional sparse SNPs.}
-#' @export
+#' @keywords internal
 simulate_sparse_effects <- function(G, h2_sparse, prop_h2_sentinel, n_other_sparse,
                                     min_sparse_effect = 0.10, sentinel_dominance = 1.5) {
   n_features <- ncol(G)
@@ -99,7 +99,7 @@ simulate_sparse_effects <- function(G, h2_sparse, prop_h2_sentinel, n_other_spar
 #'   \item{beta}{A vector of effect sizes for oligogenic effects (zeros elsewhere).}
 #'   \item{oligogenic_indices}{Indices of the oligogenic SNPs.}
 #'   \item{mixture_assignments}{A vector (indexed by SNP) of mixture component assignments.}
-#' @export
+#' @keywords internal
 simulate_oligogenic_effects <- function(G, h2_oligogenic, n_oligogenic, mixture_props,
                                         non_sparse_indices,
                                         max_oligogenic_effect = 0.12,
@@ -166,7 +166,7 @@ simulate_oligogenic_effects <- function(G, h2_oligogenic, n_oligogenic, mixture_
 #' @param h2_infinitesimal Heritability allocated to infinitesimal effects.
 #' @param infinitesimal_indices Indices of SNPs to receive infinitesimal effects.
 #' @return A vector of effect sizes (with zeros outside infinitesimal_indices).
-#' @export
+#' @keywords internal
 simulate_infinitesimal_effects <- function(G, h2_infinitesimal, infinitesimal_indices) {
   n_features <- ncol(G)
   beta <- rep(0, n_features)
@@ -193,7 +193,7 @@ simulate_infinitesimal_effects <- function(G, h2_infinitesimal, infinitesimal_in
 #' @details The function calculates power using a non-central chi-square distribution
 #'   with non-centrality parameter NCP = n * beta^2 * var(X) / residual_variance, where the
 #'   significance threshold is Bonferroni-corrected (alpha = 0.05 / p).
-#' @export
+#' @keywords internal
 is_causal_power <- function(G, beta, residual_variance, power = 0.80) {
   n <- nrow(G)
   p <- ncol(G)

R/simulate_genotype.R

@@ -1,13 +1,10 @@
-
-# sim_geno_functions.R
-
 #' Simulate Independent Genotypes
 #'
 #' @param n Sample size.
 #' @param p Number of SNPs.
 #' @param min_maf Minimum minor allele frequency.
 #' @param max_maf Maximum minor allele frequency.
-#' @param scale Logical, whether to scale the data using scale_faster.
+#' @param scale Logical, whether to scale the data.
 #' @return A matrix of genotypes.
 #' @examples
 #' sim_geno_indep(n = 100, p = 10, min_maf = 0.01, max_maf = 0.4, scale = TRUE)
@@ -22,7 +19,7 @@ sim_geno_indep <- function(n, p, min_maf = 0.01, max_maf = 0.4, scale = FALSE) {
   colnames(G) <- paste0("variant", 1:p)
 
   if (scale) {
-    G <- scale_faster(G)
+    G <- scale(G)
   }
 
   return(G)
@@ -35,13 +32,14 @@ sim_geno_indep <- function(n, p, min_maf = 0.01, max_maf = 0.4, scale = FALSE) {
 #' @param min_maf Minimum minor allele frequency.
 #' @param max_maf Maximum minor allele frequency.
 #' @param lambda Regularization parameter.
-#' @param type Type of genotype to simulate, either "continuous" or "discrete".
-#' @param scale Logical, whether to scale the data using scale_faster.
+#' @param is.discrete Logical, whether to generate discrete (TRUE) or continuous (FALSE) genotypes.
+#' @param scale Logical, whether to scale the data.
+#' @param tol Tolerance for checking diagonal elements of LD matrix.
 #' @return A matrix of genotypes.
 #' @examples
-#' sim_geno_LD(n = 100, LD = matrix(runif(100), 10, 10), min_maf = 0.01, max_maf = 0.4, lambda = 1e-3, type = "continuous", scale = TRUE)
+#' sim_geno_LD(n = 100, LD = matrix(runif(100), 10, 10), min_maf = 0.01, max_maf = 0.4, lambda = 1e-3, is.discrete = FALSE, scale = TRUE)
 #' @export
-sim_geno_LD <- function(n, LD, min_maf = 0.01, max_maf = 0.4, lambda = 1e-3, is.discrete = FALSE, scale = FALSE) {
+sim_geno_LD <- function(n, LD, min_maf = 0.01, max_maf = 0.4, lambda = 1e-3, is.discrete = FALSE, scale = FALSE, tol = sqrt(.Machine$double.eps)) {
   if (missing(n)) stop("Please provide the sample size")
   if (is.null(LD)) stop("Please provide LD matrix!")
 
@@ -58,7 +56,7 @@ sim_geno_LD <- function(n, LD, min_maf = 0.01, max_maf = 0.4, lambda = 1e-3, is.
   }
 
   if (scale) {
-    G <- scale_faster(G)
+    G <- scale(G)
   }
 
   colnames(G) <- paste0("variant", 1:p)
@@ -70,7 +68,7 @@ sim_geno_LD <- function(n, LD, min_maf = 0.01, max_maf = 0.4, lambda = 1e-3, is.
 #' @param n Sample size.
 #' @param file_path Path to the UK Biobank file.
 #' @param min_maf Minimum minor allele frequency.
-#' @param scale Logical, whether to scale the data using scale_faster.
+#' @param scale Logical, whether to scale the data.
 #' @return A matrix of genotypes.
 #' @examples
 #' sim_geno_UKB(n = 100, file_path = "path/to/UKB/file", min_maf = 0.01, scale = TRUE)
@@ -83,26 +81,31 @@ sim_geno_UKB <- function(n, file_path, min_maf = 0.01, scale = FALSE) {
   G <- process_ukb(file_path, n, min_maf)
 
   if (scale) {
-    G <- scale_faster(G)
+    G <- scale(G)
   }
 
-  colnames(G) <- paste0("variant", 1:p)
+  colnames(G) <- paste0("variant", 1:ncol(G))
   return(G)
 }
 
+# Safe multivariate normal random generation
+#' @importFrom mvtnorm rmvnorm
+#' @keywords internal
 safe_rmvnorm <- function(n, p, LD, maf, lambda = 1e-3) {
 
     var_g <- diag(sqrt(2*maf*(1-maf)))
     Sigma <- var_g %*% LD %*% var_g
-    G <- try(mvtnorm::rmvnorm(n, mean = rep(0, p), sigma = Sigma), silent = TRUE)
+    G <- try(rmvnorm(n, mean = rep(0, p), sigma = Sigma), silent = TRUE)
     if(any(class(G) == "try-error")) {
         Sigma <- get_ld_pd(Sigma, lambda = lambda)
-        G <- mvtnorm::rmvnorm(n, mean = rep(0, p), sigma = Sigma)
+        G <- rmvnorm(n, mean = rep(0, p), sigma = Sigma)
     }
     return(G)
 
 }
 
+# Convert continuous genotypes to discrete (binomial)
+#' @keywords internal
 binomialize_genotype <- function(G, p, maf) {
 
     sapply(1:p, function(i) {
@@ -116,9 +119,13 @@ binomialize_genotype <- function(G, p, maf) {
 
 }
 
+# Process UK Biobank genotype file
+#' @importFrom BEDMatrix BEDMatrix
+#' @importFrom data.table fread
+#' @keywords internal
 process_ukb <- function(file.path, n, min_maf = 0.01) {
 
-    G <- BEDMatrix::BEDMatrix(file.path)
+    G <- BEDMatrix(file.path)
     G <- data.matrix(G[1:n,])
     G <- apply(G, 2, function(g) {
         tmp <- which(is.na(g))
@@ -127,7 +134,7 @@ process_ukb <- function(file.path, n, min_maf = 0.01) {
     })
     maf <- colMeans(G, na.rm = TRUE) / 2
     snpname <- strsplit(file.path, split = "\\.bed")[[1]][1]
-    snp.names <- unlist(data.table::fread(paste0(snpname, ".bim"))[,2])
+    snp.names <- unlist(fread(paste0(snpname, ".bim"))[,2])
     colnames(G) <- snp.names
     poss <- which(maf <= min_maf)
     if (length(poss) != 0) {
@@ -138,6 +145,8 @@ process_ukb <- function(file.path, n, min_maf = 0.01) {
 }
 
 
+# Generate SNP with specified LD
+#' @keywords internal
 generate_snp_LD <- function(g, ld, maf) {
 
     # - change genotype to haplotype