Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
The table of contents is too big for display.
Diff view
Diff view
  •  
  •  
  •  
416 changes: 209 additions & 207 deletions NAMESPACE

Large diffs are not rendered by default.

466 changes: 233 additions & 233 deletions R/AllClasses.R

Large diffs are not rendered by default.

172 changes: 86 additions & 86 deletions R/AllGenerics.R

Large diffs are not rendered by default.

502 changes: 251 additions & 251 deletions R/AllMethods.R

Large diffs are not rendered by default.

1,122 changes: 561 additions & 561 deletions R/LD.R

Large diffs are not rendered by default.

269 changes: 135 additions & 134 deletions R/allele_qc.R

Large diffs are not rendered by default.

1,484 changes: 742 additions & 742 deletions R/colocboost_pipeline.R

Large diffs are not rendered by default.

133 changes: 67 additions & 66 deletions R/compute_qtl_enrichment.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,26 +9,26 @@
#' @details Uses output of \code{\link[susieR]{susie}} from the
#' \code{susieR} package.
#'
#' @param gwas_pip This is a vector of GWAS PIP, genome-wide.
#' @param susie_qtl_regions This is a list of SuSiE fitted objects per QTL unit analyzed
#' @param num_gwas This parameter is highly important if GWAS input does not contain all SNPs interrogated (e.g., in some cases, only fine-mapped geomic regions are included).
#' Then users must pick a value of total_variants and estimate pi_gwas beforehand by: sum(gwas_pip$pip)/num_gwas. If num_gwas is null, pi_gwas would be sum(gwas_pip$pip)/total_variants.
#' @param pi_qtl This parameter can be safely left to default if your input QTL data has enough regions to estimate it.
#' @param gwasPip This is a vector of GWAS PIP, genome-wide.
#' @param susieQtlRegions This is a list of SuSiE fitted objects per QTL unit analyzed
#' @param numGwas This parameter is highly important if GWAS input does not contain all SNPs interrogated (e.g., in some cases, only fine-mapped geomic regions are included).
#' Then users must pick a value of total_variants and estimate piGwas beforehand by: sum(gwasPip$pip)/numGwas. If numGwas is null, piGwas would be sum(gwasPip$pip)/total_variants.
#' @param piQtl This parameter can be safely left to default if your input QTL data has enough regions to estimate it.
#' @param lambda Similar to the shrinkage parameter used in ridge regression. It takes any non-negative value and shrinks the enrichment estimate towards 0.
#' When it is set to 0, no shrinkage will be applied. A large value indicates strong shrinkage. The default value is set to 1.0.
#' @param ImpN Rounds of multiple imputation to draw QTL from, default is 25.
#' @param num_threads Number of Simultaneous running CPU threads for multiple imputation, default is 1.
#' @param impN Rounds of multiple imputation to draw QTL from, default is 25.
#' @param numThreads Number of Simultaneous running CPU threads for multiple imputation, default is 1.
#' @return A list of enrichment parameter estimates
#'
#' @examples
#'
#' # Simulate fake data for gwas_pip
#' n_gwas_pip <- 1000
#' gwas_pip <- runif(n_gwas_pip)
#' names(gwas_pip) <- paste0("snp", 1:n_gwas_pip)
#' gwas_fit <- list(pip = gwas_pip)
#' # Simulate fake data for gwasPip
#' nGwasPip <- 1000
#' gwasPip <- runif(nGwasPip)
#' names(gwasPip) <- paste0("snp", 1:nGwasPip)
#' gwasFit <- list(pip = gwasPip)
#' # Simulate fake data for a single SuSiEFit object
#' simulate_susiefit <- function(n, p) {
#' simulateSusiefit <- function(n, p) {
#' pip <- runif(n)
#' names(pip) <- paste0("snp", 1:n)
#' alpha <- t(matrix(runif(n * p), nrow = n))
Expand All @@ -40,95 +40,96 @@
#' )
#' }
#' # Simulate multiple SuSiEFit objects
#' n_susie_fits <- 2
#' susie_fits <- replicate(n_susie_fits, simulate_susiefit(n_gwas_pip, 10), simplify = FALSE)
#' nSusieFits <- 2
#' susieFits <- replicate(nSusieFits, simulateSusiefit(nGwasPip, 10), simplify = FALSE)
#' # Add these fits to a list, providing names to each element
#' names(susie_fits) <- paste0("fit", 1:length(susie_fits))
#' names(susieFits) <- paste0("fit", 1:length(susieFits))
#' # Set other parameters
#' ImpN <- 10
#' impN <- 10
#' lambda <- 1
#' num_threads <- 1
#' numThreads <- 1
#' library(pecotmr)
#' en <- compute_qtl_enrichment(gwas_fit, susie_fits, lambda = lambda, ImpN = ImpN, num_threads = num_threads)
#' en <- computeQtlEnrichment(gwasFit, susieFits, lambda = lambda, impN = impN, numThreads = numThreads)
#'
#' @seealso \code{\link[susieR]{susie}}
#' @useDynLib pecotmr, .registration = TRUE
#' @export
#'
compute_qtl_enrichment <- function(gwas_pip, susie_qtl_regions,
num_gwas = NULL, pi_qtl = NULL,
lambda = 1.0, ImpN = 25,
double_shrinkage = FALSE,
bessel_correction = TRUE,
num_threads = 1, verbose = TRUE) {
if (is.null(num_gwas)) {
warning("num_gwas is not provided. Estimating pi_gwas from the data. Note that this estimate may be biased if the input gwas_pip does not contain genome-wide variants.")
pi_gwas <- sum(gwas_pip) / length(gwas_pip)
computeQtlEnrichment <- function(gwasPip, susieQtlRegions,
numGwas = NULL, piQtl = NULL,
lambda = 1.0, impN = 25,
doubleShrinkage = FALSE,
besselCorrection = TRUE,
numThreads = 1, verbose = TRUE) {
if (is.null(numGwas)) {
warning("numGwas is not provided. Estimating piGwas from the data. Note that this estimate may be biased if the input gwasPip does not contain genome-wide variants.")
piGwas <- sum(gwasPip) / length(gwasPip)
if (verbose) {
message(paste("Estimated pi_gwas: ", round(pi_gwas, 5), "\n"))
message(paste("Estimated piGwas: ", round(piGwas, 5), "\n"))
}
} else {
pi_gwas <- sum(gwas_pip) / num_gwas
piGwas <- sum(gwasPip) / numGwas
}

if (is.null(pi_qtl)) {
warning("pi_qtl is not provided. Estimating pi_qtl from the data. Note that this estimate may be biased if either 1) the input susie_qtl_regions does not have enough data, or 2) the single effects only include variables inside of credible sets or signal clusters.")
num_signal <- 0
num_test <- 0
for (d in susie_qtl_regions) {
num_signal <- num_signal + sum(d$pip)
num_test <- num_test + length(d$pip)
if (is.null(piQtl)) {
warning("piQtl is not provided. Estimating piQtl from the data. Note that this estimate may be biased if either 1) the input susieQtlRegions does not have enough data, or 2) the single effects only include variables inside of credible sets or signal clusters.")
numSignal <- 0
numTest <- 0
for (d in susieQtlRegions) {
numSignal <- numSignal + sum(d$pip)
numTest <- numTest + length(d$pip)
}
pi_qtl <- num_signal / num_test
piQtl <- numSignal / numTest
if (verbose) {
message(paste("Estimated pi_qtl: ", round(pi_qtl, 5), "\n"))
message(paste("Estimated piQtl: ", round(piQtl, 5), "\n"))
}
}

if (pi_gwas == 0) stop("Cannot perform enrichment analysis. No association signal found in GWAS data.")
if (pi_qtl == 0) stop("Cannot perform enrichment analysis. No QTL associated with the molecular phenotype.")
if (piGwas == 0) stop("Cannot perform enrichment analysis. No association signal found in GWAS data.")
if (piQtl == 0) stop("Cannot perform enrichment analysis. No QTL associated with the molecular phenotype.")

# Check if names of gwas_pip and susie_qtl_regions$pip are both available
if (is.null(names(gwas_pip))) {
stop("Variant names are missing in gwas_pip. Please provide named gwas_pip data.")
# Check if names of gwasPip and susieQtlRegions$pip are both available
if (is.null(names(gwasPip))) {
stop("Variant names are missing in gwasPip. Please provide named gwasPip data.")
}
if (!all(sapply(susie_qtl_regions, function(x) !is.null(names(x$pip))))) {
stop("Variant names are missing in susie_qtl_regions$pip. Please provide susie_qtl_regions with named pip data.")
if (!all(sapply(susieQtlRegions, function(x) !is.null(names(x$pip))))) {
stop("Variant names are missing in susieQtlRegions$pip. Please provide susieQtlRegions with named pip data.")
}

# Align the names of susie_qtl_regions$pip to gwas_pip names and document unmatched variants
aligned_susie_qtl_regions <- lapply(susie_qtl_regions, function(x) {
alignment_result <- align_variant_names(names(x$pip), names(gwas_pip))
names(x$pip) <- alignment_result$aligned_variants
if (length(alignment_result$unmatched_indices) > 0) {
x$unmatched_variants <- names(x$pip)[alignment_result$unmatched_indices]
# Align the names of susieQtlRegions$pip to gwasPip names and document unmatched variants
alignedSusieQtlRegions <- lapply(susieQtlRegions, function(x) {
alignmentResult <- alignVariantNames(names(x$pip), names(gwasPip))
names(x$pip) <- alignmentResult$aligned_variants
if (length(alignmentResult$unmatched_indices) > 0) {
x$unmatched_variants <- names(x$pip)[alignmentResult$unmatched_indices]
}
x
})
unmatched_variants <- lapply(aligned_susie_qtl_regions, function(x) x$unmatched_variants)
unmatchedVariants <- lapply(alignedSusieQtlRegions, function(x) x$unmatched_variants)

# Update susie_qtl_regions with the aligned variant names
susie_qtl_regions <- lapply(aligned_susie_qtl_regions, function(x) {
# Update susieQtlRegions with the aligned variant names
susieQtlRegions <- lapply(alignedSusieQtlRegions, function(x) {
x$unmatched_variants <- NULL
x
})

# cpp11 requires exact integer types for int parameters
en <- qtl_enrichment_rcpp(
r_gwas_pip = gwas_pip,
r_qtl_susie_fit = susie_qtl_regions,
pi_gwas = pi_gwas,
pi_qtl = pi_qtl,
ImpN = as.integer(ImpN),
shrinkage_lambda = lambda,
double_shrinkage = double_shrinkage,
bessel_correction = bessel_correction,
num_threads = as.integer(num_threads)
en <- qtlEnrichmentRcpp(
rGwasPip = gwasPip,
rQtlSusieFit = susieQtlRegions,
piGwas = piGwas,
piQtl = piQtl,
ImpN = as.integer(impN),
shrinkageLambda = lambda,
doubleShrinkage = doubleShrinkage,
besselCorrection = besselCorrection,
numThreads = as.integer(numThreads)
)

# Add the unmatched variants to the output
en <- list(en)
en$unused_xqtl_variants <- unmatched_variants
en$unused_xqtl_variants <- unmatchedVariants

return(en)
}

24 changes: 12 additions & 12 deletions R/cpp11.R
Original file line number Diff line number Diff line change
@@ -1,25 +1,25 @@
# Generated by cpp11: do not edit by hand

dentist_iterative_impute <- function(LD_mat_r, nSample, zScore_r, pValueThreshold, propSVD, gcControl, nIter, gPvalueThreshold, ncpus, correct_chen_et_al_bug, verbose) {
.Call(`_pecotmr_dentist_iterative_impute`, LD_mat_r, nSample, zScore_r, pValueThreshold, propSVD, gcControl, nIter, gPvalueThreshold, ncpus, correct_chen_et_al_bug, verbose)
dentistIterativeImpute <- function(ldMatR, nSample, zScoreR, pValueThreshold, propSVD, gcControl, nIter, gPvalueThreshold, ncpus, correctChenEtAlBug, verbose) {
.Call(`_pecotmr_dentistIterativeImpute`, ldMatR, nSample, zScoreR, pValueThreshold, propSVD, gcControl, nIter, gPvalueThreshold, ncpus, correctChenEtAlBug, verbose)
}

lassosum_rss_rcpp <- function(z_r, LD, lambda_r, thr, maxiter) {
.Call(`_pecotmr_lassosum_rss_rcpp`, z_r, LD, lambda_r, thr, maxiter)
lassosumRssRcpp <- function(zR, LD, lambdaR, thr, maxiter) {
.Call(`_pecotmr_lassosumRssRcpp`, zR, LD, lambdaR, thr, maxiter)
}

penalized_rss_rcpp <- function(z_r, LD, lambda_r, penalty_str, gamma, alpha, lambda0, lambda2, thr, maxiter, max_swaps) {
.Call(`_pecotmr_penalized_rss_rcpp`, z_r, LD, lambda_r, penalty_str, gamma, alpha, lambda0, lambda2, thr, maxiter, max_swaps)
penalizedRssRcpp <- function(zR, LD, lambdaR, penaltyStr, gamma, alpha, lambda0, lambda2, thr, maxiter, maxSwaps) {
.Call(`_pecotmr_penalizedRssRcpp`, zR, LD, lambdaR, penaltyStr, gamma, alpha, lambda0, lambda2, thr, maxiter, maxSwaps)
}

prs_cs_rcpp <- function(a, b, phi, bhat, maf, n, ld_blk, n_iter, n_burnin, thin, verbose, seed) {
.Call(`_pecotmr_prs_cs_rcpp`, a, b, phi, bhat, maf, n, ld_blk, n_iter, n_burnin, thin, verbose, seed)
prsCsRcpp <- function(a, b, phi, bhat, maf, n, ldBlk, nIter, nBurnin, thin, verbose, seed) {
.Call(`_pecotmr_prsCsRcpp`, a, b, phi, bhat, maf, n, ldBlk, nIter, nBurnin, thin, verbose, seed)
}

qtl_enrichment_rcpp <- function(r_gwas_pip, r_qtl_susie_fit, pi_gwas, pi_qtl, ImpN, shrinkage_lambda, double_shrinkage, bessel_correction, num_threads) {
.Call(`_pecotmr_qtl_enrichment_rcpp`, r_gwas_pip, r_qtl_susie_fit, pi_gwas, pi_qtl, ImpN, shrinkage_lambda, double_shrinkage, bessel_correction, num_threads)
qtlEnrichmentRcpp <- function(rGwasPip, rQtlSusieFit, piGwas, piQtl, ImpN, shrinkageLambda, doubleShrinkage, besselCorrection, numThreads) {
.Call(`_pecotmr_qtlEnrichmentRcpp`, rGwasPip, rQtlSusieFit, piGwas, piQtl, ImpN, shrinkageLambda, doubleShrinkage, besselCorrection, numThreads)
}

sdpr_rcpp <- function(bhat_r, LD, n, per_variant_sample_size, array, a, c, M, a0k, b0k, iter, burn, thin, n_threads, opt_llk, verbose, seed) {
.Call(`_pecotmr_sdpr_rcpp`, bhat_r, LD, n, per_variant_sample_size, array, a, c, M, a0k, b0k, iter, burn, thin, n_threads, opt_llk, verbose, seed)
sdprRcpp <- function(bhatR, LD, n, perVariantSampleSize, array, a, c, M, a0k, b0k, iter, burn, thin, nThreads, optLlk, verbose, seed) {
.Call(`_pecotmr_sdprRcpp`, bhatR, LD, n, perVariantSampleSize, array, a, c, M, a0k, b0k, iter, burn, thin, nThreads, optLlk, verbose, seed)
}
Loading