## ----setup, include = FALSE---------------------------------------------------
knitr::knit_hooks$set(optipng = knitr::hook_optipng)

## ----load-libs, message = FALSE,  warning = FALSE, results = FALSE------------
library(BulkSignalR)
library(igraph)
library(STexampleData)

## ----parallel, message=FALSE, warning=FALSE-----------------------------------
library(doParallel)
n.proc <- 1
cl <- makeCluster(n.proc)
registerDoParallel(cl)

# To add at the end of your script:
# stopCluster(cl)

## ----loading, eval=TRUE-------------------------------------------------------
data(sdc)
head(sdc)

## ----BSRDataModel, eval=TRUE,cache=FALSE--------------------------------------
bsrdm <- BSRDataModel(counts = sdc)
bsrdm

## ----learnParameters,eval=TRUE,warning=FALSE,fig.dim = c(7,3)-----------------
bsrdm <- learnParameters(bsrdm, quick=TRUE)
bsrdm

## ----BSRInference,eval=TRUE---------------------------------------------------

# We use a subset of the reference to speed up
# inference in the context of the vignette.
subset <- c("REACTOME_BASIGIN_INTERACTIONS",
"REACTOME_SYNDECAN_INTERACTIONS",
"REACTOME_ECM_PROTEOGLYCANS",
"REACTOME_CELL_JUNCTION_ORGANIZATION")

reactSubset <- BulkSignalR:::.SignalR$BulkSignalR_Reactome[
BulkSignalR:::.SignalR$BulkSignalR_Reactome$`Reactome name` %in% subset,]

resetPathways(dataframe = reactSubset,
              resourceName = "Reactome")

bsrinf <- BSRInference(bsrdm,
    min.cor = 0.3,
    reference="REACTOME")

LRinter.dataframe <- LRinter(bsrinf)

head(LRinter.dataframe[
    order(LRinter.dataframe$qval),
    c("L", "R", "LR.corr", "pw.name", "qval")])


## ----BSRInferenceTfile, eval=FALSE--------------------------------------------
# write.table(LRinter.dataframe[order(LRinter.dataframe$qval), ],
#     "./sdc_LR.tsv",
#     row.names = FALSE,
#     sep = "\t",
#     quote = FALSE
# )

## ----ReduceToPathway, eval=TRUE-----------------------------------------------
bsrinf.redP <- reduceToPathway(bsrinf)

## ----ReduceToBestPathway, eval=TRUE-------------------------------------------
bsrinf.redBP <- reduceToBestPathway(bsrinf)

## ----ReduceToLigand,eval=TRUE-------------------------------------------------
bsrinf.L <- reduceToLigand(bsrinf)
bsrinf.R <- reduceToReceptor(bsrinf)

## ----doubleReduction, eval=TRUE-----------------------------------------------
bsrinf.redP <- reduceToPathway(bsrinf)
bsrinf.redPBP <- reduceToBestPathway(bsrinf.redP)

## ----scoringLR,eval=TRUE,fig.dim = c(5,4)-------------------------------------
bsrsig.redBP <- BSRSignature(bsrinf.redBP, qval.thres=0.001)

scoresLR <- scoreLRGeneSignatures(bsrdm, bsrsig.redBP,
    name.by.pathway=FALSE
)

simpleHeatmap(scoresLR[1:20, ],
    hcl.palette="Cividis",
    pointsize=8)


## ----scoringPathway,eval=TRUE,fig.dim = c(7,3)--------------------------------

bsrsig.redPBP <- BSRSignature(bsrinf.redPBP, qval.thres=0.01)

scoresPathway <- scoreLRGeneSignatures(bsrdm, bsrsig.redPBP,
    name.by.pathway=TRUE
)

simpleHeatmap(scoresPathway,
    hcl.palette="Blue-Red 2",
    pointsize=8)

## ----heatmapMulti,eval=TRUE,fig.dim = c(7,10)---------------------------------
pathway1 <- pathways(bsrsig.redPBP)[1]
signatureHeatmaps(
    pathway = pathway1,
    bsrdm = bsrdm,
    heights = c(3,2,15),
    bsrsig = bsrsig.redPBP
)

## ----AlluvialPlot,eval=TRUE,fig.dim = c(8,3)----------------------------------
alluvialPlot(bsrinf,
    keywords = c("LAMC1"),
    type = "L",
    qval.thres = 0.01
)

## ----BubblePlot,eval=TRUE,fig.dim = c(8,4)------------------------------------

pathways <- LRinter(bsrinf)[1,c("pw.name")]
bubblePlotPathwaysLR(bsrinf,
    pathways = pathways,
    qval.thres = 0.001,
    color = "red",
    pointsize = 8
)

## ----Chordiagram,eval=TRUE,fig.dim = c(6,4.5)---------------------------------
chordDiagramLR(bsrinf,
    pw.id.filter = "R-HSA-210991",
    limit = 20,
    ligand="FYN", 
    receptor="SPN"
)

## ----network1, eval=TRUE------------------------------------------------------
# Generate a ligand-receptor network and export it in .graphML
# for Cytoscape or similar tools
gLR <- getLRNetwork(bsrinf.redBP, qval.thres = 1e-3)

# save to file
# write.graph(gLR,file="SDC-LR-network.graphml",format="graphml")

# As an alternative to Cytoscape, you can play with igraph package functions.
plot(gLR,
    edge.arrow.size = 0.1,
    vertex.label.color = "black",
    vertex.label.family = "Helvetica",
    vertex.label.cex = 0.1
)

## ----network2,  eval=TRUE-----------------------------------------------------
# You can apply other functions.


# Community detection
u.gLR <- as_undirected(gLR) # most algorithms work for undirected graphs only
comm <- cluster_edge_betweenness(u.gLR)
# plot(comm,u.gLR,
#     vertex.label.color="black",
#     vertex.label.family="Helvetica",
#     vertex.label.cex=0.1)

# Cohesive blocks
cb <- cohesive_blocks(u.gLR)
plot(cb, u.gLR,
    vertex.label.color = "black",
    vertex.label.family = "Helvetica",
    vertex.label.cex = 0.1,
    edge.color = "black"
)

## ----network3, eval=FALSE,warning=FALSE---------------------------------------
# # For the next steps, we just share the code below, but graph generation
# # functions are commented to lighten the vignette.
# 
# # Generate a ligand-receptor network complemented with intra-cellular,
# # receptor downstream pathways [computations are a bit longer here]
# #
# # You can save to a file for cystoscape or plot with igraph.
# 
# gLRintra <- getLRIntracellNetwork(bsrinf.redBP, qval.thres = 1e-3)
# 
# lay <- layout_with_kk(gLRintra)
# # plot(gLRintra,
# #     layout=lay,
# #     edge.arrow.size=0.1,
# #     vertex.label.color="black",
# #     vertex.label.family="Helvetica",
# #     vertex.label.cex=0.1)
# 
# # Reduce complexity by focusing on strongly targeted pathways
# pairs <- LRinter(bsrinf.redBP)
# top <- unique(pairs[pairs$pval <  1e-3, c("pw.id", "pw.name")])
# top
# gLRintra.res <- getLRIntracellNetwork(bsrinf.redBP,
#     qval.thres = 0.01,
#     restrict.pw = top$pw.id
# )
# lay <- layout_with_fr(gLRintra.res)
# 
# # plot(gLRintra.res,
# #     layout=lay,
# #     edge.arrow.size=0.1,
# #     vertex.label.color="black",
# #     vertex.label.family="Helvetica",
# #     vertex.label.cex=0.4)

## ----mouse,eval=TRUE,warning=FALSE--------------------------------------------
data(bodyMap.mouse)

ortholog.dict <- findOrthoGenes(
    from_organism = "mmusculus",
    from_values = rownames(bodyMap.mouse)
)

matrix.expression.human <- convertToHuman(
    counts = bodyMap.mouse,
    dictionary = ortholog.dict
)

bsrdm <- BSRDataModel(
    counts = matrix.expression.human,
    species = "mmusculus",
    conversion.dict = ortholog.dict
)

bsrdm <- learnParameters(bsrdm,quick=TRUE)

bsrinf <- BSRInference(bsrdm,reference="REACTOME")

bsrinf <- resetToInitialOrganism(bsrinf, conversion.dict = ortholog.dict)

# For example, if you want to explore L-R interactions
# you can proceed as shown above for a human dataset.

bsrinf.redBP <- reduceToBestPathway(bsrinf)
bsrsig.redBP <- BSRSignature(bsrinf.redBP, qval.thres=0.001)
scoresLR <- scoreLRGeneSignatures(bsrdm, bsrsig.redBP, name.by.pathway=FALSE)
head(LRinter(bsrinf.redBP))
#simpleHeatmap(scoresLR, column.names=TRUE,
#              pointsize=8)

## ----spatial1,eval=TRUE,message=FALSE-----------------------------------------
# load data =================================================

# We re-initialize the environment variable of Reactome
# to have all the pathways (we reduced it to a subset above to fasten
# example computations)
reactSubset <- getResource(resourceName = "Reactome",
cache = TRUE)

resetPathways(dataframe = reactSubset,
resourceName = "Reactome")

# Few steps of pre-process to subset a spatialExperiment object
# from STexampleData package ==================================

spe <- Visium_humanDLPFC()
set.seed(123)

speSubset <- spe[, colData(spe)$ground_truth%in%c("Layer1","Layer2")]

idx <- sample(ncol(speSubset), 10)
speSubset <- speSubset[, idx]

my.image.as.raster <- SpatialExperiment::imgRaster(speSubset, 
    sample_id = imgData(spe)$sample_id[1], image_id = "lowres")

colData(speSubset)$idSpatial <- paste(colData(speSubset)[[4]],
                colData(speSubset)[[5]],sep = "x")


annotation <- colData(speSubset)


## ----spatial2,eval=TRUE,warning=FALSE-----------------------------------------
# prepare data =================================================

bsrdm <- BSRDataModel(speSubset,
    min.count = 1,
    prop = 0.01,
    method = "TC",
    symbol.col = 2,
    x.col = 4,
    y.col = 5, 
    barcodeID.col = 1)

bsrdm <- learnParameters(bsrdm,
    quick = TRUE,
    min.positive = 2,
    verbose = TRUE)

bsrinf <- BSRInference(bsrdm, min.cor = -1,reference="REACTOME")

# spatial analysis ============================================

bsrinf.red <- reduceToBestPathway(bsrinf)
pairs.red <- LRinter(bsrinf.red)

thres <- 0.01
min.corr <- 0.01
pairs.red <- pairs.red[pairs.red$qval < thres & pairs.red$LR.corr > min.corr,]

head(pairs.red[
    order(pairs.red$qval),
    c("L", "R", "LR.corr", "pw.name", "qval")])

s.red  <- BSRSignature(bsrinf.red, qval.thres=thres)
scores.red <- scoreLRGeneSignatures(bsrdm,s.red)

head(scores.red)


## ----spatialPlot3,eval=TRUE,fig.dim = c(6,4.5)--------------------------------
# Visualization ============================================

# plot one specific interaction

# we have to follow the syntax with {} 
# to be compatible with the reduction functions
inter <- "{SLIT2} / {GPC1}"

# with raw tissue reference
spatialPlot(scores.red[inter, ], annotation, inter,
    ref.plot = TRUE, ref.plot.only = FALSE,
    image.raster = NULL, dot.size = 1,
    label.col = "ground_truth"
)

# or with synthetic image reference
spatialPlot(scores.red[inter, ], annotation, inter,
    ref.plot = TRUE, ref.plot.only = FALSE,
    image.raster = my.image.as.raster, dot.size = 1,
    label.col = "ground_truth"

)

## ----spatialPlot4,eval=TRUE,fig.dim = c(6,4.5)--------------------------------

separatedLRPlot(scores.red, "SLIT2", "GPC1", 
    ncounts(bsrdm), 
    annotation,
    label.col = "ground_truth")


## ----spatialPlot5,eval=TRUE---------------------------------------------------

# generate a visual index in a pdf file directly
spatialIndexPlot(scores.red, annotation,  
    label.col = "ground_truth",
    out.file="spatialIndexPlot")


## ----spatialPlot6,eval=TRUE,fig.dim = c(6,4.5)--------------------------------

# statistical association with tissue areas based on correlations
# For display purpose, we only use a subset here
assoc.bsr.corr <- spatialAssociation(scores.red[c(1:17), ],
annotation, label.col = "ground_truth",test = "Spearman")

head(assoc.bsr.corr)

spatialAssociationPlot(assoc.bsr.corr)


## ----inferCells,eval=FALSE,include=FALSE--------------------------------------
# 
# ## Additional functions
# 
# # This is not part of the main workflow for analyzing L-R interactions but
# # we offer convenience functions for inferring cell types. However
# # we recommend using other software packages specifically dedicated to this
# # purpose.
# 
# # Inferring cell types
# 
# 
# data(sdc, package = "BulkSignalR")
# bsrdm <- BSRDataModel(counts = sdc)
# bsrdm <- learnParameters(bsrdm)
# bsrinf <- BSRInference(bsrdm)
# 
# # Common TME cell type signatures
# data(immune.signatures, package = "BulkSignalR")
# unique(immune.signatures$signature)
# immune.signatures <- immune.signatures[immune.signatures$signature %in% c(
#     "B cells", "Dentritic cells", "Macrophages",
#     "NK cells", "T cells", "T regulatory cells"
# ), ]
# data("tme.signatures", package = "BulkSignalR")
# signatures <- rbind(immune.signatures,
# tme.signatures[tme.signatures$signature %in%
# c("Endothelial cells", "Fibroblasts"), ])
# tme.scores <- scoreSignatures(bsrdm, signatures)
# 
# # assign cell types to interactions
# lr2ct <- assignCellTypesToInteractions(bsrdm, bsrinf, tme.scores)
# head(lr2ct)
# 
# # cellular network computation and plot
# g.table <- cellularNetworkTable(lr2ct)
# 
# gCN <- cellularNetwork(g.table)
# 
# plot(gCN, edge.width = 5 * E(gCN)$score)
# 
# gSummary <- summarizedCellularNetwork(g.table)
# plot(gSummary, edge.width = 1 + 30 * E(gSummary)$score)
# 
# # relationship with partial EMT---
# # Should be tested HNSCC data instead of SDC!!
# 
# # find the ligands
# data(p.EMT, package = "BulkSignalR")
# gs <- p.EMT$gene
# triggers <- relateToGeneSet(bsrinf, gs)
# triggers <- triggers[triggers$n.genes > 1, ] # at least 2 target genes in the gs
# ligands.in.gs <- intersect(triggers$L, gs)
# triggers <- triggers[!(triggers$L %in% ligands.in.gs), ]
# ligands <- unique(triggers$L)
# 
# # link to cell types
# cf <- cellTypeFrequency(triggers, lr2ct, min.n.genes = 2)
# missing <- setdiff(rownames(tme.scores), names(cf$s))
# cf$s[missing] <- 0
# cf$t[missing] <- 0
# 
# op <- par(mar = c(2, 10, 2, 2))
# barplot(cf$s, col = "lightgray", horiz = T, las = 2)
# par(op)
# 
# # random selections based on random gene sets
# qval.thres <- 1e-3
# inter <- LRinter(bsrinf)
# tg <- tgGenes(bsrinf)
# tcor <- tgCorr(bsrinf)
# good <- inter$qval <= qval.thres
# inter <- inter[good, ]
# tg <- tg[good]
# tcor <- tcor[good]
# all.targets <- unique(unlist(tg))
# r.cf <- list()
# for (k in 1:100) { # should 1000 or more
#     r.gs <- sample(all.targets, length(intersect(gs, all.targets)))
#     r.triggers <- relateToGeneSet(bsrinf, r.gs, qval.thres = qval.thres)
#     r.triggers <- r.triggers[r.triggers$n.genes > 1, ]
#     r.ligands.in.gs <- intersect(r.triggers$L, r.gs)
#     r.triggers <- r.triggers[!(r.triggers$L %in% r.ligands.in.gs), ]
#     r <- cellTypeFrequency(r.triggers, lr2ct, min.n.genes = 2)
#     missing <- setdiff(rownames(tme.scores), names(r$s))
#     r$s[missing] <- 0
#     r$t[missing] <- 0
#     o <- order(names(r$t))
#     r$s <- r$s[o]
#     r$t <- r$t[o]
#     r.cf <- c(r.cf, list(r))
# }
# r.m.s <- foreach(i = seq_len(length(r.cf)), .combine = rbind) %do% {
#     r.cf[[i]]$s
# }
# 
# # plot results
# op <- par(mar = c(2, 10, 2, 2))
# boxplot(r.m.s, col = "lightgray", horizontal = T, las = 2)
# pts <- data.frame(x = as.numeric(cf$s[colnames(r.m.s)]), cty = colnames(r.m.s))
# stripchart(x ~ cty, data = pts, add = TRUE, pch = 19, col = "red")
# par(op)
# for (cty in rownames(tme.scores)) {
#     cat(cty, ": P=", sum(r.m.s[, cty] >= cf$s[cty]) / nrow(r.m.s),
#     "\n", sep = "")
# }

## ----session-info-------------------------------------------------------------
sessionInfo()