minor fix

NAMESPACE

@@ -30,6 +30,7 @@ export(dentist)
 export(dentist_single_window)
 export(enet_weights)
 export(extract_cs_info)
+export(extract_flatten_sumstats_from_nested)
 export(extract_top_pip_info)
 export(filter_invalid_summary_stat)
 export(filter_mixture_components)
@@ -52,6 +53,7 @@ export(lbf_to_alpha)
 export(load_LD_matrix)
 export(load_genotype_region)
 export(load_ld_snp_info)
+export(load_multicontext_sumstats)
 export(load_multitask_regional_data)
 export(load_multitrait_R_sumstat)
 export(load_multitrait_tensorqtl_sumstat)
@@ -67,6 +69,7 @@ export(load_twas_weights)
 export(mash_pipeline)
 export(mash_rand_null_sample)
 export(merge_mash_data)
+export(merge_sumstats_matrices)
 export(mr_analysis)
 export(mr_ash_rss)
 export(mr_ash_rss_weights)
@@ -117,6 +120,9 @@ export(univariate_analysis_pipeline)
 export(wald_test_pval)
 export(xqtl_enrichment_wrapper)
 export(z_to_pvalue)
+import(dplyr)
+import(tibble)
+import(tidyr)
 import(vroom)
 importFrom(IRanges,IRanges)
 importFrom(IRanges,end)

R/file_utils.R

@@ -1049,6 +1049,7 @@ get_filter_lbf_index <- function(susie_obj, coverage = 0.5, size_factor = 0.5) {
 
   # Return filtered lbf_variable rows (one per CS)
   return(cs_index)
+}
 
 #' Function to load LD reference data variants
 #' @export

R/twas.R

@@ -165,7 +165,7 @@ harmonize_twas <- function(twas_weights_data, ld_meta_file_path, gwas_meta_file)
           )
           weights_matrix_subset <- as.matrix(weights_matrix_qced$target_data_qced[, !colnames(weights_matrix_qced$target_data_qced) %in% c(
             "chrom",
-            "pos", "A2", "A1", "variant_id"
+            "pos", "A2", "A1", "variant_id", "variants_id_original"
           ), drop = FALSE])
           rownames(weights_matrix_subset) <- weights_matrix_qced$target_data_qced$variant_id # weight variant names are flipped/corrected
 
@@ -259,6 +259,7 @@ harmonize_gwas <- function(gwas_file, query_region, ld_variants, col_to_flip=NUL
     if (!any(gwas_data_sumstats$pos %in% gsub("\\:.*$", "", sub("^.*?\\:", "", ld_variants)))) return(NULL)
     gwas_allele_flip <- allele_qc(gwas_data_sumstats, ld_variants, col_to_flip=col_to_flip, match_min_prop = match_min_prop)
     gwas_data_sumstats <- gwas_allele_flip$target_data_qced # post-qc gwas data that is flipped and corrected - gwas study level
+    gwas_data_sumstats <- gwas_data_sumstats[!is.na(gwas_data_sumstats$z) & !is.infinite(gwas_data_sumstats$z), ]
     return(gwas_data_sumstats)
 }
 
@@ -464,9 +465,9 @@ twas_pipeline <- function(twas_weights_data,
     # Nested lapply for contexts and gwas studies
     twas_gene_results <- lapply(contexts, function(context) {
       study_results <- lapply(gwas_studies, function(study) {
-        twas_variants <- intersect(
-          rownames(twas_data_qced[[weight_db]][["weights_qced"]][[context]][[study]][["weights"]]),
-          twas_data_qced[[weight_db]][["variant_names"]][[context]][[study]]
+        twas_variants <- Reduce(intersect, list(rownames(twas_data_qced[[weight_db]][["weights_qced"]][[context]][[study]][["weights"]]), 
+          twas_data_qced[[weight_db]][["variant_names"]][[context]][[study]],
+          twas_data_qced[[weight_db]][["gwas_qced"]][[study]]$variant_id)
         )
         if (length(twas_variants) == 0) {
           return(list(twas_rs_df = data.frame(), mr_rs_df = data.frame()))