## ----setup, include = FALSE--------------------------------------------------- knitr::knit_hooks$set(optipng = knitr::hook_optipng) ## ----load-libs, message = FALSE, warning = FALSE, results = FALSE------------ library(BulkSignalR) library(igraph) library(dplyr) library(STexampleData) ## ----parallel, message=FALSE, warnings=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,width=4,height=3) ## ----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,width=6,height=2) ## ----heatmapMulti,eval=TRUE,fig.dim = c(7,14)--------------------------------- pathway1 <- pathways(bsrsig.redPBP)[1] signatureHeatmaps( pathway = pathway1, bsrdm = bsrdm, bsrsig = bsrsig.redPBP, h.width = 6, h.height = 8, fontsize = 4, show_column_names = TRUE ) ## ----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 function # # 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) # simpleHeatmap(scoresLR[1:20,],column.names=TRUE, # width=9, height=5, pointsize=8) ## ----spatial1,eval=TRUE,message=FALSE----------------------------------------- # load data ================================================= # We re-initialise the environment variable of Reactome # whith the whole set. (because we previously made a subset before) 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 reduction operations 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 visual index on disk in pdf file 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 LR interactions but # #we offer convenient functions for inferring cell types. However # #we recommend using other softwares 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()