Add epi data; Add merging functions; Add haploR annot functions; Add summary plot functions

Author: bschilder

.DS_Store

@@ -0,0 +1,13 @@
+.*\.Rproj$
+^\.Rproj\.user$
+^README.Rmd
+^\.github$
+^doc$
+^Meta$
+^codecov\.yml$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$
+
+Dockerfile
+^LICENSE\.md$

.Rbuildignore

@@ -1,39 +1,42 @@
+# R project files
+*.Rproj
+.Rproj.user
+.Ruserdata 
 # History files
 .Rhistory
 .Rapp.history
-
 # Session Data files
 .RData
-
 # User-specific files
 .Ruserdata
-
+# .DS_Store
+# find . -name .DS_Store -print0 | xargs -0 git rm -f --ignore-unmatch
+./.DS_Store
+./**/.DS_Store
+./**/**/.DS_Store
+./**/**/**/.DS_Store
+./**/**/**/**/.DS_Store
+./**/**/**/**/**/.DS_Store
+./**/**/**/**/**/**/.DS_Store
 # Example code in package build process
 *-Ex.R
-
 # Output files from R CMD build
 /*.tar.gz
-
 # Output files from R CMD check
 /*.Rcheck/
-
 # RStudio files
 .Rproj.user/
-
 # produced vignettes
 vignettes/*.html
 vignettes/*.pdf
-
+vignettes/*.R
 # OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
 .httr-oauth
-
 # knitr and R markdown default cache directories
 *_cache/
 /cache/
-
 # Temporary files created by R markdown
 *.utf8.md
 *.knit.md
-
 # R Environment Variables
-.Renviron
+.Renviron

.gitignore

@@ -0,0 +1,68 @@
+Package: echoannot
+Type: Package
+Title: echoverse module: Annotate fine-mapping results
+Version: 0.99.1
+Authors@R: 
+    c(person(given = "Brian",
+           family = "Schilder",
+           role = c("cre"),
+           email = "brian_schilder@alumni.brown.edu",
+           comment = c(ORCID = "0000-0001-5949-2191")),
+    person(given = "Jack",
+           family = "Humphrey",
+           role = c("aut"),
+           email = "Jack.Humphrey@mssm.edu",
+           comment = c(ORCID = "0000-0002-6274-6620")),
+    person(given = "Towfique",
+           family = "Raj",
+           role = c("aut"),
+           email = "towfique.raj@mssm.edu",
+           comment = c(ORCID = "0000-0002-9355-5704"))
+           ) 
+Description: echoverse module: Annotate fine-mapping results.
+URL: https://github.com/RajLabMSSM/echoannot
+BugReports: https://github.com/RajLabMSSM/echoannot/issues
+Encoding: UTF-8
+LazyData: true
+Depends: R (>= 4.1)
+SystemRequirements: Python (>= 3.7.0)
+biocViews:
+Imports:  
+    echodata,
+    echotabix,
+    dplyr,
+    data.table,
+    utils, 
+    stats, 
+    tidyr,
+    parallel,
+    haploR,
+    ggplot2,
+    patchwork,
+    ggbio,
+    RColorBrewer,
+    scales,
+    GenomicRanges,
+    DescTools,
+    pheatmap,
+    grDevices,
+    rtracklayer,
+    S4Vectors,
+    GenomeInfoDb,
+    biomaRt
+Suggests:
+    markdown,
+    rmarkdown,
+    remotes,
+    knitr,
+    BiocStyle,
+    covr,
+    testthat (>= 3.0.0),
+    corrplot
+Remotes:
+    github::RajLabMSSM/echodata,
+    github::RajLabMSSM/echotabix
+RoxygenNote: 7.1.2
+VignetteBuilder: knitr
+License: GPL-3
+Config/testthat/edition: 3

DESCRIPTION

@@ -0,0 +1,70 @@
+# Generated by roxygen2: do not edit by hand
+
+export(CORCES_2020.get_HiChIP_FitHiChIP_overlap)
+export(CS_bin_plot)
+export(CS_counts_plot)
+export(merge_finemapping_results)
+export(peak_overlap_plot)
+export(plot_dataset_overlap)
+export(super_summary_plot)
+import(ggplot2)
+importFrom(DescTools,Divisors)
+importFrom(GenomeInfoDb,seqlevelsStyle)
+importFrom(GenomicRanges,GRanges)
+importFrom(GenomicRanges,elementMetadata)
+importFrom(GenomicRanges,end)
+importFrom(GenomicRanges,findOverlaps)
+importFrom(GenomicRanges,makeGRangesFromDataFrame)
+importFrom(GenomicRanges,mcols)
+importFrom(GenomicRanges,start)
+importFrom(IRanges,IRanges)
+importFrom(RColorBrewer,brewer.pal)
+importFrom(S4Vectors,queryHits)
+importFrom(S4Vectors,subjectHits)
+importFrom(biomaRt,getBM)
+importFrom(biomaRt,useMart)
+importFrom(data.table,as.data.table)
+importFrom(data.table,data.table)
+importFrom(data.table,fread)
+importFrom(data.table,fwrite)
+importFrom(data.table,melt.data.table)
+importFrom(data.table,merge.data.table)
+importFrom(data.table,rbindlist)
+importFrom(dplyr,"%>%")
+importFrom(dplyr,arrange)
+importFrom(dplyr,case_when)
+importFrom(dplyr,desc)
+importFrom(dplyr,funs)
+importFrom(dplyr,group_by)
+importFrom(dplyr,mutate)
+importFrom(dplyr,n_distinct)
+importFrom(dplyr,rename)
+importFrom(dplyr,select)
+importFrom(dplyr,slice)
+importFrom(dplyr,summarise)
+importFrom(dplyr,summarise_at)
+importFrom(dplyr,top_n)
+importFrom(dplyr,vars)
+importFrom(ggbio,autoplot)
+importFrom(ggbio,geom_arch)
+importFrom(ggbio,geom_rect)
+importFrom(ggbio,ggbio)
+importFrom(ggbio,ggsave)
+importFrom(ggbio,plotGrandLinear)
+importFrom(ggbio,scale_x_sequnit)
+importFrom(ggbio,theme_genome)
+importFrom(grDevices,dev.off)
+importFrom(grDevices,png)
+importFrom(haploR,queryRegulome)
+importFrom(parallel,mclapply)
+importFrom(patchwork,plot_layout)
+importFrom(patchwork,plot_spacer)
+importFrom(pheatmap,pheatmap)
+importFrom(rtracklayer,import.bw)
+importFrom(scales,alpha)
+importFrom(stats,as.formula)
+importFrom(stats,formula)
+importFrom(stats,median)
+importFrom(stats,setNames)
+importFrom(tidyr,separate)
+importFrom(utils,head)

NAMESPACE

@@ -0,0 +1,98 @@
+#' Get overlap between datatable of SNPs and scATAC peaks
+#'
+#' Can optionally add \code{Cicero} coaccessibility scores,
+#' which are also derived from scATAC-seq data.
+#'
+#' @family CORCES_2020
+#' @source \url{https://doi.org/10.1038/s41588-020-00721-x}
+CORCES_2020.get_ATAC_peak_overlap <- function(finemap_dat,
+                                              FDR_filter = NULL,
+                                              add_cicero = TRUE,
+                                              cell_type_specific = TRUE,
+                                              verbose = TRUE) {
+    FDR <- Cicero <- NULL
+
+    if (cell_type_specific) {
+        messager("CORCES_2020:: Extracting overlapping",
+            "cell-type-specific scATAC-seq peaks",
+            v = verbose
+        )
+        dat <- echoannot::CORCES_2020.scATACseq_celltype_peaks
+        Assay <- "scATAC"
+    } else {
+        messager("CORCES_2020:: Extracting overlapping",
+            "bulkATAC-seq peaks from brain tissue",
+            v = verbose
+        )
+        dat <- echoannot::CORCES_2020.bulkATACseq_peaks
+        Assay <- "bulkATAC"
+    }
+    gr.peaks_lifted <- echotabix::liftover(
+        sumstats_dt = dat,
+        convert_ref_genome = "hg19",
+        ref_genome = "hg38",
+        chrom_col = "hg38_Chromosome",
+        start_col = "hg38_Start",
+        end_col = "hg38_Stop",
+        as_granges = TRUE,
+        style = "NCBI",
+        verbose = FALSE
+    )
+    # Get overlap with PEAKS
+    gr.hits <- granges_overlap(
+        dat1 = finemap_dat,
+        chrom_col.1 = "CHR",
+        start_col.1 = "POS",
+        end_col.1 = "POS",
+        dat2 = gr.peaks_lifted
+    )
+    gr.hits$Assay <- Assay
+    if (!is.null(FDR_filter)) {
+        gr.hits <- subset(gr.hits, FDR < FDR_filter)
+    }
+
+    if (add_cicero & cell_type_specific) {
+        try({
+            # Pretty sure the Peak_IDs are shared between the
+            # sc-ATACseq data and cicero,
+            # because Cicero derives coaccess from sc-ATAC-seq data:
+            # http://www.cell.com/molecular-cell/retrieve/pii/S1097276518305471?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS1097276518305471%3Fshowall%3Dtrue
+            ## Also pretty sure that checking for cicero overlap
+            # only in the scATACseq gr.hits object is ok
+            # bc you can only test for coaccessibility if there's a
+            # peak to begin with.
+            cicero <- echoannot::CORCES_2020.cicero_coaccessibility
+            cicero_dict <- c(
+                stats::setNames(
+                    cicero$Coaccessibility,
+                    cicero$Peak_ID_Peak1
+                ),
+                stats::setNames(
+                    cicero$Coaccessibility,
+                    cicero$Peak_ID_Peak2
+                )
+            )
+            if (any(!is.na(cicero_dict[gr.hits$Peak_ID]))) {
+                gr.hits$Cicero <- cicero_dict[gr.hits$Peak_ID]
+                gr.cicero <- subset(gr.hits, !is.na(Cicero))
+                gr.cicero$Assay <- "Cicero"
+                messager(
+                    "+ CORCES_2020:: Cicero coaccessibility scores",
+                    "identified for",
+                    length(gr.cicero), "/", length(gr.hits), "peak hits.",
+                    v = verbose
+                )
+                gr.hits <- rbind_granges(
+                    gr1 = gr.hits,
+                    gr2 = gr.cicero
+                )
+            } else {
+                messager("+ CORCES_2020:: No Cicero hits found.", v = verbose)
+            }
+        })
+    }
+    if (cell_type_specific == FALSE) {
+        gr.hits$brain <- 1
+    }
+    return(gr.hits)
+}

R/CORCES_2020.get_ATAC_peak_overlap.R

@@ -0,0 +1,68 @@
+#' Get overlap between data table of SNPs and HiChIP_FitHiChIP
+#' coaccessibility anchors
+#'
+#' Anchors are the genomic regions that have evidence of being
+#' functionally connected to one another (coaccessible),
+#'  e.g. enhancer-promoter interactions.
+#'
+#' @param finemap_dat Fine-mapping results.
+#' @param verbose Print messages.
+#'
+#' @family CORCES_2020
+#' @source \url{https://doi.org/10.1038/s41588-020-00721-x}
+#' @export
+#'
+CORCES_2020.get_HiChIP_FitHiChIP_overlap <- function(finemap_dat,
+                                                     verbose = TRUE) {
+    loops <- echoannot::CORCES_2020.HiChIP_FitHiChIP_loop_calls
+    # Anchor 1
+    gr.anchor1 <- echotabix::liftover(
+        sumstats_dt = loops,
+        convert_ref_genome = "hg19",
+        chrom_col = "hg38_Chromosome_Anchor1",
+        start_col = "hg38_Start_Anchor1",
+        end_col = "hg38_Stop_Anchor1",
+        verbose = FALSE,
+        as_granges = TRUE,
+        style = "NCBI"
+    )
+    gr.anchor1_hits <- granges_overlap(
+        dat1 = finemap_dat,
+        chrom_col.1 = "CHR",
+        start_col.1 = "POS",
+        end_col.1 = "POS",
+        dat2 = gr.anchor1
+    )
+    gr.anchor1_hits$Anchor <- 1
+
+    # Anchor 2
+    gr.anchor2 <- echotabix::liftover(
+        sumstats_dt = loops,
+        convert_ref_genome = "hg19",
+        chrom_col = "hg38_Chromosome_Anchor2",
+        start_col = "hg38_Start_Anchor2",
+        end_col = "hg38_Stop_Anchor2",
+        verbose = FALSE,
+        as_granges = TRUE,
+        style = "NCBI"
+    )
+    gr.anchor2_hits <- granges_overlap(
+        dat1 = finemap_dat,
+        chrom_col.1 = "CHR",
+        start_col.1 = "POS",
+        end_col.1 = "POS",
+        dat2 = gr.anchor2
+    )
+    gr.anchor2_hits$Anchor <- 2
+    # Merge and report
+    gr.anchor <- rbind_granges(gr.anchor1_hits, gr.anchor2_hits)
+    gr.anchor$Assay <- "HiChIP_FitHiChIP"
+    # Have to make a pseudo cell-type col bc (i think)
+    # this analysis was done on bulk data
+    gr.anchor$brain <- 1
+    messager("+ CORCES_2020:: Found", length(gr.anchor),
+        "hits with HiChIP_FitHiChIP coaccessibility loop anchors.",
+        v = verbose
+    )
+    return(gr.anchor)
+}