diff --git a/R/analysis.R b/R/analysis.R index b2d9ada..f70963f 100644 --- a/R/analysis.R +++ b/R/analysis.R @@ -1,4 +1,3 @@ - #' Compute and visualize the contribution of each ligand-receptor pair in the overall signaling pathways #' #' @param object CellChat object @@ -1029,14 +1028,14 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m measure <- match.arg(measure) mode <- match.arg(mode) options(warn = -1) - object.names <- names(methods::slot(object, slot.name)) + object1 <- methods::slot(object, slot.name) + object.names <- names(object1) if (measure == "weight") { ylabel = "Information flow" } else if (measure == "count") { ylabel = "Number of interactions" } if (mode == "single") { - object1 <- methods::slot(object, slot.name) prob = object1$prob prob[object1$pval > thresh] <- 0 if (measure == "count") { @@ -1071,12 +1070,8 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m pSum <- apply(prob, 3, sum) pSum.original <- pSum if (measure == "weight") { - pSum <- -1/log(pSum) + pSum <- log(1+pSum) pSum[is.na(pSum)] <- 0 - idx1 <- which(is.infinite(pSum) | pSum < 0) - values.assign <- seq(max(pSum)*1.1, max(pSum)*1.5, length.out = length(idx1)) - position <- sort(pSum.original[idx1], index.return = TRUE)$ix - pSum[idx1] <- values.assign[match(1:length(idx1), position)] } else if (measure == "count") { pSum <- pSum.original } @@ -1087,12 +1082,7 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m idx <- with(df, order(df$contribution)) df <- df[idx, ] df$name <- factor(df$name, levels = as.character(df$name)) - for (i in 1:length(pair.name)) { - df.t <- df[df$name == pair.name[i], "contribution"] - if (sum(df.t) == 0) { - df <- df[-which(df$name == pair.name[i]), ] - } - } + df <- df[ave(df$contribution, df$name, FUN = sum) != 0, ] if (!is.null(signaling.type)) { LR <- subset(object@DB$interaction, annotation %in% signaling.type) @@ -1163,38 +1153,32 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m } prob.list[[i]] <- prob pSum.original[[i]] <- apply(prob, 3, sum) - if (measure == "weight") { - pSum[[i]] <- -1/log(pSum.original[[i]]) + } + # Step 1: Flatten both lists to a numeric vector + all_values <- unlist(pSum.original, use.names = FALSE) + + # Step 2: Compute global min and max + global_min <- min(all_values, na.rm = TRUE) + global_max <- max(all_values, na.rm = TRUE) + + for (i in 1:length(comparison)) { + if (measure == "weight") { + pSum[[i]] <- log(1+pSum.original[[i]]) pSum[[i]][is.na(pSum[[i]])] <- 0 - idx[[i]] <- which(is.infinite(pSum[[i]]) | pSum[[i]] < 0) - pSum.original.all <- c(pSum.original.all, pSum.original[[i]][idx[[i]]]) } else if (measure == "count") { pSum[[i]] <- pSum.original[[i]] # the prob is already binarized in line 1136 } pair.name[[i]] <- names(pSum.original[[i]]) object.names.comparison <- c(object.names.comparison, object.names[comparison[i]]) } - if (measure == "weight") { - values.assign <- seq(max(unlist(pSum))*1.1, max(unlist(pSum))*1.5, length.out = length(unlist(idx))) - position <- sort(pSum.original.all, index.return = TRUE)$ix - for (i in 1:length(comparison)) { - if (i == 1) { - pSum[[i]][idx[[i]]] <- values.assign[match(1:length(idx[[i]]), position)] - } else { - pSum[[i]][idx[[i]]] <- values.assign[match(length(unlist(idx[1:i-1]))+1:length(unlist(idx[1:i])), position)] - } - } - } - - pair.name.all <- as.character(unique(unlist(pair.name))) - df <- list() - for (i in 1:length(comparison)) { - df[[i]] <- data.frame(name = pair.name.all, contribution = 0, contribution.scaled = 0, group = object.names[comparison[i]], row.names = pair.name.all) - df[[i]][pair.name[[i]],3] <- pSum[[i]] - df[[i]][pair.name[[i]],2] <- pSum.original[[i]] - } + df <- lapply(seq_along(comparison), function(i) { + d <- data.frame(name = pair.name.all, contribution = 0, contribution.scaled = 0, group = object.names[comparison[i]], row.names = pair.name.all) + d[pair.name[[i]], "contribution.scaled"] <- pSum[[i]] + d[pair.name[[i]], "contribution"] <- pSum.original[[i]] + d + }) # contribution.relative <- as.numeric(format(df[[length(comparison)]]$contribution/abs(df[[1]]$contribution), digits=1)) @@ -1214,23 +1198,13 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m contribution.relative[[i]][is.na(contribution.relative[[i]])] <- 0 } names(contribution.relative) <- paste0("contribution.relative.", 1:length(contribution.relative)) - for (i in 1:length(comparison)) { - for (j in 1:length(contribution.relative)) { - df[[i]][[names(contribution.relative)[j]]] <- contribution.relative[[j]] - } + for (j in seq_along(contribution.relative)) { + df <- lapply(df, function(x) { x[[names(contribution.relative)[j]]] <- contribution.relative[[j]]; x }) } df[[1]]$contribution.data2 <- df[[length(comparison)]]$contribution - if (length(comparison) == 2) { - idx <- with(df[[1]], order(-contribution.relative.1, contribution, -contribution.data2)) - } else if (length(comparison) == 3) { - idx <- with(df[[1]], order(-contribution.relative.1, -contribution.relative.2,contribution, -contribution.data2)) - } else if (length(comparison) == 4) { - idx <- with(df[[1]], order(-contribution.relative.1, -contribution.relative.2, -contribution.relative.3, contribution, -contribution.data2)) - } else { - idx <- with(df[[1]], order(-contribution.relative.1, -contribution.relative.2, -contribution.relative.3, -contribution.relative.4, contribution, -contribution.data2)) - } - - + order_args <- c(paste0("-contribution.relative.", seq_along(contribution.relative)), "contribution", "-contribution.data2") + order_expr <- paste(order_args, collapse = ", ") + idx <- eval(parse(text = paste0("with(df[[1]], order(", order_expr, "))"))) for (i in 1:length(comparison)) { df[[i]] <- df[[i]][idx, ] @@ -1242,13 +1216,9 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m df$group <- factor(df$group, levels = object.names.comparison) if (is.null(color.use)) { - color.use = ggPalette(length(comparison)) + color.use <- ggPalette(length(levels(df$group))) } - # https://stackoverflow.com/questions/49448497/coord-flip-changes-ordering-of-bars-within-groups-in-grouped-bar-plot - df$group <- factor(df$group, levels = rev(levels(df$group))) - color.use <- rev(color.use) - # perform statistical analysis # if (do.stat) { # pvalues <- c() @@ -1272,30 +1242,30 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m # df$pvalues <- pvalues # } if (do.stat & length(comparison) == 2) { - for (i in 1:length(pair.name.all)) { - if (nrow(prob.list[[j]]) != nrow(prob.list[[1]])) { - if (paired.test) { - stop("Paired test is not applicable to datasets with different cellular compositions! Please set `do.stat = FALSE` or `paired.test = FALSE`! \n") - } + if (nrow(prob.list[[1]]) != nrow(prob.list[[2]])) { + if (paired.test) { + stop("Paired test is not applicable to datasets with different cellular compositions! Please set `do.stat = FALSE` or `paired.test = FALSE`! \n") + } } - prob.values <- matrix(0, nrow = nrow(prob.list[[1]]) * nrow(prob.list[[1]]), ncol = length(comparison)) - for (j in 1:length(comparison)) { - if (pair.name.all[i] %in% pair.name[[j]]) { - prob.values[, j] <- as.vector(prob.list[[j]][ , , pair.name.all[i]]) + for (i in 1:length(pair.name.all)) { + prob.values <- matrix(0, nrow = nrow(prob.list[[1]]) * nrow(prob.list[[1]]), ncol = length(comparison)) + for (j in 1:length(comparison)) { + if (pair.name.all[i] %in% pair.name[[j]]) { + prob.values[, j] <- as.vector(prob.list[[j]][ , , pair.name.all[i]]) + } else { + prob.values[, j] <- NA + } + } + prob.values <- prob.values[rowSums(prob.values, na.rm = TRUE) != 0, , drop = FALSE] + if (nrow(prob.values) >3 & sum(is.na(prob.values)) == 0) { + pvalues <- wilcox.test(prob.values[ ,1], prob.values[ ,2], paired = paired.test)$p.value } else { - prob.values[, j] <- NA + pvalues <- 0 } + pvalues[is.na(pvalues)] <- 0 + df$pvalues[df$name == pair.name.all[i]] <- pvalues } - prob.values <- prob.values[rowSums(prob.values, na.rm = TRUE) != 0, , drop = FALSE] - if (nrow(prob.values) >3 & sum(is.na(prob.values)) == 0) { - pvalues <- wilcox.test(prob.values[ ,1], prob.values[ ,2], paired = paired.test)$p.value - } else { - pvalues <- 0 - } - pvalues[is.na(pvalues)] <- 0 - df$pvalues[df$name == pair.name.all[i]] <- pvalues } - } if (length(comparison) == 2) { @@ -1309,12 +1279,7 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m colors.text = NULL } - for (i in 1:length(pair.name.all)) { - df.t <- df[df$name == pair.name.all[i], "contribution"] - if (sum(df.t) == 0) { - df <- df[-which(df$name == pair.name.all[i]), ] - } - } + df <- df[ave(df$contribution, df$name, FUN = sum) != 0, ] if ((slot.name == "netP") && (!is.null(signaling))) { df <- subset(df, name %in% signaling) @@ -1372,7 +1337,7 @@ rankNet <- function(object, slot.name = "netP", measure = c("weight","count"), m } gg <- gg + theme(axis.text=element_text(size=font.size), axis.title.y = element_text(size=font.size)) - gg <- gg + scale_fill_manual(name = "", values = color.use) + gg <- gg + scale_fill_manual(name = "", values = setNames(color.use, rev(levels(df$group)))) gg <- gg + guides(fill = guide_legend(reverse = TRUE)) gg <- gg + theme(axis.text.x = element_text(angle = x.angle, hjust=x.hjust), axis.text.y = element_text(angle = y.angle, hjust=y.hjust)) diff --git a/R/visualization.R b/R/visualization.R index 96abb61..d3987db 100644 --- a/R/visualization.R +++ b/R/visualization.R @@ -1900,9 +1900,8 @@ netVisual_heatmap <- function(object, comparison = c(1,2), measure = c("count", } df.net <- subset(df.net, target %in% targets.use) } - cells.level <- rownames(net.diff) - df.net$source <- factor(df.net$source, levels = cells.level) - df.net$target <- factor(df.net$target, levels = cells.level) + df.net$source <- factor(df.net$source, levels = sources.use) + df.net$target <- factor(df.net$target, levels = targets.use) df.net$value[is.na(df.net$value)] <- 0 net <- tapply(df.net[["value"]], list(df.net[["source"]], df.net[["target"]]), sum) } @@ -2294,15 +2293,16 @@ netVisual_bubble <- function(object, sources.use = NULL, targets.use = NULL, sig group.names <- paste0(group.names0, " (", dataset.name[comparison[ii]], ")") if (nrow(df.net) > 0) { - df.net$pval[df.net$pval > 0.05] = 1 - df.net$pval[df.net$pval > 0.01 & df.net$pval <= 0.05] = 2 - df.net$pval[df.net$pval <= 0.01] = 3 + df.net$pval_categorie <- 1 + df.net$pval_categorie[df.net$pval > 0.05] = 1 + df.net$pval_categorie[df.net$pval > 0.01 & df.net$pval <= 0.05] = 2 + df.net$pval_categorie[df.net$pval <= 0.01] = 3 df.net$prob[df.net$prob == 0] <- NA df.net$prob.original <- df.net$prob df.net$prob <- -1/log(df.net$prob) } else { df.net <- as.data.frame(matrix(NA, nrow = length(group.names), ncol = 5)) - colnames(df.net) <- c("interaction_name_2","source.target","prob","pval","prob.original") + colnames(df.net) <- c("interaction_name_2","source.target","prob","pval","prob.original", "pval_categorie") df.net$source.target <- group.names0 } # df.net$group.names <- sub(paste0(' \\(',dataset.name[comparison[ii]],'\\)'),'',as.character(df.net$source.target)) @@ -2428,7 +2428,7 @@ netVisual_bubble <- function(object, sources.use = NULL, targets.use = NULL, sig } } - g <- ggplot(df, aes(x = source.target, y = interaction_name_2, color = prob, size = pval)) + + g <- ggplot(df, aes(x = source.target, y = interaction_name_2, color = prob, size = pval_categorie)) + geom_point(pch = 16) + theme_linedraw() + theme(panel.grid.major = element_blank()) + theme(axis.text.x = element_text(angle = angle.x, hjust= hjust.x, vjust = vjust.x), @@ -2438,12 +2438,12 @@ netVisual_bubble <- function(object, sources.use = NULL, targets.use = NULL, sig values <- c(1,2,3); names(values) <- c("p > 0.05", "0.01 < p < 0.05","p < 0.01") if (is.null(dot.size.max)) { - dot.size.max = max(df$pval) + dot.size.max = max(df$pval_categorie) } if (is.null(dot.size.min)) { - dot.size.min = min(df$pval) + dot.size.min = min(df$pval_categorie) } - g <- g + scale_radius(range = c(dot.size.min, dot.size.max), breaks = sort(unique(df$pval)),labels = names(values)[values %in% sort(unique(df$pval))], name = "p-value") + g <- g + scale_radius(range = c(dot.size.min, dot.size.max), breaks = sort(unique(df$pval_categorie)),labels = names(values)[values %in% sort(unique(df$pval_categorie))], name = "p-value") #g <- g + scale_radius(range = c(1,3), breaks = values,labels = names(values), name = "p-value") if (min(df$prob, na.rm = T) != max(df$prob, na.rm = T)) { g <- g + scale_colour_gradientn(colors = colorRampPalette(color.use)(99), na.value = "white", limits=c(quantile(df$prob, 0,na.rm= T), quantile(df$prob, 1,na.rm= T)), @@ -2927,37 +2927,26 @@ netVisual_chord_gene <- function(object, slot.name = "net", color.use = NULL, } df$id <- 1:nrow(df) - # deal with duplicated sector names - ligand.uni <- unique(df$ligand) - for (i in 1:length(ligand.uni)) { - df.i <- df[df$ligand == ligand.uni[i], ] - source.uni <- unique(df.i$source) - for (j in 1:length(source.uni)) { - df.i.j <- df.i[df.i$source == source.uni[j], ] - df.i.j$ligand <- paste0(df.i.j$ligand, paste(rep(' ',j-1),collapse = '')) - df$ligand[df$id %in% df.i.j$id] <- df.i.j$ligand - } - } - receptor.uni <- unique(df$receptor) - for (i in 1:length(receptor.uni)) { - df.i <- df[df$receptor == receptor.uni[i], ] - target.uni <- unique(df.i$target) - for (j in 1:length(target.uni)) { - df.i.j <- df.i[df.i$target == target.uni[j], ] - df.i.j$receptor <- paste0(df.i.j$receptor, paste(rep(' ',j-1),collapse = '')) - df$receptor[df$id %in% df.i.j$id] <- df.i.j$receptor - } - } cell.order.sources <- levels(object@idents)[levels(object@idents) %in% sources.use] cell.order.targets <- levels(object@idents)[levels(object@idents) %in% targets.use] df$source <- factor(df$source, levels = cell.order.sources) df$target <- factor(df$target, levels = cell.order.targets) - # df.ordered.source <- df[with(df, order(source, target, -prob)), ] - # df.ordered.target <- df[with(df, order(target, source, -prob)), ] - df.ordered.source <- df[with(df, order(source, -prob)), ] - df.ordered.target <- df[with(df, order(target, -prob)), ] + + df$ligand_std <- toupper(gsub("_", "", df$ligand)) + df$receptor_std <- toupper(gsub("_", "", df$receptor)) + + df.ordered.source <- df %>% + mutate( + ligand_std = toupper(gsub("_", "", ligand)), + receptor_std = toupper(gsub("_", "", receptor)) + ) %>% + arrange(ligand_std, source, desc(prob)) + + df.ordered.target <- df.ordered.source %>% + arrange(receptor_std, target, desc(prob)) + df.ordered.target <- df.ordered.source %>% arrange(receptor, target, desc(prob)) order.source <- unique(df.ordered.source[ ,c('ligand','source')]) order.target <- unique(df.ordered.target[ ,c('receptor','target')]) @@ -2987,54 +2976,182 @@ netVisual_chord_gene <- function(object, slot.name = "net", color.use = NULL, grid.col <- c(as.character(grid.col.ligand), as.character(grid.col.receptor)) names(grid.col) <- order.sector - df.plot <- df.ordered.source[ ,c('ligand','receptor','prob')] - if (directional == 2) { link.arr.type = "triangle" } else { link.arr.type = "big.arrow" } + + df_source <- df.ordered.source %>% + mutate(ligand_std = toupper(gsub("_", "", ligand))) %>% + group_by(ligand, source, receptor, target) %>% + summarise(cluster_value = sum(prob), .groups = "drop") %>% + group_by(ligand) %>% + mutate(fraction = cluster_value / sum(cluster_value)) %>% + ungroup() %>% + arrange(toupper(gsub("_", "", ligand)), receptor, source, target) + + df_target <- df.ordered.source %>% + mutate(receptor_std = toupper(gsub("_", "", receptor))) %>% + group_by(ligand, source, receptor, target) %>% + summarise(cluster_value = sum(prob), .groups = "drop") %>% + group_by(receptor) %>% + mutate(fraction = cluster_value / sum(cluster_value)) %>% + ungroup() %>% + arrange(toupper(gsub("_", "", receptor)), ligand, source, target) + + df.plot <- df.ordered.source[ ,c('source', 'target', 'ligand','receptor','prob')] %>% + arrange(ligand, source, receptor, target) + + + #sector order + order.source <- unique(df.ordered.source[, c("ligand", "source")]) + order.target <- unique(df.ordered.target[, c("receptor", "target")]) + + # This defines order grouped by source/target clusters + ligand_order <- order.source$ligand + receptor_order <- order.target$receptor + + order.sector <- c(ligand_order, receptor_order) + circos.clear() - chordDiagram(df.plot, + + grid::grid.newpage() + vp_main <- grid::viewport(x = 0.5, y = 0.5, width = 1, height = 1, name = "main_vp") + grid::pushViewport(vp_main) + + # Helper to define layout positions + #vplayout <- function(x, y) grid::viewport(layout.pos.row = x, layout.pos.col = y) + + if (!is.null(title.name)) { + grid::grid.text( + label = title.name, + x = unit(0.5, "npc"), + y = unit(0.01, "npc"), + gp = grid::gpar(fontsize = 14, fontface = "bold") + ) + } + + n <- length(unique(names(grid.col))) + + # Small gaps within group, large gap between + small.gap <- 1 + big.gap <- 8 # visibly separate source/target + + gap.after <- c( + rep(small.gap, length(ligand_order) - 1), + big.gap, # gap between ligand and receptor + rep(small.gap, length(receptor_order) - 1), + big.gap # closing the circle + ) + + circos.par( + start.degree = 0, # or whatever layout you prefer + #gap.after = gap.after, # wider gaps at group breaks if needed + track.margin = c(0.01, 0.01), + canvas.xlim = c(-1.2, 1.2), + canvas.ylim = c(-1.2, 1.2) # expand canvas to give room + ) + + chordDiagram(df.plot[ ,c('ligand','receptor','prob')], order = order.sector, col = edge.color, - grid.col = grid.col, - transparency = transparency, + #grid.col = grid.col, + transparency = 0.4, link.border = link.border, directional = directional, - direction.type = c("diffHeight","arrows"), + direction.type = c("arrows"), link.arr.type = link.arr.type, annotationTrack = "grid", - annotationTrackHeight = annotationTrackHeight, - preAllocateTracks = list(track.height = max(strwidth(order.sector))), - small.gap = small.gap, - big.gap = big.gap, - link.visible = link.visible, + preAllocateTracks = list(track.height = 0.05), + link.visible = TRUE, scale = scale, - link.target.prop = link.target.prop, - reduce = reduce, - ...) + link.target.prop = TRUE, + reduce = -1) + + circos.track( + ylim = c(0, 1), + track.index = 2, # ← key change: force higher track level + bg.border = NA, + panel.fun = function(x, y) { + sector_name <- get.cell.meta.data("sector.index") + xlim <- get.cell.meta.data("xlim") + + if ((xlim[2] - xlim[1]) < 1e-4) return(NULL) + + if (sector_name %in% df_target$receptor) { + df_sector <- df_target[df_target$receptor == sector_name, ] %>% arrange(desc(row_number())) + x_start <- xlim[1] + x_end <- xlim[2] + width <- x_end - x_start + + df_sector$fraction <- df_sector$fraction / sum(df_sector$fraction) + + for (i in seq_len(nrow(df_sector))) { + xleft <- x_start + xright <- min(x_start + width * df_sector$fraction[i], x_end) + circos.rect(xleft, 0, xright, 1, + col = custom_colors$discrete[as.character(df_sector$target[i])], + border = NA) + x_start <- xright + } + } + } + ) + + circos.track( + ylim = c(0, 1), + track.index = 2, # ← key change: force higher track level + bg.border = NA, + panel.fun = function(x, y) { + sector_name <- get.cell.meta.data("sector.index") + xlim <- get.cell.meta.data("xlim") + + if ((xlim[2] - xlim[1]) < 1e-4) return(NULL) + + if (sector_name %in% df_source$ligand) { + df_sector <- df_source[df_source$ligand == sector_name, ] %>% arrange(desc(receptor)) + x_start <- xlim[1] + x_end <- xlim[2] + width <- x_end - x_start + + df_sector$fraction <- df_sector$fraction / sum(df_sector$fraction) + + for (i in seq_len(nrow(df_sector))) { + xleft <- x_start + xright <- min(x_start + width * df_sector$fraction[i], x_end) + circos.rect(xleft, 0, xright, 1, + col = custom_colors$discrete[as.character(df_sector$source[i])], + border = NA) + x_start <- xright + } + } + } + ) circos.track(track.index = 1, panel.fun = function(x, y) { - xlim = get.cell.meta.data("xlim") - xplot = get.cell.meta.data("xplot") - ylim = get.cell.meta.data("ylim") - sector.name = get.cell.meta.data("sector.index") - circos.text(mean(xlim), ylim[1], sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5),cex = lab.cex) - }, bg.border = NA) + xlim = get.cell.meta.data("xlim") + xplot = get.cell.meta.data("xplot") + ylim = get.cell.meta.data("ylim") + sector.name = get.cell.meta.data("sector.index") + circos.text(mean(xlim), ylim[1], sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5),cex = lab.cex) + }, bg.border = NA) - # https://jokergoo.github.io/circlize_book/book/legends.html if (show.legend) { - lgd <- ComplexHeatmap::Legend(at = names(color.use), type = "grid", legend_gp = grid::gpar(fill = color.use), title = "Cell State") - ComplexHeatmap::draw(lgd, x = unit(1, "npc")-unit(legend.pos.x, "mm"), y = unit(legend.pos.y, "mm"), just = c("right", "bottom")) - } - + lgd <- ComplexHeatmap::Legend( + at = names(color.use), + type = "grid", + legend_gp = grid::gpar(fill = color.use), + title = "Cell State" + ) + ComplexHeatmap::draw( + lgd, + x = unit(1, "npc") - unit(5, "mm"), # right side + y = unit(0.5, "npc"), # centered vertically + just = c("right", "center") + ) + } circos.clear() - if(!is.null(title.name)){ - text(-0, 1.02, title.name, cex=1) - } - gg <- recordPlot() - return(gg) }