update tests and documentations

xinhe-lab · Jul 11, 2024 · 3edb693 · 3edb693
1 parent a859278
commit 3edb693
Show file tree

Hide file tree

Showing 26 changed files with 346 additions and 613 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -4,7 +4,7 @@ Type: Package
 Title: Adjusting for Genetic Confounders in Transcriptome-Wide Association
   Studies Improves Discovery of Risk Genes of Complex Traits
 Date: 2024-06-14
-Version: 0.3.5
+Version: 0.3.6
 Authors@R: c(person("Siming","Zhao",role="aut",email="siming.zhao06@gmail.com"),
              person("Wesley","Crouse",role="aut"),
              person("Sheng","Qian",role="aut",email="shengqian@uchicago.edu"),

diff --git a/R/ctwas_finemapping.R b/R/ctwas_finemapping.R
@@ -111,6 +111,8 @@ finemap_region <- function(region_data,
   rm(res)
 
   if (!use_LD) {
+    # set L = 1 in no-LD version
+    L <- 1
     # use an identity matrix as R in no-LD version
     R <- diag(length(z))
     # do not include cs_index in no-LD version

diff --git a/R/ctwas_merge_regions.R b/R/ctwas_merge_regions.R
@@ -45,6 +45,11 @@
 #'
 #' @param cor_dir a string, the directory to store correlation (R) matrices
 #'
+#' @param LD_format file format for LD matrix. If "custom", use a user defined
+#' \code{LD_loader_fun()} function to load LD matrix.
+#'
+#' @param LD_loader_fun a user defined function to load LD matrix when \code{LD_format = "custom"}.
+#'
 #' @param ncore The number of cores used to parallelize computation over regions
 #'
 #' @param verbose TRUE/FALSE. If TRUE, print detail messages

diff --git a/R/ctwas_summarize_finemap_res.R b/R/ctwas_summarize_finemap_res.R
@@ -70,7 +70,7 @@ anno_finemap_res <- function(finemap_res,
   }
 
   # extract gene ids
-  finemap_gene_res <- finemap_res[finemap_res$type!="SNP",]
+  finemap_gene_res <- finemap_res[finemap_res$group!="SNP",]
 
   if (is.null(finemap_gene_res$gene_id)) {
     finemap_gene_res$gene_id <- sapply(strsplit(finemap_gene_res$id, split = "[|]"), "[[", 1)
@@ -115,7 +115,7 @@ anno_finemap_res <- function(finemap_res,
 
   # add SNP positions
   loginfo("add SNP positions from snp_info")
-  finemap_snp_res <- finemap_res[finemap_res$type=="SNP",]
+  finemap_snp_res <- finemap_res[finemap_res$group=="SNP",]
   snp_idx <- match(finemap_snp_res$id, snp_info$id)
   finemap_snp_res$chrom <- snp_info$chrom[snp_idx]
   finemap_snp_res$chrom <- parse_number(as.character(finemap_snp_res$chrom))
@@ -142,6 +142,7 @@ anno_finemap_res <- function(finemap_res,
 #'
 #' @export
 get_gene_annot_from_ens_db <- function(ens_db, gene_ids) {
+
   gene_ids <- unique(na.omit(gene_ids))
   if (any(grep("[.]", gene_ids))) {
     gene_ids_trimmed <- sapply(strsplit(gene_ids, split = "[.]"), "[[", 1)

diff --git a/README.md b/README.md
@@ -26,29 +26,30 @@ Running a cTWAS analysis involves four main steps:
 
 1. Preparing the input data. 
 
-2. Computing the gene z-scores. 
+2. Computing associations of molecular traits with the phenotype (Z-scores). 
 
 3. Estimating the model parameters. 
 
-4. Fine-mapping. 
+4. Fine-mapping causal molecular traits 
 
 The outputs of cTWAS are posterior inclusion probabilities (PIPs) for all variants and molecular traits.
 
 To learn more about the ctwas R package, we recommend starting with this introductory tutorial: 
 
-[A minimal tutorial of how to run cTWAS without LD](https://xinhe-lab.github.io/multigroup_ctwas/articles/simple_ctwas_tutorial.html) 
+[A minimal tutorial of how to run cTWAS without LD](https://xinhe-lab.github.io/singlegroup_ctwas
+/articles/minimal_ctwas_tutorial.html) 
 
 To run the full cTWAS, a few more tutorials including:
 
-- [Preparing input data](https://xinhe-lab.github.io/multigroup_ctwas/articles/preparing_ctwas_input_data.html) 
+- [Preparing input data](https://xinhe-lab.github.io/singlegroup_ctwas/articles/preparing_ctwas_input_data.html) 
 
-- [Running cTWAS main function](https://xinhe-lab.github.io/multigroup_ctwas/articles/ctwas_main_function.html)
+- [Running cTWAS main function](https://xinhe-lab.github.io/singlegroup_ctwas/articles/ctwas_main_function.html)
 
-- [Running cTWAS by modules](https://xinhe-lab.github.io/multigroup_ctwas/articles/ctwas_modules.html)
+- [Running cTWAS by modules](https://xinhe-lab.github.io/singlegroup_ctwas/articles/ctwas_modules.html)
 
-- [Summarizing and visualizing cTWAS results](https://xinhe-lab.github.io/multigroup_ctwas/articles/summarizing_ctwas_results.html)
+- [Summarizing and visualizing cTWAS results](https://xinhe-lab.github.io/singlegroup_ctwas/articles/summarizing_ctwas_results.html)
 
-- [Post-processing cTWAS results](https://xinhe-lab.github.io/multigroup_ctwas/articles/postprocessing_ctwas_results.html)
+- [Post-processing cTWAS results](https://xinhe-lab.github.io/singlegroup_ctwas/articles/postprocessing_ctwas_results.html)
 
 You can [browse source code](https://github.com/xinhe-lab/ctwas/tree/single_group) and [report a bug](https://github.com/xinhe-lab/ctwas/issues) here. 
 

diff --git a/...testthat/LDL_example.ctwas_parameters.RDS → ...ple_data/LDL_example.ctwas_parameters.RDS b/...testthat/LDL_example.ctwas_parameters.RDS → ...ple_data/LDL_example.ctwas_parameters.RDS
diff --git a/...t/LDL_example.ctwas_sumstats_noLD_res.RDS → ...a/LDL_example.ctwas_sumstats_noLD_res.RDS b/...t/LDL_example.ctwas_sumstats_noLD_res.RDS → ...a/LDL_example.ctwas_sumstats_noLD_res.RDS
diff --git a/...stthat/LDL_example.ctwas_sumstats_res.RDS → ...e_data/LDL_example.ctwas_sumstats_res.RDS b/...stthat/LDL_example.ctwas_sumstats_res.RDS → ...e_data/LDL_example.ctwas_sumstats_res.RDS
diff --git a/man/merge_finemap_regions.Rd b/man/merge_finemap_regions.Rd
diff --git a/tests/testthat/LDL_example.annotated_finemap_res.RDS b/tests/testthat/LDL_example.annotated_finemap_res.RDS
diff --git a/tests/testthat/test-ctwas_est_parameters.R b/tests/testthat/test-ctwas_est_parameters.R
@@ -1,6 +1,7 @@
 test_that("est_param works", {
 
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_noLD_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
+
   region_data <- ctwas_res$region_data
   expected_param <- ctwas_res$param
 

diff --git a/tests/testthat/test-ctwas_finemapping.R b/tests/testthat/test-ctwas_finemapping.R
@@ -2,7 +2,8 @@ test_that("finemap_regions works with no-LD", {
 
   # finemap a single region with no-LD version
   weights <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.weights.RDS", package = "ctwas"))
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_noLD_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
+
   precomputed_finemap_res <- ctwas_res$finemap_res
   screened_region_data <- ctwas_res$screened_region_data
   param <- ctwas_res$param
@@ -34,7 +35,8 @@ test_that("finemap_regions works with LD", {
   snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas"))
   weights <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.weights.RDS", package = "ctwas"))
 
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas"))
+
   precomputed_finemap_res <- ctwas_res$finemap_res
   screened_region_data <- ctwas_res$screened_region_data
   L <- ctwas_res$L

diff --git a/tests/testthat/test-ctwas_get_region_cor.R b/tests/testthat/test-ctwas_get_region_cor.R
@@ -24,7 +24,8 @@ test_that("get_region_cor correctly computes correlation matrices", {
   region_info <- readRDS(system.file("extdata/sample_data", "LDL_example.region_info.RDS", package = "ctwas"))
   snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas"))
   weights <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.weights.RDS", package = "ctwas"))
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas"))
+
   screened_region_data <- ctwas_res$screened_region_data
   # test computing the correlation matrices
   cor_res <- get_region_cor(region_id,

diff --git a/tests/testthat/test-ctwas_preprocess_regions.R b/tests/testthat/test-ctwas_preprocess_regions.R
@@ -16,9 +16,7 @@ test_that("preprocess_region_LD_snp_info (no-LD version) works", {
   region_info <- res$region_info
   snp_info <- res$snp_info
 
-  # Load region info
   expected_region_info <- readRDS(system.file("extdata/sample_data", "LDL_example.region_info.RDS", package = "ctwas"))
-  # Load SNP info
   expected_snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas"))
 
   expect_equal(region_info, expected_region_info)

diff --git a/tests/testthat/test-ctwas_region_data.R b/tests/testthat/test-ctwas_region_data.R
@@ -11,7 +11,8 @@ test_that("assemble_region_data works", {
   # Load gene z-scores
   z_gene <- readRDS(system.file("extdata/sample_data", "LDL_example.z_gene.RDS", package = "ctwas"))
 
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_noLD_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
+
   preassembled_region_data <- ctwas_res$region_data
   preassembled_boundary_genes <- ctwas_res$boundary_genes
 

diff --git a/tests/testthat/test-ctwas_screen_regions.R b/tests/testthat/test-ctwas_screen_regions.R
@@ -5,7 +5,7 @@ test_that("screen_regions (no-LD version) work", {
   z_snp <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.z_snp.RDS", package = "ctwas"))
   z_gene <- readRDS(system.file("extdata/sample_data", "LDL_example.z_gene.RDS", package = "ctwas"))
 
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_noLD_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
   region_data <- ctwas_res$region_data
   param <- ctwas_res$param
   group_prior <- param$group_prior
@@ -45,7 +45,8 @@ test_that("screen_regions (LD version) work", {
   z_snp <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.z_snp.RDS", package = "ctwas"))
   z_gene <- readRDS(system.file("extdata/sample_data", "LDL_example.z_gene.RDS", package = "ctwas"))
 
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas"))
+
   region_data <- ctwas_res$region_data
   param <- ctwas_res$param
   group_prior <- param$group_prior

diff --git a/tests/testthat/test-ctwas_summarize_finemap_res.R b/tests/testthat/test-ctwas_summarize_finemap_res.R
@@ -1,19 +1,19 @@
-# test_that("anno_finemap_res works", {
-#
-#   snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas"))
-#   ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas"))
-#   finemap_res <- ctwas_res$finemap_res
-#   gene_annot <- readRDS(system.file("extdata/sample_data", "LDL_example.gene_annot.RDS", package = "ctwas"))
-#
-#   expected_annotated_finemap_res <- readRDS("LDL_example.annotated_finemap_res.RDS")
-#
-#   capture.output({
-#     annotated_finemap_res <- anno_finemap_res(finemap_res,
-#                                               snp_info,
-#                                               gene_annot,
-#                                               use_gene_pos = "mid")
-#   })
-#
-#   expect_equal(annotated_finemap_res, expected_annotated_finemap_res)
-#
-# })
+test_that("anno_finemap_res works", {
+
+  snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas"))
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
+  finemap_res <- ctwas_res$finemap_res
+  gene_annot <- readRDS(system.file("extdata/sample_data", "LDL_example.gene_annot.RDS", package = "ctwas"))
+
+  expected_annotated_finemap_res <- readRDS("LDL_example.annotated_finemap_res.RDS")
+
+  capture.output({
+    annotated_finemap_res <- anno_finemap_res(finemap_res,
+                                              snp_info,
+                                              gene_annot,
+                                              use_gene_pos = "mid")
+  })
+
+  expect_equal(annotated_finemap_res, expected_annotated_finemap_res)
+
+})
diff --git a/tests/testthat/test-ctwas_summarize_parameters.R b/tests/testthat/test-ctwas_summarize_parameters.R
@@ -1,11 +1,13 @@
 test_that("summarize_param works", {
 
   gwas_n <- 343621
-  ctwas_res <- readRDS("LDL_example.ctwas_sumstats_noLD_res.RDS")
+  ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
   param <- ctwas_res$param
-  precomputed_ctwas_parameters <- readRDS("LDL_example.ctwas_parameters.RDS")
+
+  precomputed_ctwas_parameters <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_parameters.RDS", package = "ctwas"))
 
   ctwas_parameters <- summarize_param(param, gwas_n)
+
   expect_equal(ctwas_parameters, precomputed_ctwas_parameters)
 
 })
diff --git a/tests/testthat/test-ctwas_sumstats.R b/tests/testthat/test-ctwas_sumstats.R
@@ -11,7 +11,7 @@ test_that("ctwas_sumstats works", {
   region_info <- readRDS(system.file("extdata/sample_data", "LDL_example.region_info.RDS", package = "ctwas"))
   snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas"))
 
-  precomputed_ctwas_res <- readRDS("LDL_example.ctwas_sumstats_res.RDS")
+  precomputed_ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas"))
 
   capture.output({
     ctwas_res <- ctwas_sumstats(z_snp,

diff --git a/tests/testthat/test-ctwas_sumstats_noLD.R b/tests/testthat/test-ctwas_sumstats_noLD.R
@@ -5,7 +5,7 @@ test_that("ctwas_sumstats_noLD works", {
   z_snp <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.z_snp.RDS", package = "ctwas"))
   weights <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.weights.RDS", package = "ctwas"))
 
-  precomputed_ctwas_res <- readRDS("LDL_example.ctwas_sumstats_noLD_res.RDS")
+  precomputed_ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas"))
 
   capture.output({
     ctwas_res <- ctwas_sumstats_noLD(z_snp,