From 295443470e37324c9716fd5b694747848317f2c9 Mon Sep 17 00:00:00 2001 From: Anjing Date: Sun, 6 Apr 2025 06:31:59 -0400 Subject: [PATCH 1/3] minor fix in qvalue: if all NA, then NA --- R/misc.R | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/R/misc.R b/R/misc.R index 05d3632e..bf014853 100644 --- a/R/misc.R +++ b/R/misc.R @@ -55,6 +55,10 @@ compute_qvalues <- function(pvalues) { if (!requireNamespace("qvalue", quietly = TRUE)) { stop("To use this function, please install qvalue: https://www.bioconductor.org/packages/release/bioc/html/qvalue.html") } + if (all(is.na(pvalues))) { + message("All p-values are NA. Returning NA vector.") + return(rep(NA_real_, length(pvalues))) + } tryCatch( { if (length(pvalues) < 2) { From 728b58157d8d31f3236537b6e9a9b57bba3e09df Mon Sep 17 00:00:00 2001 From: Anjing Date: Sun, 6 Apr 2025 06:34:19 -0400 Subject: [PATCH 2/3] twas:harmonize_twas, fixed ld_info_data col names issue, fixed missing results[[molecular_id]][["variant_names"]] structure issue --- R/twas.R | 34 ++++++++++++++++++++++++++++------ 1 file changed, 28 insertions(+), 6 deletions(-) diff --git a/R/twas.R b/R/twas.R index 9cb219b0..bce5c3ac 100644 --- a/R/twas.R +++ b/R/twas.R @@ -103,7 +103,8 @@ harmonize_twas <- function(twas_weights_data, ld_meta_file_path, gwas_meta_file) # function to extract LD variance for the query region query_variance <- function(ld_variant_info_data, extract_coordinates) { ld_info_data <- do.call(rbind, ld_variant_info_data) - ld_info_data_filtered <- ld_info_data[ld_info_data$V4 %in% extract_coordinates$pos, , drop = FALSE] + ld_pos <- as.numeric(ld_info_data[, 4]) + ld_info_data_filtered <- ld_info_data[ld_pos %in% extract_coordinates$pos, , drop = FALSE] variance_df <- ld_info_data_filtered[, c(1, 4, 5:7)] # Extract the variance column (7th column) colnames(variance_df) <- c("chrom", "pos", "A1", "A2", "variance") return(variance_df) @@ -154,6 +155,15 @@ harmonize_twas <- function(twas_weights_data, ld_meta_file_path, gwas_meta_file) for (molecular_id in molecular_ids) { results[[molecular_id]][["chrom"]] <- chrom results[[molecular_id]][["data_type"]] <- if ("data_type" %in% names(twas_weights_data[[molecular_id]])) twas_weights_data[[molecular_id]]$data_type + results[[molecular_id]][["variant_names"]] <- list() + for (context in names(twas_weights_data[[molecular_id]][["weights"]])) { + for (study in names(gwas_files)) { + if (is.null(results[[molecular_id]][["variant_names"]][[context]])) { + results[[molecular_id]][["variant_names"]][[context]] <- list() + } + results[[molecular_id]][["variant_names"]][[context]][[study]] <- rownames(twas_weights_data[[molecular_id]][["weights"]][[context]]) + } + } # group contexts based on the variant position context_clusters <- group_contexts_by_region(twas_weights_data[[molecular_id]], molecular_id, chrom, tolerance = 5000) @@ -226,7 +236,17 @@ harmonize_twas <- function(twas_weights_data, ld_meta_file_path, gwas_meta_file) paste0("chr", variant_qc$target_data_qced$variant_id[variant_qc$target_data_qced$variant_id %in% postqc_weight_variants]) }) } - rownames(weights_matrix_subset) <- if (!grepl("^chr", rownames(weights_matrix_subset)[1])) paste0("chr", rownames(weights_matrix_subset)) else rownames(weights_matrix_subset) + #rownames(weights_matrix_subset) <- if (!grepl("^chr", rownames(weights_matrix_subset)[1])) paste0("chr", rownames(weights_matrix_subset)) else rownames(weights_matrix_subset) + if (nrow(weights_matrix_subset) > 0) { + rownames(weights_matrix_subset) <- if (!grepl("^chr", rownames(weights_matrix_subset)[1])) { + paste0("chr", rownames(weights_matrix_subset)) + } else { + rownames(weights_matrix_subset) + } + } else { + warning("weights_matrix_subset is empty. Skipping this context.") + next + } results[[molecular_id]][["variant_names"]][[context]][[study]] <- rownames(weights_matrix_subset) # Step 6: scale weights by variance @@ -694,13 +714,15 @@ twas_analysis <- function(weights_matrix, gwas_sumstats_db, LD_matrix, extract_v gwas_sumstats_subset <- gwas_sumstats_db[match(extract_variants_objs, gwas_sumstats_db$variant_id), ] # Validate that the GWAS subset is not empty if (nrow(gwas_sumstats_subset) == 0 | all(is.na(gwas_sumstats_subset))) { - stop("No GWAS summary statistics found for the specified variants.") - } + warning("No GWAS summary statistics found for the specified variants.") + return(NULL) + } # Check if extract_variants_objs are in the rownames of LD_matrix valid_indices <- extract_variants_objs %in% rownames(LD_matrix) if (!any(valid_indices)) { - stop("None of the specified variants are present in the LD matrix.") - } + warning("None of the specified variants are present in the LD matrix. Skipping this context.") + return(NULL) + } # Extract only the valid indices from extract_variants_objs valid_variants_objs <- extract_variants_objs[valid_indices] # Extract LD_matrix subset using valid indices From 9d18931c8b10d598d4813207719b10234a1e6d5e Mon Sep 17 00:00:00 2001 From: Anjing Date: Tue, 8 Apr 2025 00:21:31 -0400 Subject: [PATCH 3/3] removed duplicate code variant_names --- R/twas.R | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/R/twas.R b/R/twas.R index bce5c3ac..bd386c01 100644 --- a/R/twas.R +++ b/R/twas.R @@ -156,14 +156,7 @@ harmonize_twas <- function(twas_weights_data, ld_meta_file_path, gwas_meta_file) results[[molecular_id]][["chrom"]] <- chrom results[[molecular_id]][["data_type"]] <- if ("data_type" %in% names(twas_weights_data[[molecular_id]])) twas_weights_data[[molecular_id]]$data_type results[[molecular_id]][["variant_names"]] <- list() - for (context in names(twas_weights_data[[molecular_id]][["weights"]])) { - for (study in names(gwas_files)) { - if (is.null(results[[molecular_id]][["variant_names"]][[context]])) { - results[[molecular_id]][["variant_names"]][[context]] <- list() - } - results[[molecular_id]][["variant_names"]][[context]][[study]] <- rownames(twas_weights_data[[molecular_id]][["weights"]][[context]]) - } - } + # group contexts based on the variant position context_clusters <- group_contexts_by_region(twas_weights_data[[molecular_id]], molecular_id, chrom, tolerance = 5000)