Skip to content
This repository has been archived by the owner on Jan 24, 2024. It is now read-only.

netAnalysis_river for subset of pathways and order them based on similarity #255

Open
jv-20 opened this issue Jul 27, 2021 · 2 comments
Open

Comments

@jv-20
Copy link

jv-20 commented Jul 27, 2021

Hi @sqjin

Quick question on netAnalysis_river function.

I am trying to generate river plot for a set of pathways of interest and would like to order them based on similarity.

netAnalysis_river(cellchat.obj, pattern = "incoming", do.order=TRUE,signaling = pathways_interest)

However, I receive the following error when using do.order=TRUE.
Error in hclust(d, "ave") : NA/NaN/Inf in foreign function call (arg 10)

Would you be able to give some insights on this please?

Thanks

@sqjin
Copy link
Owner

sqjin commented Aug 3, 2021

@jv-20 Have you figure it out? If not, you can source the following codes.

#' River plot showing the associations of latent patterns with cell groups and ligand-receptor pairs or signaling pathways
#'
#' River (alluvial) plot shows the correspondence between the inferred latent patterns and cell groups as well as ligand-receptor pairs or signaling pathways.
#'
#' The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. The height of each pattern is proportional to the number of its associated cell groups or signaling pathways.
#'
#' Outgoing patterns reveal how the sender cells coordinate with each other as well as how they coordinate with certain signaling pathways to drive communication.
#'
#' Incoming patterns show how the target cells coordinate with each other as well as how they coordinate with certain signaling pathways to respond to incoming signaling.
#'
#' @param object CellChat object
#' @param slot.name the slot name of object that is used to compute centrality measures of signaling networks
#' @param pattern "outgoing" or "incoming"
#' @param cutoff the threshold for filtering out weak links
#' @param sources.use a vector giving the index or the name of source cell groups of interest
#' @param targets.use a vector giving the index or the name of target cell groups of interest
#' @param signaling a character vector giving the name of signaling pathways of interest
#' @param color.use the character vector defining the color of each cell group
#' @param color.use.pattern the character vector defining the color of each pattern
#' @param color.use.signaling the character vector defining the color of each signaling
#' @param do.order whether reorder the cell groups or signaling according to their similarity
#' @param main.title the title of plot
#' @param font.size font size of the text
#' @param font.size.title font size of the title
#' @importFrom methods slot
#' @importFrom stats cutree dist hclust
#' @importFrom grDevices colorRampPalette
#' @importFrom RColorBrewer brewer.pal
#' @import ggalluvial

#' @importFrom ggalluvial geom_stratum geom_flow to_lodes_form

#' @importFrom ggplot2 geom_text scale_x_discrete scale_fill_manual theme ggtitle
#' @importFrom cowplot plot_grid ggdraw draw_label
#' @return
#' @export
#'
#' @examples
netAnalysis_river <- function(object, slot.name = "netP", pattern = c("outgoing","incoming"), cutoff = 0.5,
sources.use = NULL, targets.use = NULL, signaling = NULL,
color.use = NULL, color.use.pattern = NULL, color.use.signaling = "grey50",
do.order = FALSE, main.title = NULL,
font.size = 2.5, font.size.title = 12){
message("Please make sure you have load library(ggalluvial) when running this function")
requireNamespace("ggalluvial")

suppressMessages(require(ggalluvial))

res.pattern <- methods::slot(object, slot.name)$pattern[[pattern]]
data1 = res.pattern$pattern$cell
data2 = res.pattern$pattern$signaling
if (is.null(color.use.pattern)) {
nPatterns <- length(unique(data1$Pattern))
if (pattern == "outgoing") {
color.use.pattern = ggPalette(nPatterns2)[seq(1,nPatterns2, by = 2)]
} else if (pattern == "incoming") {
color.use.pattern = ggPalette(nPatterns2)[seq(2,nPatterns2, by = 2)]
}
}
if (is.null(main.title)) {
if (pattern == "outgoing") {
main.title = "Outgoing communication patterns of secreting cells"
} else if (pattern == "incoming") {
main.title = "Incoming communication patterns of target cells"
}
}

if (is.null(data2)) {
data1$Contribution[data1$Contribution < cutoff] <- 0
plot.data <- data1
nPatterns<-length(unique(plot.data$Pattern))
nCellGroup<-length(unique(plot.data$CellGroup))
if (is.null(color.use)) {
color.use <- scPalette(nCellGroup)
}
if (is.null(color.use.pattern)){
color.use.pattern <- ggPalette(nPatterns)
}

plot.data.long <- to_lodes_form(plot.data, axes = 1:2, id = "connection")
if (do.order) {
  mat = tapply(plot.data[["Contribution"]], list(plot.data[["CellGroup"]], plot.data[["Pattern"]]), sum)
  d <- dist(as.matrix(mat))
  hc <- hclust(d, "ave")
  k <- length(unique(grep("Pattern", plot.data.long$stratum[plot.data.long$Contribution != 0], value = T)))
  cluster <- hc %>% cutree(k)
  order.name <- order(cluster)
  plot.data.long$stratum <- factor(plot.data.long$stratum, levels = c(names(cluster)[order.name], colnames(mat)))
  color.use <- color.use[order.name]
}
color.use.all <- c(color.use, color.use.pattern)
gg <- ggplot(plot.data.long,aes(x = factor(x, levels = c("CellGroup", "Pattern")),y=Contribution,
                                stratum = stratum, alluvium = connection,
                                fill = stratum, label = stratum)) +
  geom_flow(width = 1/3,aes.flow = "backward") +
  geom_stratum(width=1/3,size=0.1,color="black", alpha = 0.8, linetype = 1) +
  geom_text(stat = "stratum", size = font.size) +
  scale_x_discrete(limits = c(),  labels=c("Cell groups", "Patterns")) +
  scale_fill_manual(values = alpha(color.use.all, alpha = 0.8), drop = FALSE) +
  theme_bw()+
  theme(legend.position = "none",
        axis.title = element_blank(),
        axis.text.y= element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor  = element_blank(),
        panel.border = element_blank(),
        axis.ticks = element_blank(),axis.text=element_text(size=10))+
  ggtitle(main.title)

} else {
data1$Contribution[data1$Contribution < cutoff] <- 0
plot.data <- data1
nPatterns<-length(unique(plot.data$Pattern))
nCellGroup<-length(unique(plot.data$CellGroup))
cells.level = levels(object@idents)
if (is.null(color.use)) {
color.use <- scPalette(length(cells.level))[cells.level %in% unique(plot.data$CellGroup)]
}
if (is.null(color.use.pattern)){
color.use.pattern <- ggPalette(nPatterns)
}
if (!is.null(sources.use)) {
if (is.numeric(sources.use)) {
sources.use <- cells.level[sources.use]
}
plot.data <- subset(plot.data, CellGroup %in% sources.use)
}
if (!is.null(targets.use)) {
if (is.numeric(targets.use)) {
targets.use <- cells.level[targets.use]
}
plot.data <- subset(plot.data, CellGroup %in% targets.use)
}
## connect cell groups with patterns
plot.data.long <- to_lodes_form(plot.data, axes = 1:2, id = "connection")
if (do.order) {
mat = tapply(plot.data[["Contribution"]], list(plot.data[["CellGroup"]], plot.data[["Pattern"]]), sum)
d <- dist(as.matrix(mat))
hc <- hclust(d, "ave")
k <- length(unique(grep("Pattern", plot.data.long$stratum[plot.data.long$Contribution != 0], value = T)))
cluster <- hc %>% cutree(k)
order.name <- order(cluster)
plot.data.long$stratum <- factor(plot.data.long$stratum, levels = c(names(cluster)[order.name], colnames(mat)))
color.use <- color.use[order.name]
}
color.use.all <- c(color.use, color.use.pattern)
StatStratum <- ggalluvial::StatStratum
gg1 <- ggplot(plot.data.long,aes(x = factor(x, levels = c("CellGroup", "Pattern")),y=Contribution,
stratum = stratum, alluvium = connection,
fill = stratum, label = stratum)) +
geom_flow(width = 1/3,aes.flow = "backward") +
geom_stratum(width=1/3,size=0.1,color="black", alpha = 0.8, linetype = 1) +
geom_text(stat = "stratum", size = font.size) +
scale_x_discrete(limits = c(), labels=c("Cell groups", "Patterns")) +
scale_fill_manual(values = alpha(color.use.all, alpha = 0.8), drop = FALSE) +
theme_bw()+
theme(legend.position = "none",
axis.title = element_blank(),
axis.text.y= element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.ticks = element_blank(),axis.text=element_text(size=10)) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

## connect patterns with signaling
data2$Contribution[data2$Contribution < cutoff] <- 0
plot.data <- data2
nPatterns<-length(unique(plot.data$Pattern))
nSignaling<-length(unique(plot.data$Signaling))
if (length(color.use.signaling) == 1) {
  color.use.all <- c(color.use.pattern, rep(color.use.signaling, nSignaling))
} else {
  color.use.all <- c(color.use.pattern, color.use.signaling)
}

if (!is.null(signaling)) {
  plot.data <- plot.data[plot.data$Signaling %in% signaling, ]
}

plot.data.long <- ggalluvial::to_lodes_form(plot.data, axes = 1:2, id = "connection")
if (do.order) {
  mat = tapply(plot.data[["Contribution"]], list(plot.data[["Signaling"]], plot.data[["Pattern"]]), sum)
  mat[is.na(mat)] <- 0; mat <- mat[-which(rowSums(mat) == 0), ]
  d <- dist(as.matrix(mat))
  hc <- hclust(d, "ave")
  k <- length(unique(grep("Pattern", plot.data.long$stratum[plot.data.long$Contribution != 0], value = T)))
  cluster <- hc %>% cutree(k)
  order.name <- order(cluster)
  plot.data.long$stratum <- factor(plot.data.long$stratum, levels = c(colnames(mat),names(cluster)[order.name]))
}

gg2 <- ggplot(plot.data.long,aes(x = factor(x, levels = c("Pattern", "Signaling")),y= Contribution,
                                 stratum = stratum, alluvium = connection,
                                 fill = stratum, label = stratum)) +
  geom_flow(width = 1/3,aes.flow = "forward") +
  geom_stratum(width=1/3,size=0.1,color="black", alpha = 0.8, linetype = 1) +
  geom_text(stat = "stratum", size = font.size) + # 2.5
  scale_x_discrete(limits = c(),  labels=c("Patterns", "Signaling")) +
  scale_fill_manual(values = alpha(color.use.all, alpha = 0.8), drop = FALSE) +
  theme_bw()+
  theme(legend.position = "none",
        axis.title = element_blank(),
        axis.text.y= element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor  = element_blank(),
        panel.border = element_blank(),
        axis.ticks = element_blank(),axis.text=element_text(size= 10))+
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

gg <- cowplot::plot_grid(gg1, gg2,align = "h", nrow = 1)
title <- cowplot::ggdraw() + cowplot::draw_label(main.title,size = font.size.title)
gg <- cowplot::plot_grid(title, gg, ncol=1, rel_heights=c(0.1, 1))

}
return(gg)
}

@jv-20
Copy link
Author

jv-20 commented Aug 4, 2021

Hi @sqjin , Many thanks for this. I haven't figured out yet. I will try this and get back.

Thanks again

Kind regards

Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants