inst/doc/TxRegInfra.R

In TxRegInfra: Metadata management for multiomic specification of transcriptional regulatory networks

## ----setup,echo=FALSE,results="hide", eval=TRUE-------------------------------
suppressPackageStartupMessages({
library(TxRegInfra)
library(GenomicFiles)
library(TFutils)
})

## ----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

## ----lkpar, eval=TRUE---------------------------------------------------------
parent.env(con1)$orig

## ----getl, eval=TRUE----------------------------------------------------------
if (verifyHasMongoCmd()) {
  head(c1 <- listAllCollections(url=URL_txregInAWS(), db="txregnet"))
  }

## ----getl2, eval=TRUE---------------------------------------------------------
mongo(url=URL_txregInAWS(), db="txregnet", 
   collection="Adipose_Subcutaneous_allpairs_v7_eQTL")$find(limit=1)

## ----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")

## ----lkagggg, eval=TRUE-------------------------------------------------------
head(m1$aggregate(newagg))

## ----getcold, eval=TRUE-------------------------------------------------------
# cd = makeColData() # works when mongo does
cd = TxRegInfra::basicColData
head(cd,2)

## ----domor1, eval=TRUE--------------------------------------------------------
rme0 = RaggedMongoExpt(con1, colData=cd)
rme1 = rme0[, which(cd$type=="FP")]

## ----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")]

## ----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)

## ----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))

## ----lkeeee-------------------------------------------------------------------
names(mcols(s1))
head(s1[, c("gene_id", "variant_id", "maf", "pval_nominal")])

## ----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)

## ----lkgr---------------------------------------------------------------------
library(graph)
g1 = sbov_to_graphNEL(demo_eQTL_granges)
g1

## ----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

## ----doov2--------------------------------------------------------------------
fo2 = findOverlaps(demo_eQTL_granges, sbov_output_FP)
fo2 
eq_by_fp = split(demo_eQTL_granges[queryHits(fo2)],
   subjectHits(fo2))
eq_by_fp

## ----dotfs--------------------------------------------------------------------
library(TFutils)
data(demo_fimo_granges)
seqlevelsStyle(demo_eQTL_granges) = "UCSC"
lapply(demo_fimo_granges, lapply, function(x) 
   subsetByOverlaps(demo_eQTL_granges, x))