### R code from vignette source 'vignettes/rTRM/inst/doc/rTRM_Introduction.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: rTRM_Introduction.Rnw:21-30
###################################################
library(knitr)
library(rTRM)
library(org.Mm.eg.db)

opts_chunk$set(fig.align = 'center', fig.width = 3, fig.height = 3, fig.show = 'hold', tidy = TRUE, comment = "", highlight = FALSE, prompt = TRUE)
options(replace.assign=TRUE,width = 80)
knit_hooks$set(no.mar = function(before, options, envir) {
    if (before) par(mar = rep(0,4))  # no margins.
})


###################################################
### code chunk number 2: rTRM_Introduction.Rnw:53-76
###################################################
# load the rTRM package
library(rTRM)

# load network example.
load(system.file(package = "rTRM", "extra/example.rda"))

# plot network
plotGraph(g, vertex.cex = 5)

# define target and query nodes:
target = "N6"
query = c("N7", "N12", "N28")

# find TRM:
s = findTRM(g, target = target, query = query, method = "nsa", max.bridge = 1)

# annotate nodes:
V(s)$color = "white"
V(s)[ name %in% query]$color = "steelblue2"
V(s)[ name %in% target]$color = "steelblue4"

# plot:
plotGraph(s, mar = 5, vertex.cex = 5)


###################################################
### code chunk number 3: rTRM_Introduction.Rnw:88-90
###################################################
pwm = getMatrices()
head(pwm, 1)


###################################################
### code chunk number 4: rTRM_Introduction.Rnw:94-96
###################################################
ann = getAnnotations()
head(ann)


###################################################
### code chunk number 5: rTRM_Introduction.Rnw:100-102
###################################################
map = getMaps()
head(map)


###################################################
### code chunk number 6: rTRM_Introduction.Rnw:106-108
###################################################
o = getOrthologs(organism = "mouse")
head(o)


###################################################
### code chunk number 7: rTRM_Introduction.Rnw:113-115
###################################################
getOrthologFromMatrix("MA0009.1", organism="human")
getOrthologFromMatrix("MA0009.1", organism="mouse")


###################################################
### code chunk number 8: rTRM_Introduction.Rnw:123-127
###################################################
# check statistics about the network.
biogrid_mm()
# load mouse PPI network:
data(biogrid_mm)


###################################################
### code chunk number 9: rTRM_Introduction.Rnw:132-139 (eval = FALSE)
###################################################
## # obtain dataset (currently need to state the version number explicitly, see
## # http://www.thebiogrid.org to identify the latest version:
## db = getBiogridData("3.2.96")
## 
## # process PPI data for different organisms (currently supported human and mouse):
## biogrid_hs = processBiogrid(db, org = "human")
## biogrid_mm = processBiogrid(db, org = "mouse")


###################################################
### code chunk number 10: rTRM_Introduction.Rnw:148-168 (eval = FALSE)
###################################################
## library(PSICQUIC)
## psicquic <- PSICQUIC()
## providers(psicquic)
## 
## # obtain BioGrid human PPIs (as data.frame):
## tbl = interactions(psicquic, species="9606",provider="BioGrid")
## 
## # the target and source node information needs to be polished (i.e. must be Entrez gene id only)
## biogrid_hs = data.frame(source=tbl$A,target=tbl$B)
## biogrid_hs$source = sub(".*locuslink:(.*)\\|BIOGRID:.*","\\1", biogrid_hs$source)
## biogrid_hs$target = sub(".*locuslink:(.*)\\|BIOGRID:.*","\\1", biogrid_hs$target)
## 
## # create graph.
## library(igraph)
## biogrid_hs=graph.data.frame(biogrid_hs,directed=FALSE)
## biogrid_hs=simplify(biogrid_hs)
## 
## # annotate with symbols.
## library(org.Hs.eg.db)
## V(biogrid_hs)$label=select(org.Hs.eg.db,keys=V(biogrid_hs)$name,columns=c("SYMBOL"))$SYMBOL


###################################################
### code chunk number 11: rTRM_Introduction.Rnw:177-184
###################################################
#library(org.Mm.eg.db)

# read motif enrichment results.
motif_file = system.file("extra/sox2_motif_list.rda", package = "rTRM")
load(motif_file)
length(motif_list)
head(motif_list)


###################################################
### code chunk number 12: rTRM_Introduction.Rnw:189-194
###################################################
# get the corresponding gene.
tfs_list = getOrthologFromMatrix(motif_list, organism = "mouse")
tfs_list = unique(unlist(tfs_list, use.names = FALSE))
length(tfs_list)
head(tfs_list)


###################################################
### code chunk number 13: rTRM_Introduction.Rnw:199-208
###################################################
# load expression data.
eg_esc_file = system.file("extra/ESC-expressed.txt", package = "rTRM")
eg_esc = scan(eg_esc_file, what = "")
length(eg_esc)
head(eg_esc)

tfs_list_esc = tfs_list[tfs_list %in% eg_esc]
length(tfs_list_esc)
head(tfs_list_esc)


###################################################
### code chunk number 14: rTRM_Introduction.Rnw:212-235
###################################################
# load and process PPI data.
biogrid_mm()
data(biogrid_mm)
ppi = biogrid_mm
vcount(ppi)
ecount(ppi)

# remove outliers.
f = c("Ubc", "Sumo1", "Sumo2", "Sumo3")
f = unlist(mget(f, revmap(org.Mm.egSYMBOL)))
ppi = removeVertices(ppi, f)
vcount(ppi)
ecount(ppi)

# filter by expression.
ppi_esc = induced.subgraph(ppi, V(ppi)[ name %in% eg_esc ])
vcount(ppi_esc)
ecount(ppi_esc)

# ensure a single component.
ppi_esc = getLargestComp(ppi_esc)
vcount(ppi_esc)
ecount(ppi_esc)


###################################################
### code chunk number 15: rTRM_Introduction.Rnw:240-248
###################################################
# define target.
target = unlist(mget("Sox2", org.Mm.egSYMBOL2EG))
target

# find TRM.
s = findTRM(ppi_esc, target, tfs_list_esc, method = "nsa", max.bridge = 1)
vcount(s)
ecount(s)


###################################################
### code chunk number 16: rTRM_Introduction.Rnw:252-259
###################################################
# generate layout (order by cluster, then label)
cl = getConcentricList(s, target, tfs_list_esc)
l = layout.concentric(s, cl, order = "label")

# plot TRM.
plotTRM(s, layout = l, vertex.cex = 15, label.cex = .8)
plotTRMlegend(s, title = "ESC Sox2 TRM", cex = .8)


###################################################
### code chunk number 17: rTRM_Introduction.Rnw:269-278
###################################################
library(rTRM)
library(BSgenome.Mmusculus.UCSC.mm8) # Sox2 peaks found against mm8
library(PWMEnrich)
registerCoresPWMEnrich(1) # register number of cores for parallelization in PWMEnrich
library(MotifDb)

# select mouse PWMs:
sel.mm = values(MotifDb)$organism %in% c("Mmusculus")
pwm.mm = MotifDb[sel.mm]


###################################################
### code chunk number 18: rTRM_Introduction.Rnw:283-287
###################################################
# generate logn background model of PWMs:
p = as.list(pwm.mm)
p = lapply(p, function(x) round(x * 100))
p = lapply(p, function(x) t(apply(x, 1, as.integer)))


###################################################
### code chunk number 19: rTRM_Introduction.Rnw:292-293 (eval = FALSE)
###################################################
## pwm_logn = makeBackground(p, Mmusculus, type = "logn")


###################################################
### code chunk number 20: rTRM_Introduction.Rnw:295-296
###################################################
load(system.file("extra/pwm_mm_logn.rda", package = "rTRM"))


###################################################
### code chunk number 21: rTRM_Introduction.Rnw:301-309
###################################################
sox2_bed = read.table(system.file("extra/ESC_Sox2_peaks.txt", package = "rTRM"))

colnames(sox2_bed) = c("chr", "start", "end")

sox2_seq = getSequencesFromGenome(sox2_bed, Mmusculus, append.id="Sox2")
# PWMEnrich throws an error if the sequences are shorter than the motifs so we filter those sequences.
min.width = max(sapply(p, ncol))
sox2_seq_filter = sox2_seq[width(sox2_seq) >= min.width]


###################################################
### code chunk number 22: rTRM_Introduction.Rnw:314-316 (eval = FALSE)
###################################################
## # find enrichment:
## sox2_enr = motifEnrichment(sox2_seq_filter, pwms=pwm_logn, group.only=TRUE)


###################################################
### code chunk number 23: rTRM_Introduction.Rnw:319-320
###################################################
load(system.file("extra/sox2_enr.rda", package = "rTRM"))


###################################################
### code chunk number 24: rTRM_Introduction.Rnw:325-358
###################################################
res = motifRankingForGroup(sox2_enr)
res.gene = unique(values(MotifDb[names(res[res < 0.05])])$geneId)
res.gene = res.gene[!is.na(res.gene)]

data(biogrid_mm)
ppi = biogrid_mm
vcount(ppi)
ecount(ppi)

f = c("Ubc", "Sumo1", "Sumo2", "Sumo3")
f = unlist(mget(f, revmap(org.Mm.egSYMBOL)))
ppi = removeVertices(ppi, f)
vcount(ppi)
ecount(ppi)

# filter by expression.
eg_esc = scan(system.file("extra/ESC-expressed.txt", package = "rTRM"), what = "")
ppi_esc = induced.subgraph(ppi, V(ppi)[ name %in% eg_esc ])
vcount(ppi_esc)
ecount(ppi_esc)

# ensure a single component.
ppi_esc = getLargestComp(ppi_esc)
vcount(ppi_esc)
ecount(ppi_esc)

sox2.gene = unlist(mget("Sox2", org.Mm.egSYMBOL2EG))
sox2_trm = findTRM(ppi_esc, target=sox2.gene, query = res.gene)

cl = getConcentricList(sox2_trm, t=sox2.gene,e=res.gene)
l = layout.concentric(sox2_trm, concentric=cl, order="label")
plotTRM(sox2_trm, layout = l, vertex.cex = 15, label.cex = .8)
plotTRMlegend(sox2_trm, title = "ESC Sox2 TRM", cex = .8)


###################################################
### code chunk number 25: rTRM_Introduction.Rnw:365-368
###################################################
plotTRM(sox2_trm, layout = l, vertex.cex = 15, label.cex = .7)
l=layout.arc(sox2_trm,target=sox2.gene,query=res.gene)
plotTRM(sox2_trm, layout=l,vertex.cex=15,label.cex=.7)


###################################################
### code chunk number 26: rTRM_Introduction.Rnw:375-376
###################################################
citation(package="rTRM")


###################################################
### code chunk number 27: rTRM_Introduction.Rnw:381-382
###################################################
sessionInfo()