search

neurogenomics / MAGMA_Celltyping

Find causal cell-types underlying complex trait genetics
https://neurogenomics.github.io/MAGMA_Celltyping
71 stars 31 forks source link

`plot_celltype_associations`: Fix adding full method to results #121

Closed AMCalejandro closed 1 year ago

AMCalejandro commented 2 years ago

Hey Brian,

This sum will always retrieve 0, as far as I am aware. I think what you wanted to do is what I just added?

Also, note that MAGMA.Celltyping does not load ggplot2 into the workspace. Probably the user should do this!! ( my bad) But maybe it is safer to do ggplot2::theme_set

AMCalejandro commented 2 years ago

I was revisitting this and I thought I would be clearer showing the input data.

I do see you were expecting an unnammed list to store the results after running your example code. However, when I perform the association test, results are stored underl "levels" elements in the list

See here:

# Map GWAS snps to genes # This is basically the out file from
genesOutPath_intelligence <- MAGMA.Celltyping::map_snps_to_genes(
  path_formatted = gwas_sumstats_path,
  genome_build = "GRCh37",
  population = "eur",
  N = 111114,force_new = T)

# Loading transcriptome data from the Karolinska instititute
ctd <- MAGMA.Celltyping::get_ctd("ctd_allKI")
ctd_quant = 
  MAGMA.Celltyping::prepare_quantile_groups(ctd = ctd,
                                            non121_strategy = "drop_both_species",
                                            standardise = TRUE,
                                            input_species = "mouse",
                                            output_species = "human",
                                            numberOfBins = 40)

ctAssocsLinear <- MAGMA.Celltyping::calculate_celltype_associations(
  ctd = ctd_quant,
  gwas_sumstats_path = gwas_sumstats_path,
  magma_dir = NULL,
  genesOutCOND = NA,
  ctd_species = infer_ctd_species(ctd),
  force_new = FALSE)

> names(ctAssocsLinear)
[1] "level1"                         "level2"                         "total_baseline_tests_performed" "gwas_sumstats_path"            
[5] "analysis_name"                  "upstream_kb"                    "downstream_kb"        

 str(ctAssocsLinear$level1)
List of 1
 $ results:'data.frame':    24 obs. of  14 variables:
  ..$ Celltype      : chr [1:24] "astrocytes_ependymal" "Dopaminergic_Adult" "Dopaminergic_Neuroblast" "Embryonic_Dopaminergic_Neuron" ...
  ..$ OBS_GENES     : chr [1:24] "71" "71" "71" "71" ...
  ..$ BETA          : num [1:24] -0.005079 -0.004212 -0.000215 0.00292 -0.005705 ...
  ..$ BETA_STD      : num [1:24] -0.05943 -0.04793 -0.00273 0.03446 -0.0622 ...
  ..$ SE            : num [1:24] 0.00846 0.0085 0.00748 0.0082 0.00951 ...
  ..$ P             : num [1:24] 0.725 0.689 0.511 0.362 0.724 ...
  ..$ level         : int [1:24] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ Method        : chr [1:24] "MAGMA" "MAGMA" "MAGMA" "MAGMA" ...
  ..$ GCOV_FILE     : chr [1:24] "educational_attainment.tsv.35UP.10DOWN.level1.MainRun.gsa.out" "educational_attainment.tsv.35UP.10DOWN.level1.MainRun.gsa.out" "educational_attainment.tsv.35UP.10DOWN.level1.MainRun.gsa.out" "educational_attainment.tsv.35UP.10DOWN.level1.MainRun.gsa.out" ...
  ..$ CONTROL       : chr [1:24] "BASELINE" "BASELINE" "BASELINE" "BASELINE" ...
  ..$ CONTROL_label : chr [1:24] "BASELINE" "BASELINE" "BASELINE" "BASELINE" ...
  ..$ log10p        : num [1:24] -0.14 -0.162 -0.291 -0.442 -0.14 ...
  ..$ genesOutCOND  : chr [1:24] "NA" "NA" "NA" "NA" ...
  ..$ EnrichmentMode: chr [1:24] "Linear" "Linear" "Linear" "Linear" ...
bschilder commented 1 year ago

Hey Brian,

This sum will always retrieve 0, as far as I am aware. I think what you wanted to do is what I just added?

Based on the code that was edited, it looks like you're talking about the plot_celltype_associations function. I've updated the title of this PR to clarify.

Regarding the issue, looks like your PR does indeed fix this. Thanks for this! I'll add a unit test as well to ensure this gets checked.

Also, note that MAGMA.Celltyping does not load ggplot2 into the workspace. Probably the user should do this!! ( my bad) But maybe it is safer to do ggplot2::theme_set

Thanks for catching this, I'll make a fix.

Reprex

I've extended your code to make a fully reproducible example (eg one that i can run all the way through by copying and pasting it into my R console)

gwas_sumstats_path <-  MAGMA.Celltyping::get_example_gwas()

# Map GWAS snps to genes # This is basically the out file from
genesOutPath_intelligence <- MAGMA.Celltyping::map_snps_to_genes(
    path_formatted = gwas_sumstats_path,
    genome_build = "GRCh37",
    population = "eur",
    N = 111114,force_new = T)

# Loading transcriptome data from the Karolinska instititute
ctd <- MAGMA.Celltyping::get_ctd("ctd_allKI")
ctd_quant = 
  MAGMA.Celltyping::prepare_quantile_groups(ctd = ctd,
                                            non121_strategy = "drop_both_species",
                                            standardise = TRUE,
                                            input_species = "mouse",
                                            output_species = "human",
                                            numberOfBins = 40)

ctAssocsLinear <- MAGMA.Celltyping::calculate_celltype_associations(
  ctd = ctd_quant,
  gwas_sumstats_path = gwas_sumstats_path,
  magma_dir = NULL,
  genesOutCOND = NA,
  ctd_species = infer_ctd_species(ctd),
  force_new = FALSE)

### Passing in ctAssocsLinear directly ####
figs <- MAGMA.Celltyping::plot_celltype_associations(
    ctAssocs = ctAssocsLinear,
    ctd = ctd)
I

## Indeed this returns an empty list 
figs <- MAGMA.Celltyping::plot_celltype_associations(
    ctAssocs = ctAssocsLinear,
    ctd = ctd)

Session info

``` R version 4.2.1 (2022-06-23) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Monterey 12.4 Matrix products: default LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] MAGMA.Celltyping_2.0.8 loaded via a namespace (and not attached): [1] utf8_1.2.3 R.utils_2.12.2 tidyselect_1.2.0 [4] lme4_1.1-31 RSQLite_2.2.20 AnnotationDbi_1.60.0 [7] htmlwidgets_1.6.1 grid_4.2.1 BiocParallel_1.32.5 [10] munsell_0.5.0 codetools_0.2-19 withr_2.5.0 [13] colorspace_2.1-0 Biobase_2.58.0 filelock_1.0.2 [16] knitr_1.42 rstudioapi_0.14 orthogene_1.4.1 [19] stats4_4.2.1 SingleCellExperiment_1.20.0 ggsignif_0.6.4 [22] gitcreds_0.1.2 MatrixGenerics_1.10.0 GenomeInfoDbData_1.2.9 [25] bit64_4.0.5 rprojroot_2.0.3 vctrs_0.5.2 [28] treeio_1.22.0 generics_0.1.3 xfun_0.37 [31] timechange_0.2.0 BiocFileCache_2.6.0 R6_2.5.1 [34] GenomeInfoDb_1.34.9 bitops_1.0-7 cachem_1.0.6 [37] gridGraphics_0.5-1 DelayedArray_0.24.0 assertthat_0.2.1 [40] promises_1.2.0.1 BiocIO_1.8.0 scales_1.2.1 [43] gtable_0.3.1 rlang_1.0.6 MungeSumstats_1.6.0 [46] splines_4.2.1 rtracklayer_1.58.0 rstatix_0.7.2 [49] lazyeval_0.2.2 gargle_1.3.0 broom_1.0.3 [52] BiocManager_1.30.19 yaml_2.3.7 reshape2_1.4.4 [55] abind_1.4-5 GenomicFeatures_1.50.4 backports_1.4.1 [58] httpuv_1.6.8 tools_4.2.1 ggplotify_0.1.0 [61] ggplot2_3.4.0 ellipsis_0.3.2 RColorBrewer_1.1-3 [64] ggdendro_0.1.23 BiocGenerics_0.44.0 Rcpp_1.0.10 [67] plyr_1.8.8 progress_1.2.2 zlibbioc_1.44.0 [70] purrr_1.0.1 RCurl_1.98-1.10 prettyunits_1.1.1 [73] ggpubr_0.5.0 S4Vectors_0.36.1 SummarizedExperiment_1.28.0 [76] grr_0.9.5 fs_1.6.1 magrittr_2.0.3 [79] data.table_1.14.6 gh_1.3.1 matrixStats_0.63.0 [82] pkgload_1.3.2 hms_1.1.2 patchwork_1.1.2 [85] mime_0.12 evaluate_0.20 xtable_1.8-4 [88] XML_3.99-0.13 EWCE_1.6.0 IRanges_2.32.0 [91] gridExtra_2.3 compiler_4.2.1 biomaRt_2.54.0 [94] tibble_3.1.8 templateR_0.99.0 crayon_1.5.2 [97] minqa_1.2.5 R.oo_1.25.0 htmltools_0.5.4 [100] ggfun_0.0.9 later_1.3.0 tidyr_1.3.0 [103] aplot_0.1.9 lubridate_1.9.1 DBI_1.1.3 [106] ExperimentHub_2.6.0 gprofiler2_0.2.1 dbplyr_2.3.0 [109] MASS_7.3-58.2 rappdirs_0.3.3 boot_1.3-28.1 [112] babelgene_22.9 dlstats_0.1.6 Matrix_1.5-3 [115] car_3.1-1 piggyback_0.1.4 cli_3.6.0 [118] R.methodsS3_1.8.2 parallel_4.2.1 GenomicRanges_1.50.2 [121] pkgconfig_2.0.3 GenomicAlignments_1.34.0 plotly_4.10.1 [124] xml2_1.3.3 ggtree_3.6.2 XVector_0.38.0 [127] yulab.utils_0.0.6 stringr_1.5.0 VariantAnnotation_1.44.0 [130] digest_0.6.31 Biostrings_2.66.0 rmarkdown_2.20 [133] HGNChelper_0.8.1 tidytree_0.4.2 restfulr_0.0.15 [136] curl_5.0.0 shiny_1.7.4 Rsamtools_2.14.0 [139] rjson_0.2.21 nloptr_2.0.3 lifecycle_1.0.3 [142] nlme_3.1-162 jsonlite_1.8.4 carData_3.0-5 [145] desc_1.4.2 viridisLite_0.4.1 limma_3.54.1 [148] BSgenome_1.66.2 fansi_1.0.4 pillar_1.8.1 [151] lattice_0.20-45 homologene_1.4.68.19.3.27 KEGGREST_1.38.0 [154] fastmap_1.1.0 httr_1.4.4 googleAuthR_2.0.0 [157] interactiveDisplayBase_1.36.0 glue_1.6.2 RNOmni_1.0.1 [160] png_0.1-8 ewceData_1.6.0 BiocVersion_3.16.0 [163] bit_4.0.5 stringi_1.7.12 blob_1.2.3 [166] AnnotationHub_3.6.0 memoise_2.0.1 dplyr_1.1.0 [169] ape_5.6-2 ```