sqjin / CellChat

R toolkit for inference, visualization and analysis of cell-cell communication from single-cell data
GNU General Public License v3.0
630 stars 143 forks source link

Arrow missing from netVisual_chord_gene #266

Open aaronkwong opened 3 years ago

aaronkwong commented 3 years ago

Hi CellChat Team,

I am encountering an issue where an arrow appears to be missing when I run netVisual_chord_gene on the "CX3C" pathway. Running other pathways produce chord plots as expected without issue.

When I create the chord diagram, it appears that space is made for the arrow, but one isn't drawn. netVisual_chord_gene(cc_post, signaling = c("CX3C"),legend.pos.x = 8) image

I believe that two interactions should be drawn: df.net_post <- subsetCommunication(cc_post) sig_pathways_post<-df.net_post[df.net_post$pathway_name=="CX3C",] image

Do you have any advice on this issue?

Thank you, AW

Session Info: R version 4.0.3 (2020-10-10) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 18.04.1 LTS

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale: [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] parallel stats graphics grDevices utils datasets methods base

other attached packages: [1] CellChat_1.1.3 Biobase_2.50.0 BiocGenerics_0.36.0 ggplot2_3.3.5 igraph_1.2.6 dplyr_1.0.7 trqwe_0.1

loaded via a namespace (and not attached): [1] matrixStats_0.59.0 fs_1.5.0 usethis_2.0.1 devtools_2.4.2 doParallel_1.0.16 RColorBrewer_1.1-2 rprojroot_2.0.2
[8] tools_4.0.3 irlba_2.3.3 utf8_1.2.2 R6_2.5.0 colorspace_2.0-2 GetoptLong_1.0.5 withr_2.4.2
[15] tidyselect_1.1.1 prettyunits_1.1.1 processx_3.5.2 curl_4.3.2 compiler_4.0.3 cli_3.0.1 Cairo_1.5-12.2
[22] network_1.17.1 desc_1.3.0 pkgmaker_0.32.2 scales_1.1.1 pbapply_1.4-3 callr_3.7.0 NMF_0.23.0
[29] systemfonts_1.0.2 stringr_1.4.0 digest_0.6.27 svglite_2.0.0 pkgconfig_2.0.3 parallelly_1.27.0 sessioninfo_1.1.1
[36] fastmap_1.1.0 rlang_0.4.11 GlobalOptions_0.1.2 rstudioapi_0.13 FNN_1.1.3 shape_1.4.6 generics_0.1.0
[43] jsonlite_1.7.2 statnet.common_4.5.0 magrittr_2.0.1 patchwork_1.1.1 Matrix_1.3-4 Rcpp_1.0.7 munsell_0.5.0
[50] S4Vectors_0.28.1 fansi_0.5.0 reticulate_1.20 lifecycle_1.0.0 stringi_1.7.3 ggalluvial_0.12.3 pkgbuild_1.2.0
[57] plyr_1.8.6 grid_4.0.3 ggrepel_0.9.1 listenv_0.8.0 crayon_1.4.1 lattice_0.20-41 cowplot_1.1.1
[64] circlize_0.4.13 sna_2.6 ComplexHeatmap_2.6.2 ps_1.5.0 pillar_1.6.2 rjson_0.2.20 rngtools_1.5
[71] future.apply_1.8.1 reshape2_1.4.4 codetools_0.2-18 stats4_4.0.3 pkgload_1.2.1 glue_1.4.2 remotes_2.4.0
[78] png_0.1-7 vctrs_0.3.8 foreach_1.5.1 testthat_3.0.4 gtable_0.3.0 purrr_0.3.4 clue_0.3-59
[85] future_1.21.0 assertthat_0.2.1 cachem_1.0.5 gridBase_0.4-7 xtable_1.8-4 RSpectra_0.16-0 coda_0.19-4
[92] tibble_3.1.3 iterators_1.0.13 registry_0.5-1 memoise_2.0.0 IRanges_2.24.1 cluster_2.1.0 globals_0.14.0
[99] ellipsis_0.3.2

sqjin commented 3 years ago

@aaronkwong Interesting! It looks like another target is already drawn. Can you run the source codes of netVisual_chord_gene to check the input for the core plotting function chordDiagram?

aaronkwong commented 3 years ago

@sqjin Thank you for the quick response. I added some print statements before the call to chordDiagram

  cat("\ndf.plot:\n")
  print(df.plot)
  cat("\norder.sector:\n")
  print(order.sector)
  cat("\nedge.color\n")
  print(edge.color)
  cat("\ngrid.col\n")
  print(grid.col)
  cat("\ntransparency\n")
  print(transparency)
  cat("\nlink.border\n")
  print(link.border)
  cat("\ndirectional\n")
  print(directional)
  cat("\nlink.arr.type\n")
  print(link.arr.type)
  cat("\nannotationTrackHeight\n")
  print(annotationTrackHeight)
  cat("\nlist(track.height = max(strwidth(order.sector)))\n")
  print(list(track.height = max(strwidth(order.sector))))
  cat("\nsmall.gap\n")
  print(small.gap)
  cat("\nbig.gap\n")
  print(big.gap)
  cat("\nlink.visible\n")
  print(link.visible)
  cat("\nscale\n")
  print(scale)
  cat("\nlink.target.prop\n")
  print(link.target.prop)
  cat("\nreduce\n")
  print(reduce)
  chordDiagram(df.plot,
               order = order.sector,
               col = edge.color,
               grid.col = grid.col,
               transparency = transparency,
               link.border = link.border,
               directional = directional,
               direction.type = c("diffHeight","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,
               scale = scale,
               link.target.prop = link.target.prop,
               reduce = reduce,
               ...)

Here is the result:

df.plot: ligand receptor prob CX3CL1_CX3CR1.396 CX3CL1 CX3CR1 0.005746460 CX3CL1_CX3CR1.258 CX3CL1 CX3CR1 0.001922853

order.sector: [1] "CX3CL1" "CX3CR1" "CX3CR1 "

edge.color Capillary Capillary "#F781BF" "#F781BF"

grid.col CX3CL1 CX3CR1 CX3CR1
"#F781BF" "#B2DF8A" "#A6CEE3"

transparency [1] 0.4

link.border [1] NA

directional [1] 1

link.arr.type [1] "big.arrow"

annotationTrackHeight [1] 0.03

list(track.height = max(strwidth(order.sector))) $track.height [1] 0.34855

small.gap [1] 1

big.gap [1] 10

link.visible [1] TRUE

scale [1] FALSE

link.target.prop [1] TRUE

reduce [1] -1

sqjin commented 3 years ago

@aaronkwong I am sorry I cannot see any issue.

pedriniedoardo commented 2 years ago

@sqjin I have encountered the same issue. It seems to me it is triggered when there are two targets for one source.

Here is an example using a specific LR pair:

netVisual_chord_gene(object.list[[2]],pairLR.use = data.frame('interaction_name' = "SPP1_CD44"),
                     slot.name = 'net',
                     net = net.up, lab.cex = 0.8, small.gap = 3.5,
                     title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))

image

like in @aaronkwong case, the input matrix suggest there should be two edges.

             source target interaction_name       prob interaction_name_2 pathway_name ligand receptor         annotation       evidence
SPP1_CD44.1     IMM    AST        SPP1_CD44 0.14315963        SPP1 - CD44         SPP1   SPP1     CD44 Secreted Signaling PMID: 21907263
SPP1_CD44.15    IMM    LYM        SPP1_CD44 0.05779804        SPP1 - CD44         SPP1   SPP1     CD44 Secreted Signaling PMID: 21907263

interestingly If I plot the more general signalling pathway (not the specific LR pair) the edge is correctly plotted.

netVisual_chord_gene(object.list[[2]],signaling = "SPP1",
                     slot.name = 'net',
                     net = net.up, lab.cex = 0.8, small.gap = 3.5,
                     title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))

image

But If I try to filter on the source, targets or pairs the issue is triggered.

levels(object.list[[2]]@idents)

[1] "AST"   "IMM"   "LYM"   "NEU"   "OLIGO" "OPC"   "VAS" 

netVisual_chord_gene(object.list[[2]],signaling = "SPP1",
                     targets.use = c(1,3),
                     slot.name = 'net',
                     net = net.up, lab.cex = 0.8, small.gap = 3.5,
                     title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))

image

The same issue is triggered if the target is not that same receptor.

netVisual_chord_gene(object.list[[2]],signaling = "SPP1",
                     targets.use = c(1,2),
                     slot.name = 'net',
                     net = net.up, lab.cex = 0.8, small.gap = 3.5,
                     title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))

image

But If I have more than two targets the issue is not triggered.

netVisual_chord_gene(object.list[[2]],signaling = "SPP1",
                     targets.use = c(1:3),
                     slot.name = 'net',
                     net = net.up, lab.cex = 0.8, small.gap = 3.5,
                     title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))

image

pedriniedoardo commented 2 years ago

@aaronkwong @sqjin after looking into that I think I have found a possible solution. I tried to simplify the example and started from the circlize::chordDiagram function.

I was surprised to see that the issue was also reproducible in the simple call of the function.

df2 = data.frame(from = c("SPP1"), to = c("CD44","CD44 "), value = c(1,2))
circlize::chordDiagram(df2, directional = 1,
                       direction.type = c("diffHeight", "arrows"),
                       link.arr.type = "big.arrow")

image

interestingly also in this case the issue was not presented if the number of target is three

df1 = data.frame(from = c("SPP1"), to = c("CD44","CD44 ","CD44"), value = c(1,2,0))
circlize::chordDiagram(df1, directional = 1,
                       direction.type = c("diffHeight", "arrows"),
                       link.arr.type = "big.arrow")

image

Aftre reading the documentation of the circlize::chordDiagramFromDataFrame function I noticed that the argument direction.type is passed as "diffHeight+arrows" rather than a vector of c("diffHeight", "arrows"). I therefore tried to change it and it seems to be solved.

df2 = data.frame(from = c("SPP1"), to = c("CD44","CD44 "), value = c(1,2))
circlize::chordDiagram(df2, 
                       directional = 1,
                       direction.type = c("diffHeight+arrows"),
                       link.arr.type = "big.arrow")

image

I therefore changed the implementation fo the netVisual_chord_gene function changing the direction.type accordingly and implemented it as netVisual_chord_gene2

netVisual_chord_gene2 <- function (object, slot.name = "net", color.use = NULL, signaling = NULL, 
          pairLR.use = NULL, net = NULL, sources.use = NULL, targets.use = NULL, 
          lab.cex = 0.8, small.gap = 1, big.gap = 10, annotationTrackHeight = c(0.03), 
          link.visible = TRUE, scale = FALSE, directional = 1, link.target.prop = TRUE, 
          reduce = -1, transparency = 0.4, link.border = NA, title.name = NULL, 
          legend.pos.x = 20, legend.pos.y = 20, show.legend = TRUE, 
          thresh = 0.05, ...) 
{
  if (!is.null(pairLR.use)) {
    if (!is.data.frame(pairLR.use)) {
      stop("pairLR.use should be a data frame with a signle column named either 'interaction_name' or 'pathway_name' ")
    }
    else if ("pathway_name" %in% colnames(pairLR.use)) {
      message("slot.name is set to be 'netP' when pairLR.use contains signaling pathways")
      slot.name = "netP"
    }
  }
  if (!is.null(pairLR.use) & !is.null(signaling)) {
    stop("Please do not assign values to 'signaling' when using 'pairLR.use'")
  }
  if (is.null(net)) {
    prob <- slot(object, "net")$prob
    pval <- slot(object, "net")$pval
    prob[pval > thresh] <- 0
    net <- reshape2::melt(prob, value.name = "prob")
    colnames(net)[1:3] <- c("source", "target", "interaction_name")
    pairLR = dplyr::select(object@LR$LRsig, c("interaction_name_2", 
                                              "pathway_name", "ligand", "receptor", "annotation", 
                                              "evidence"))
    idx <- match(net$interaction_name, rownames(pairLR))
    temp <- pairLR[idx, ]
    net <- cbind(net, temp)
  }
  if (!is.null(signaling)) {
    pairLR.use <- data.frame()
    for (i in 1:length(signaling)) {
      pairLR.use.i <- searchPair(signaling = signaling[i], 
                                 pairLR.use = object@LR$LRsig, key = "pathway_name", 
                                 matching.exact = T, pair.only = T)
      pairLR.use <- rbind(pairLR.use, pairLR.use.i)
    }
  }
  if (!is.null(pairLR.use)) {
    if ("interaction_name" %in% colnames(pairLR.use)) {
      net <- subset(net, interaction_name %in% pairLR.use$interaction_name)
    }
    else if ("pathway_name" %in% colnames(pairLR.use)) {
      net <- subset(net, pathway_name %in% as.character(pairLR.use$pathway_name))
    }
  }
  if (slot.name == "netP") {
    net <- dplyr::select(net, c("source", "target", "pathway_name", 
                                "prob"))
    net$source_target <- paste(net$source, net$target, sep = "sourceTotarget")
    net <- net %>% dplyr::group_by(source_target, pathway_name) %>% 
      dplyr::summarize(prob = sum(prob))
    a <- stringr::str_split(net$source_target, "sourceTotarget", 
                            simplify = T)
    net$source <- as.character(a[, 1])
    net$target <- as.character(a[, 2])
    net$ligand <- net$pathway_name
    net$receptor <- " "
  }
  if (!is.null(sources.use)) {
    if (is.numeric(sources.use)) {
      sources.use <- levels(object@idents)[sources.use]
    }
    net <- subset(net, source %in% sources.use)
  }
  else {
    sources.use <- levels(object@idents)
  }
  if (!is.null(targets.use)) {
    if (is.numeric(targets.use)) {
      targets.use <- levels(object@idents)[targets.use]
    }
    net <- subset(net, target %in% targets.use)
  }
  else {
    targets.use <- levels(object@idents)
  }
  df <- subset(net, prob > 0)
  if (nrow(df) == 0) {
    stop("No signaling links are inferred! ")
  }
  if (length(unique(net$ligand)) == 1) {
    message("You may try the function `netVisual_chord_cell` for visualizing individual signaling pathway")
  }
  df$id <- 1:nrow(df)
  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, -prob)), ]
  df.ordered.target <- df[with(df, order(target, -prob)), ]
  order.source <- unique(df.ordered.source[, c("ligand", "source")])
  order.target <- unique(df.ordered.target[, c("receptor", 
                                               "target")])
  order.sector <- c(order.source$ligand, order.target$receptor)
  if (is.null(color.use)) {
    color.use = scPalette(nlevels(object@idents))
    names(color.use) <- levels(object@idents)
    color.use <- color.use[levels(object@idents) %in% as.character(union(df$source, 
                                                                         df$target))]
  }
  else if (is.null(names(color.use))) {
    names(color.use) <- levels(object@idents)
    color.use <- color.use[levels(object@idents) %in% as.character(union(df$source, 
                                                                         df$target))]
  }
  edge.color <- color.use[as.character(df.ordered.source$source)]
  names(edge.color) <- as.character(df.ordered.source$source)
  grid.col.ligand <- color.use[as.character(order.source$source)]
  names(grid.col.ligand) <- as.character(order.source$source)
  grid.col.receptor <- color.use[as.character(order.target$target)]
  names(grid.col.receptor) <- as.character(order.target$target)
  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"
  }
  circos.clear()
  chordDiagram(df.plot, order = order.sector, col = edge.color, 
               grid.col = grid.col, transparency = transparency, link.border = link.border, 
               directional = directional, direction.type = c("diffHeight+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, 
               scale = scale, link.target.prop = link.target.prop, reduce = reduce, 
               ...)
  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)
  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"))
  }
  circos.clear()
  if (!is.null(title.name)) {
    text(-0, 1.02, title.name, cex = 1)
  }
  gg <- recordPlot()
  return(gg)
}

now it is working as expected

netVisual_chord_gene2(object.list[[2]],
                     slot.name = 'net',
                     net = net.up %>%
                       filter(ligand=="SPP1",
                              receptor=="CD44"), lab.cex = 0.8, small.gap = 3.5,
                     title.name = paste0("significant signaling in ", names(object.list)[2]))

image

DHelix commented 1 year ago

@pedriniedoardo, thanks a lot!!!