rajlabmssm / echoplot Goto Github PK
View Code? Open in Web Editor NEWechoverse module: Locus plot creation for fine-mapping and colocalization studies.
Home Page: https://rajlabmssm.github.io/echoplot
echoverse module: Locus plot creation for fine-mapping and colocalization studies.
Home Page: https://rajlabmssm.github.io/echoplot
I think there is a bug somewhere trying to plot Nott's epigenetic data.
I guess somewhere I did not manage to find there is a "RUE" instead of a "TRUE"
finemap_loci(# GENERAL ARGUMENTS
topSNPs = topSNPs,
results_dir = fullRS_path,
loci = topSNPs$Locus,
dataset_name = "PiD_sURVIVAL_NOTT",
dataset_type = "GWAS",
force_new_subset = TRUE,
force_new_LD = FALSE,
force_new_finemap = FALSE,
remove_tmps = FALSE,
finemap_methods = c("SUSIE"),
# Munge full sumstats first
munged = FALSE,
colmap = columnsnames,
# SUMMARY STATS ARGUMENTS
fullSS_path = newSS_name_colmap,
fullSS_genome_build = "hg19",
query_by ="tabix",
bp_distance = 500000*2,
min_MAF = 0.001,
trim_gene_limits = FALSE,
case_control = TRUE,
# FINE-MAPPING ARGUMENTS
## General
n_causal = 5,
credset_thresh = .95,
consensus_thresh = 2,
# LD ARGUMENTS
LD_reference = "1KGphase3",#"UKB",
superpopulation = "EUR",
download_method = "axel",
LD_genome_build = "hg19",
leadSNP_LD_block = FALSE,
#### PLotting args ####
plot_types = c("fancy"),
show_plot = TRUE,
zoom = c("1x", "3x", "4x", "5x" ,"6x", "7x", "8x", "20x"),
#zoom = "1x",
tx_biotypes = NULL,
nott_epigenome = TRUE,
nott_show_placseq = TRUE,
nott_binwidth = 200,
nott_bigwig_dir = NULL,
#xgr_libnames =c("ENCODE_TFBS_ClusteredV3_CellTypes", "TFBS_Conserved", "Uniform_TFBS"),
#### General args ####
seed = 2022,
nThread = 20,
verbose = TRUE
)
────────────────────────────────────────────────────────────────────────────────
── Step 5 ▶▶▶ Plot 📈 ──────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────
+-------- Locus Plot: CSMD1 --------+
+ support_thresh = 2
+ Calculating mean Posterior Probability (mean.PP)...
+ 1 fine-mapping methods used.
+ 0 Credible Set SNPs identified.
+ 0 Consensus SNPs identified.
+ Filling NAs in CS cols with 0.
+ Filling NAs in PP cols with 0.
LD_matrix detected. Coloring SNPs by LD with lead SNP.
++ echoplot:: GWAS full window track
++ echoplot:: GWAS track
++ echoplot:: Merged fine-mapping track
Melting PP and CS from 2 fine-mapping methods.
+ echoplot:: Constructing SNP labels.
Adding SNP group labels to locus plot.
++ echoplot:: Adding Gene model track.
Converting dat to GRanges object.
max_transcripts= 1 .
2 transcripts from 2 genes returned.
Fetching data...OK
Parsing exons...OK
Defining introns...OK
Defining UTRs...OK
Defining CDS...OK
aggregating...
Done
Constructing graphics...
NOTT2019:: Creating epigenomic histograms plot
+ Inferring genomic limits for window: 1x
Constructing GRanges query using min/max ranges across one or more chromosomes.
Downloading data from UCSC.
Importing... [1] exvivo_H3K27ac_tbp
Importing... [2] microglia_H3K27ac
Importing... [3] neurons_H3K27ac
Importing... [4] oligodendrocytes_H3K27ac
Importing... [5] astrocytes_H3K27ac
Importing... [6] exvivo_atac_tbp
Importing... [7] microglia_atac
Importing... [8] neurons_atac
Importing... [9] oligodendrocytes_atac
Importing... [10] astrocytes_atac
Importing... [11] microglia_H3K4me3
Importing... [12] neurons_H3K4me3
Importing... [13] oligodendrocytes_H3K4me3
Importing... [14] astrocytes_H3K4me3
Saving bigwig query ==> /tmp/RtmpR6G9c0/CSMD1_Nott2019_bigwig.rds
Importing previously downloaded files: /home/rstudio/.cache/R/echoannot/NOTT2019_epigenomic_peaks.rds
++ NOTT2019:: 634,540 ranges retrieved.
dat is already a GRanges object.
361 query SNP(s) detected with reference overlap.
+ Calculating max histogram height
+ Converting label units to Mb.
NOTT2019:: Creating PLAC-seq interactome plot
++ NOTT2019:: Getting promoter cell-type-specific data.
++ NOTT2019:: Getting interactome data.
++ NOTT2019:: Getting regulatory regions data.
Importing Astrocyte enhancers ...
Importing Astrocyte promoters ...
Importing Neuronal enhancers ...
Importing Neuronal promoters ...
Importing Oligo enhancers ...
Importing Oligo promoters ...
Importing Microglia enhancers ...
Importing Microglia promoters ...
Converting dat to GRanges object.
++ NOTT2019:: Getting interaction anchors data.
Importing Microglia interactome ...
Importing Neuronal interactome ...
Importing Oligo interactome ...
Converting dat to GRanges object.
27 query SNP(s) detected with reference overlap.
Converting dat to GRanges object.
14 query SNP(s) detected with reference overlap.
Converting dat to GRanges object.
Preparing data for highlighting PLAC-seq interactions that overlap with SNP subset: Support>0
No target SNPs overlapped with PLAC-seq anchors.Locus CSMD1 complete in: 5.1 min
────────────────────────────────────────────────────────────────────────────────
── Step 6 ▶▶▶ Postprocess data 🎁 ──────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────
Returning results as nested list.
All loci done in: 9.08 min
$NLGN1
NULL
$KATNAL2
NULL
$CSMD1
NULL
$merged_dat
Null data.table (0 rows and 0 cols)
Warning messages:
1: In .Call(.make_vcf_geno, filename, fixed, names(geno), as.list(geno), :
converting NULL pointer to R NULL
2: In .local(x, ...) : variants with >1 ALT allele are set to NA
3: In .Call(.make_vcf_geno, filename, fixed, names(geno), as.list(geno), :
converting NULL pointer to R NULL
4: In .local(x, ...) : variants with >1 ALT allele are set to NA
5: In .Call(.make_vcf_geno, filename, fixed, names(geno), as.list(geno), :
converting NULL pointer to R NULL
6: In .local(x, ...) : variants with >1 ALT allele are set to NA
Error: object 'RUE' not found
Execution halted
> sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8
[8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] snpStats_1.48.0 Matrix_1.4-1 survival_3.3-1 SNPlocs.Hsapiens.dbSNP155.GRCh37_0.99.22
[5] susieR_0.12.28 SNPlocs.Hsapiens.dbSNP144.GRCh37_0.99.20 BSgenome_1.66.1 rtracklayer_1.58.0
[9] Biostrings_2.66.0 XVector_0.38.0 GenomicRanges_1.50.1 GenomeInfoDb_1.34.3
[13] IRanges_2.32.0 S4Vectors_0.36.0 BiocGenerics_0.44.0 data.table_1.14.2
[17] forcats_0.5.2 stringr_1.4.1 dplyr_1.0.10 purrr_0.3.4
[21] readr_2.1.2 tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6
[25] tidyverse_1.3.2 reticulate_1.26 coloc_5.1.0 echolocatoR_2.0.3
loaded via a namespace (and not attached):
[1] rappdirs_0.3.3 GGally_2.1.2 R.methodsS3_1.8.2 ragg_1.2.2 echoLD_0.99.8 bit64_4.0.5
[7] knitr_1.40 irlba_2.3.5 DelayedArray_0.24.0 R.utils_2.12.0 rpart_4.1.16 KEGGREST_1.38.0
[13] RCurl_1.98-1.8 AnnotationFilter_1.22.0 generics_0.1.3 GenomicFeatures_1.50.2 RSQLite_2.2.16 proxy_0.4-27
[19] bit_4.0.4 tzdb_0.3.0 xml2_1.3.3 lubridate_1.8.0 SummarizedExperiment_1.28.0 assertthat_0.2.1
[25] viridis_0.6.2 gargle_1.2.0 xfun_0.32 hms_1.1.2 fansi_1.0.3 restfulr_0.0.15
[31] progress_1.2.2 dbplyr_2.2.1 readxl_1.4.1 Rgraphviz_2.41.1 igraph_1.3.4 DBI_1.1.3
[37] htmlwidgets_1.5.4 reshape_0.8.9 downloadR_0.99.5 googledrive_2.0.0 ellipsis_0.3.2 ggnewscale_0.4.7
[43] backports_1.4.1 biomaRt_2.54.0 deldir_1.0-6 MatrixGenerics_1.10.0 MungeSumstats_1.7.1 vctrs_0.4.1
[49] Biobase_2.58.0 ensembldb_2.22.0 cachem_1.0.6 withr_2.5.0 checkmate_2.1.0 GenomicAlignments_1.34.0
[55] prettyunits_1.1.1 cluster_2.1.3 ape_5.6-2 dir.expiry_1.6.0 lazyeval_0.2.2 crayon_1.5.1
[61] basilisk.utils_1.10.0 crul_1.2.0 labeling_0.4.2 pkgconfig_2.0.3 nlme_3.1-159 ProtGenerics_1.30.0
[67] XGR_1.1.8 nnet_7.3-17 pals_1.7 rlang_1.0.5 lifecycle_1.0.1 filelock_1.0.2
[73] httpcode_0.3.0 BiocFileCache_2.6.0 modelr_0.1.9 echotabix_0.99.8 dichromat_2.0-0.1 cellranger_1.1.0
[79] matrixStats_0.62.0 graph_1.76.0 osfr_0.2.8 boot_1.3-28 reprex_2.0.2 base64enc_0.1-3
[85] googlesheets4_1.0.1 png_0.1-7 viridisLite_0.4.1 rjson_0.2.21 rootSolve_1.8.2.3 bitops_1.0-7
[91] R.oo_1.25.0 ggnetwork_0.5.10 blob_1.2.3 mixsqp_0.3-43 echoplot_0.99.6 dnet_1.1.7
[97] jpeg_0.1-9 echodata_0.99.16 scales_1.2.1 memoise_2.0.1 magrittr_2.0.3 plyr_1.8.7
[103] hexbin_1.28.2 zlibbioc_1.44.0 compiler_4.2.0 echoconda_0.99.8 BiocIO_1.8.0 RColorBrewer_1.1-3
[109] catalogueR_1.0.0 EnsDb.Hsapiens.v75_2.99.0 Rsamtools_2.14.0 cli_3.3.0 echoannot_0.99.10 patchwork_1.1.2
[115] htmlTable_2.4.1 Formula_1.2-4 MASS_7.3-58.1 tidyselect_1.1.2 stringi_1.7.8 textshaping_0.3.6
[121] yaml_2.3.5 supraHex_1.35.0 latticeExtra_0.6-30 ggrepel_0.9.1 grid_4.2.0 VariantAnnotation_1.44.0
[127] tools_4.2.0 lmom_2.9 parallel_4.2.0 rstudioapi_0.14 foreign_0.8-82 piggyback_0.1.3
[133] gridExtra_2.3 gld_2.6.5 farver_2.1.1 RcppZiggurat_0.1.6 digest_0.6.29 BiocManager_1.30.18
[139] Rcpp_1.0.9 broom_1.0.1 OrganismDbi_1.40.0 httr_1.4.4 AnnotationDbi_1.60.0 RCircos_1.2.2
[145] ggbio_1.46.0 biovizBase_1.46.0 colorspace_2.0-3 rvest_1.0.3 XML_3.99-0.10 fs_1.5.2
[151] splines_4.2.0 RBGL_1.74.0 expm_0.999-6 echofinemap_0.99.4 basilisk_1.9.12 Exact_3.1
[157] mapproj_1.2.8 systemfonts_1.0.4 jsonlite_1.8.0 Rfast_2.0.6 R6_2.5.1 Hmisc_4.7-1
[163] pillar_1.8.1 htmltools_0.5.3 glue_1.6.2 fastmap_1.1.0 DT_0.24 BiocParallel_1.32.1
[169] class_7.3-20 codetools_0.2-18 maps_3.4.0 mvtnorm_1.1-3 utf8_1.2.2 lattice_0.20-45
[175] curl_4.3.2 DescTools_0.99.46 zip_2.2.0 openxlsx_4.2.5 interp_1.1-3 googleAuthR_2.0.0
This was happening when the number of transcript retrieved was 0.
dat <- echodata::MEX3C[100:101,]
gene_track <- echoplot:: transcript_model_track(dat=dat)
Converting dat to GRanges object.
max_transcripts= 1 .
0 transcripts from 0 genes returned.
Fetching data...
Error: near ")": syntax error
Add in the manual all the different options we can pass to tx_biotypes
A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]
Describe the solution you'd like
Somewhere in manual show this
> unique(txdb_transcripts@elementMetadata@listData$tx_biotype)
[1] "processed_transcript" "transcribed_unprocessed_pseudogene" "unprocessed_pseudogene" "lincRNA"
[5] "miRNA" "protein_coding" "processed_pseudogene" "antisense"
[9] "snRNA" "pseudogene" "retained_intron" "nonsense_mediated_decay"
[13] "LRG_gene" "sense_intronic" "misc_RNA" "transcribed_processed_pseudogene"
[17] "snoRNA" "non_stop_decay" "rRNA" "unitary_pseudogene"
[21] "3prime_overlapping_ncrna" "polymorphic_pseudogene" "sense_overlapping" "IG_V_pseudogene"
[25] "TR_V_gene" "TR_V_pseudogene" "TR_D_gene" "TR_J_gene"
[29] "TR_C_gene" "TR_J_pseudogene" "translated_processed_pseudogene" "IG_C_gene"
[33] "IG_C_pseudogene" "IG_J_gene" "IG_J_pseudogene" "IG_D_gene"
[37] "IG_V_gene" "Mt_tRNA" "Mt_rRNA"
Corces_2020 single-cell epigenomic data is currently only used in multi-locus summary plots. Add as an option to PLOT.locus()
.
Example from vignette:
trk_plot.roadmap <- echoplot::plot_locus(dat = dat,
LD_matrix = LD_matrix,
LD_reference = LD_reference,
locus_dir = locus_dir,
roadmap = TRUE,
roadmap_query = "monocyte",
zoom = "5x")
+-------- Locus Plot: BST1 --------+
+ support_thresh = 2
+ Calculating mean Posterior Probability (mean.PP)...
+ 4 fine-mapping methods used.
+ 7 Credible Set SNPs identified.
+ 3 Consensus SNPs identified.
+ Filling NAs in CS cols with 0.
+ Filling NAs in PP cols with 0.
LD_matrix detected. Coloring SNPs by LD with lead SNP.
Filling r/r2 NAs with 0
++ PLOT:: GWAS full window track
++ PLOT:: GWAS track
++ PLOT:: Merged fine-mapping track
Melting PP and CS from 5 fine-mapping methods.
++ PLOT:: Adding Gene model track.
Converting dat to GRanges object.
max_transcripts= 1 .
15 transcripts from 15 genes returned.
Fetching data...OK
Parsing exons...OK
Defining introns...OK
Defining UTRs...OK
Defining CDS...OK
aggregating...
Done
Constructing graphics...
echoannot:: Plotting ROADMAP annotations.
Converting dat to GRanges object.
dat is already a GRanges object.
+ ROADMAP:: 2 annotation(s) identified that match `keyword_query`.
+ ROADMAP:: Querying subset from Roadmap API: E124 - 1/2
Downloading Roadmap Chromatin Marks: E124
Constructing GRanges query using min/max ranges within a single chromosome.
========= echotabix::query =========
Warning: The selected method (rsamtools) depends on Rhtslib, and the version you have installed (1.99.5) contains known bugs. Please install Rhtslib >=1.99.2 via Bioconductor >=3.16:
BiocManager::install(version='devel')Selecting valid alternative method: conda
query_dat is already a GRanges object. Returning directly.
Inferred format: 'table'
Selecting valid alternative method: seqminer
WARNING: 'seqminer' cannot query remote files. Switching query_method to 'rsamtools'.
Querying tabular tabix file using: Rsamtools.
Checking query chromosome style is correct.
[E::bgzf_read] Read block operation failed with error -1 after 0 of 8 bytes
Chromosome format: chr1
Retrieving data.
[E::bgzf_read] Read block operation failed with error -1 after 0 of 8 bytes
Converting query results to data.table.
Processing query: chr4:14737349-16737284
<simpleError in data.table::fread(text = paste(query, collapse = "\n"), header = FALSE, sep = sep, fill = TRUE): Input is empty or only contains BOM or terminal control characters>
Adding 'query' column to results.
Retrieved data with 0 rows
Saving query ==> /var/folders/zq/h7mtybc533b1qzkys_ttgpth0000gn/T//RtmpPBFcAL/filec8f848c6e4ee.gz
Warning: Input has no columns; creating an empty file at '/var/folders/zq/h7mtybc533b1qzkys_ttgpth0000gn/T//RtmpPBFcAL/filec8f848c6e4ee.gz' and exiting.Error in `[.data.table`(dat, , paste0("V", seq(1, 4))) :
column(s) not found: V1, V2, V3, V4
+ ROADMAP:: Querying subset from Roadmap API: E029 - 2/2
Downloading Roadmap Chromatin Marks: E029
Constructing GRanges query using min/max ranges within a single chromosome.
========= echotabix::query =========
Warning: The selected method (rsamtools) depends on Rhtslib, and the version you have installed (1.99.5) contains known bugs. Please install Rhtslib >=1.99.2 via Bioconductor >=3.16:
BiocManager::install(version='devel')Selecting valid alternative method: conda
query_dat is already a GRanges object. Returning directly.
Inferred format: 'table'
Selecting valid alternative method: seqminer
WARNING: 'seqminer' cannot query remote files. Switching query_method to 'rsamtools'.
Querying tabular tabix file using: Rsamtools.
Checking query chromosome style is correct.
[E::bgzf_read] Read block operation failed with error -1 after 0 of 8 bytes
Chromosome format: chr1
Retrieving data.
[E::bgzf_read] Read block operation failed with error -1 after 0 of 8 bytes
Converting query results to data.table.
Processing query: chr4:14737349-16737284
<simpleError in data.table::fread(text = paste(query, collapse = "\n"), header = FALSE, sep = sep, fill = TRUE): Input is empty or only contains BOM or terminal control characters>
Adding 'query' column to results.
Retrieved data with 0 rows
Saving query ==> /var/folders/zq/h7mtybc533b1qzkys_ttgpth0000gn/T//RtmpPBFcAL/filec8f879dfea95.gz
Warning: Input has no columns; creating an empty file at '/var/folders/zq/h7mtybc533b1qzkys_ttgpth0000gn/T//RtmpPBFcAL/filec8f879dfea95.gz' and exiting.Error in `[.data.table`(dat, , paste0("V", seq(1, 4))) :
column(s) not found: V1, V2, V3, V4
ROADMAP:: All downloads complete
12.3
Merging and processing ROADMAP annotations.
Generating ROADMAP track plot.
Warning: Ignoring unknown parameters: which, facetsWarning: no non-missing arguments to max; returning -Inf+ Adding vertical lines to highlight SNP groups...
Error in `combine_vars()`:
! Faceting variables must have at least one value
Backtrace:
1. echoplot::plot_locus(...)
2. echoplot:::add_multitrack_lines(...)
3. echoplot:::guess_genomic_units(gg = trk)
5. ggplot2:::ggplot_build.ggplot(gg)
6. layout$setup(data, plot$data, plot$plot_env)
7. ggplot2 (local) f(..., self = self)
8. self$facet$compute_layout(data, self$facet_params)
9. ggplot2 (local) f(...)
10. ggplot2::combine_vars(data, params$plot_env, rows, drop = params$drop)
Happens across all platforms:
Quitting from lines 69-76 (echoplot.Rmd)
Error: Error: processing vignette 'echoplot.Rmd' failed with diagnostics:
Can't combine `..1$Method` <factor<2f980>> and `..4$Method` <ordered<4e1b3>>.
Different error, but similar situation to this issue: RajLabMSSM/echolocatoR#129
This error happens when you try to add a plot with one color scale (e.g. SNP groups) to one that already has a color scale (e.g. LD)
Error in `ggplot2::geom_point()`:
! Problem while converting geom to grob.
ℹ Error occurred in the 8th layer.
Caused by error:
! colours encodes as numbers must be positive
ggnewscale
fixes this perfectly! Specifically, adding ggnewscale::new_scale(new_aes = "color") +
https://eliocamp.github.io/codigo-r/2018/09/multiple-color-and-fill-scales-with-ggplot2/
Hi Brian,
Sorry for the bad post but see an image of the finemapping plots.
I do not know if there is something going on in the selection of the top SNPs to highlight.
At the beginning I thought it could be related with the new names in the multifinemap df for the new version of FINEMAP but this is actually happening for susie as well.
The locus in the image is a case in which there is no consensus SNP, in case you think it might be related...
I cannot share the data for this but happy to have a chat.
If you try to reproduce with the roxygen example, please not the the toy data of locus plot is generated using a previous version of FINEMAP. Maybe worth getting the multifinemap dat with the latest version of FINEMAP, in case the issue is relared to this
This is happening on GHA because ensembldb
is currently failing checks on Bioc 3.17:
https://bioconductor.org/packages/devel/bioc/html/ensembldb.html
Only solutions are to provide the remote to the github, or wait until it is fixed:
https://github.com/jorainer/ensembldb
Is your feature request related to a problem? Please describe.
By exploring the filtering process in the txdb transcripts df derived from ensembldb, I thought of enhancement for this plot
I have noticed the transcripts plot gets overpopulated and tough to ready ( ie gene names overlapping)
Describe alternatives you've considered
Adding a flat to PLOT.locus( ..., prot_coding = T, ...), could allow the user to retrieve protein coding variants only.
The filtering could be performed while querying ensembldb object
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.