--- title: "TxRegInfra: support for TxRegQuery" author: "Vincent J. Carey, stvjc at channing.harvard.edu" date: "`r format(Sys.time(), '%B %d, %Y')`" vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{TxRegInfra -- classes and methods for TxRegQuery} %\VignetteEncoding{UTF-8} output: BiocStyle::html_document: highlight: pygments number_sections: yes theme: united toc: yes --- ```{r setup,echo=FALSE,results="hide", eval=TRUE} suppressPackageStartupMessages({ library(TxRegInfra) library(GenomicFiles) library(TFutils) }) ``` # Introduction TxRegQuery addresses exploration of transcriptional regulatory networks by integrating data on eQTL, digital genomic footprinting (DGF), DnaseI hypersensitivity binding data (DHS), and transcription factor binding site (TFBS) data. Owing to the volume of emerging tissue-specific data, special data modalities are used. # Managing bed file content with mongodb ## Querying the `txregnet` database We have a long-running server that will respond to queries. We focus on `r CRANpkg("mongolite")` as the interface. ### The connection ```{r lkmong, eval=TRUE} suppressPackageStartupMessages({ library(TxRegInfra) library(mongolite) library(Gviz) library(EnsDb.Hsapiens.v75) library(BiocParallel) register(SerialParam()) }) con1 = mongo(url=URL_txregInAWS(), db="txregnet") con1 ``` We will write methods that work with the 'fields' of this object. There is not much explicit reflectance in the mongolite API. The following is improvised and may be fragile: ```{r lkpar, eval=TRUE} parent.env(con1)$orig ``` ### Queries and aggregation If the `mongo` utility is available as a system command, we can get a list of collections in the database as follows. ```{r getl, eval=TRUE} if (verifyHasMongoCmd()) { head(c1 <- listAllCollections(url=URL_txregInAWS(), db="txregnet")) } ``` Otherwise, as long as `r CRANpkg("mongolite")` is installed, as long as we know the collection names of interest, we can use them as noted throughout this vignette. We can get a record from a given collection: ```{r getl2, eval=TRUE} mongo(url=URL_txregInAWS(), db="txregnet", collection="Adipose_Subcutaneous_allpairs_v7_eQTL")$find(limit=1) ``` Queries can be composed using JSON. We have a tool to generate queries that employ the mongodb aggregation method. Here we demonstrate this by computing, for each chromosome, the count and minimum values of the footprint statistic on CD14 cells. ```{r doagg, eval=TRUE} m1 = mongo(url = URL_txregInAWS(), db = "txregnet", collection="CD14_DS17215_hg19_FP") newagg = makeAggregator( by="chr", vbl="stat", op="$min", opname="min") ``` The JSON layout of this aggregating query is ``` [ { "$group": { "_id": ["$chr"], "count": { "$sum": [1] }, "min": { "$min": ["$stat"] } } } ] ``` Invocation returns a data frame: ```{r lkagggg, eval=TRUE} head(m1$aggregate(newagg)) ``` # An integrative container We need to bind the metadata and information about the mongodb. ## Sample metadata The following turns a very ad hoc filtering of the collection names into a DataFrame. ```{r getcold, eval=TRUE} # cd = makeColData() # works when mongo does cd = TxRegInfra::basicColData head(cd,2) ``` ## Extended RaggedExperiment ```{r domor1, eval=TRUE} rme0 = RaggedMongoExpt(con1, colData=cd) rme1 = rme0[, which(cd$type=="FP")] ``` A key method in development is subsetting the archive by genomic coordinates. ```{r lksb, cache=TRUE, eval=TRUE} si = GenomeInfoDb::Seqinfo(genome="hg19")["chr17"] # to fix query genome myg = GRanges("chr17", IRanges(38.07e6,38.09e6), seqinfo=si) s1 = sbov(rme1, myg, simplify=FALSE) s1 dim(sa <- sparseAssay(s1, 3)) # compact gives segfault sa[953:956,c("fLung_DS14724_hg19_FP", "fMuscle_arm_DS17765_hg19_FP")] ``` # Visualizing coincidence ```{r mym, eval=TRUE} ormm = txmodels("ORMDL3", plot=FALSE, name="ORMDL3") sar = strsplit(rownames(sa), ":|-") an = as.numeric gr = GRanges(seqnames(ormm)[1], IRanges(an(sapply(sar,"[", 2)), an(sapply(sar,"[", 3)))) gr1 = gr gr1$score = 1-sa[,1] gr2 = gr gr2$score = 1-sa[,2] sc1 = DataTrack(gr1, name="Lung FP") sc2 = DataTrack(gr2, name="Musc/Arm FP") plotTracks(list(GenomeAxisTrack(), sc1, sc2, ormm), showId=TRUE) ``` # Higher-level work with `sbov` ## Building annotated GRanges for a selected target interval We begin with three 'single-concept' assays with relevance to lung genomics. The v7 GTEx lung eQTL data, an encode DnaseI narrowPeak report on lung fibroblasts, and a digital genomic footprint report for fetal lung. ```{r lksbovs} lname_eqtl = "Lung_allpairs_v7_eQTL" lname_dhs = "ENCFF001SSA_hg19_HS" # see dnmeta, fibroblast of lung lname_fp = "fLung_DS14724_hg19_FP" si17 = GenomeInfoDb::Seqinfo(genome="hg19")["chr17"] si17n = si17 GenomeInfoDb::seqlevelsStyle(si17n) = "NCBI" s1 = sbov(rme0[,lname_eqtl], GRanges("17", IRanges(38.06e6, 38.15e6), seqinfo=si17n)) s2 = sbov(rme0[,lname_dhs], GRanges("chr17", IRanges(38.06e6, 38.15e6), seqinfo=si17)) s3 = sbov(rme0[,lname_fp], GRanges("chr17", IRanges(38.06e6, 38.15e6), seqinfo=si17)) ``` Now we have annotated GRanges for each assay. The eQTL data in part are: ```{r lkeeee} names(mcols(s1)) head(s1[, c("gene_id", "variant_id", "maf", "pval_nominal")]) ``` The names of genes and variants used here are cumbersome -- symbols and rsids are preferable. ```{r doadd} addsyms = function(x, EnsDb=EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75) { ensids = gsub("\\..*", "", x$gene_id) # remove post period gns = genes(EnsDb) x$symbol = gns[ensids]$symbol x } s1 = addsyms(s1) ``` Note that it is possible to retrieve rsids for the SNPs by address. But this is a slow operation involving a huge SNPlocs package that we do not want to work with directly for this vignette. ``` > snpsByOverlaps(SNPlocs.Hsapiens.dbSNP144.GRCh37, s1b) UnstitchedGPos object with 265 positions and 2 metadata columns: seqnames pos strand | RefSNP_id alleles_as_ambig | [1] 17 38061054 * | rs36049276 R [2] 17 38061439 * | rs4795399 Y [3] 17 38062196 * | rs2305480 R [4] 17 38062217 * | rs2305479 Y [5] 17 38062503 * | rs35104165 Y ... ... ... ... . ... ... [261] 17 38149258 * | rs58212353 K [262] 17 38149350 * | rs8073254 V [263] 17 38149411 * | rs34648856 R [264] 17 38149724 * | rs3785549 Y [265] 17 38149727 * | rs3785550 H ------- seqinfo: 25 sequences (1 circular) from GRCh37.p13 genome ``` ## A bipartite graph for eQTL-gene relationships The object `s1` computed above is available as `demo_eQTL_granges`. We convert it to a graph via ```{r lkgr} library(graph) g1 = sbov_to_graphNEL(demo_eQTL_granges) g1 ``` Nodes are SNPs and genes, edges are present when the resource (in this case the GTEx lung study) declares an association (in this case, an FDR for SNP-gene association not exceeding 0.10.) The `r Biocpkg("graph")` library includes functions for creation of incidence matrices from graphs, and vice versa. ## Connecting eQTL-SNPs via DHS and DGF Given the GRanges representations for `sbov` results, we can use overlap computations to conveniently identify relationships between eQTL SNPs, genes, and hypersensitivity or footprint regions. We use `sbov_output_HS` as a persistent instance of `s2` computed above. ```{r doov} seqlevelsStyle(demo_eQTL_granges) = "UCSC" fo1 = findOverlaps(demo_eQTL_granges, sbov_output_HS) fo1 eq_by_hs = split(demo_eQTL_granges[queryHits(fo1)], subjectHits(fo1)) eq_by_hs ``` This shows that there are two DHS sites that overlap with SNPs showing eQTL associations with various genes. For the footprint data, we have: ```{r doov2} fo2 = findOverlaps(demo_eQTL_granges, sbov_output_FP) fo2 eq_by_fp = split(demo_eQTL_granges[queryHits(fo2)], subjectHits(fo2)) eq_by_fp ``` ## Relationships to FIMO-based TFBS We have a small number of cloud-resident FIMO search results through the `r Biocpkg("TFutils")` package. ```{r dotfs} library(TFutils) data(demo_fimo_granges) seqlevelsStyle(demo_eQTL_granges) = "UCSC" lapply(demo_fimo_granges, lapply, function(x) subsetByOverlaps(demo_eQTL_granges, x)) ```