In saeyslab/nichenetr: NicheNet: Modeling Intercellular Communication by Linking Ligands to Target Genes

knitr::opts_chunk$set(
  collapse = TRUE,
  # comment = "#>",
  warning = FALSE,
  message = FALSE
)

This vignette shows how to visualize the output of a NicheNet analysis in a circos plot. This is an additional demonstration to Seurat Wrapper + circos visualization. Here, we will use a normal matrix as input instead of a Seurat object.

We will use the same dataset as from NicheNet's ligand activity analysis on a gene set of interest @puram_single-cell_2017. In contrast to the basic vignette, we will look at communication between multiple cell types. More specifically, we will predict which ligands expressed by both CAFs and endothelial cells can induce the p-EMT program in neighboring malignant cells.

Load packages

library(nichenetr)
library(tidyverse)
library(circlize)

Read in NicheNet's networks

lr_network <- readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))
ligand_target_matrix <- readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds"))
weighted_networks <- readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds"))

lr_network <- lr_network %>% distinct(from, to)

Read in the expression data of interacting cells

This is publicly available single-cell data from CAF and malignant cells from HNSCC tumors.

hnscc_expression <- readRDS(url("https://zenodo.org/record/3260758/files/hnscc_expression.rds"))
expression <- hnscc_expression$expression
sample_info <- hnscc_expression$sample_info # contains meta-information about the cells

# Convert gene names
colnames(expression) <- convert_alias_to_symbols(colnames(expression), "human", verbose = FALSE)

Define a set of potential ligands

tumors_remove <- c("HN10","HN","HN12", "HN13", "HN24", "HN7", "HN8","HN23")

CAF_ids <- sample_info %>%
  filter(`Lymph node` == 0 & !(tumor %in% tumors_remove) &
           `non-cancer cell type` == "CAF") %>% pull(cell)
endothelial_ids <- sample_info %>%
  filter(`Lymph node` == 0 & !(tumor %in% tumors_remove) &
           `non-cancer cell type` == "Endothelial") %>% pull(cell)
malignant_ids <- sample_info %>% filter(`Lymph node` == 0 &
                                          !(tumor %in% tumors_remove) &
                                          `classified  as cancer cell` == 1) %>% pull(cell)


# Define expressed genes in CAFs, endothelial cells and malignant cells
expressed_genes_CAFs <- expression[CAF_ids,] %>%
  apply(2,function(x){10*(2**x - 1)}) %>%
  apply(2,function(x){log2(mean(x) + 1)}) %>% .[. >= 4] %>%
  names()

expressed_genes_endothelial <- expression[endothelial_ids,] %>%
  apply(2,function(x){10*(2**x - 1)}) %>%
  apply(2,function(x){log2(mean(x) + 1)}) %>% .[. >= 4] %>%
  names()

expressed_genes_malignant <- expression[malignant_ids,] %>%
  apply(2,function(x){10*(2**x - 1)}) %>%
  apply(2,function(x){log2(mean(x) + 1)}) %>% .[. >= 4] %>%
  names()

# Define expressed ligands and receptors
ligands <- lr_network %>% pull(from) %>% unique()
expressed_ligands_CAFs <- intersect(ligands,expressed_genes_CAFs)
expressed_ligands_endothelial <- intersect(ligands,expressed_genes_endothelial)
expressed_ligands <- union(expressed_ligands_CAFs, expressed_genes_endothelial)

receptors <- lr_network %>% pull(to) %>% unique()
expressed_receptors <- intersect(receptors,expressed_genes_malignant)

# Define potential ligands
potential_ligands <- lr_network %>%
  filter(from %in% expressed_ligands & to %in% expressed_receptors) %>%
  pull(from) %>% unique()

Define the gene set of interest and background of genes

We will use the p-EMT gene set defined by Puram et al. as gene set of interest and use all genes expressed in malignant cells as background of genes.

pemt_geneset <- readr::read_tsv(url("https://zenodo.org/record/3260758/files/pemt_signature.txt"),
                                col_names = "gene") %>%
  pull(gene) %>% .[. %in% rownames(ligand_target_matrix)]

background_expressed_genes <- expressed_genes_malignant %>%
  .[. %in% rownames(ligand_target_matrix)]

Perform NicheNet's ligand activity analysis

With the ligand activity analysis, we assess how well each potential ligand can predict the p-EMT gene set compared to the background of expressed genes. Note that we combine the ligands from CAFs and endothelial cells in one ligand activity analysis now. Later on, we will look which of the top-ranked ligands is mainly expressed by which of both cell types.

ligand_activities <-  predict_ligand_activities(geneset = pemt_geneset,
                                                background_expressed_genes = background_expressed_genes,
                                                ligand_target_matrix = ligand_target_matrix,
                                                potential_ligands = potential_ligands)

Now, we want to rank the ligands based on their ligand activity (AUPR). We will choose the top 20 ligands here - as opposed to the top 30 in the main vignette - to avoid overcrowding the circos plot.

ligand_activities %>% arrange(-aupr_corrected) 
best_upstream_ligands <- ligand_activities %>%
  top_n(20, aupr_corrected) %>%
  arrange(-aupr_corrected) %>%
  pull(test_ligand)

head(best_upstream_ligands)

Determine now which prioritized ligands are expressed by CAFs and or endothelial cells.

# There is a lot of overlap between both cell types in terms of expressed ligands
intersect(best_upstream_ligands %>% intersect(expressed_ligands_CAFs),
          best_upstream_ligands %>% intersect(expressed_ligands_endothelial))

# Therefore, determine which ligands are more strongly expressed in which of the two
# Calculate average expression of each ligand in CAFs and endothelial cells
ligand_expression_tbl <- tibble(
  ligand = best_upstream_ligands, 
  CAF = expression[CAF_ids,best_upstream_ligands] %>%
    apply(2,function(x){10*(2**x - 1)}) %>%
    apply(2,function(x){log2(mean(x) + 1)}),
  endothelial = expression[endothelial_ids,best_upstream_ligands] %>%
    apply(2,function(x){10*(2**x - 1)}) %>%
    apply(2,function(x){log2(mean(x) + 1)}))

head(ligand_expression_tbl)

# Assign ligand to a cell type based on which cell type has much higher average expression
# If ligand is not expressed more highly enough in either cell type, assign to general category
CAF_specific_ligands <- ligand_expression_tbl %>% filter(CAF > endothelial + 2) %>% pull(ligand)
endothelial_specific_ligands <- ligand_expression_tbl %>% filter(endothelial > CAF + 2) %>% pull(ligand)
general_ligands <- setdiff(best_upstream_ligands,c(CAF_specific_ligands,endothelial_specific_ligands))

ligand_type_indication_df <- tibble(
  ligand_type = c(rep("CAF-specific", times = CAF_specific_ligands %>% length()),
                  rep("General", times = general_ligands %>% length()),
                  rep("Endothelial-specific", times = endothelial_specific_ligands %>% length())),
  ligand = c(CAF_specific_ligands, general_ligands, endothelial_specific_ligands))

head(ligand_type_indication_df)

Infer target genes of top-ranked ligands and visualize in a circos plot

First, get the active ligand-target links by looking which of the p-EMT genes are among the top-predicted target genes for the prioritized ligands:

active_ligand_target_links_df <- best_upstream_ligands %>%
  lapply(get_weighted_ligand_target_links,
         geneset = pemt_geneset,
         ligand_target_matrix = ligand_target_matrix,
         n = 250) %>% bind_rows()

active_ligand_target_links_df <- active_ligand_target_links_df %>%
  mutate(target_type = "p_emt")

Note that you can make a circos plot for multiple gene sets by combining different dataframes and differentiating which target belongs to which gene set via the "target type" column.

To avoid making a circos plots with too many ligand-target links, we will show only links with a weight higher than a predefined cutoff: links belonging to the 66% of lowest scores were removed. Not that this cutoffs and other cutoffs used for this visualization can be changed according to the user's needs.

circos_links <- get_ligand_target_links_oi(ligand_type_indication_df,
                                           active_ligand_target_links_df,
                                           cutoff = 0.66) 

Prepare the circos visualization by giving each segment of ligands and targets a specific color and order, as well as gaps between different cell types. By default, cell types are ordered alphabetically, followed by "General" (then they are drawn counter-clockwise). Users can give a specific order to the cell types by providing a vector of cell types to the argument celltype_order. The gaps between the different segments can also be defined by providing a named list to the argument widths.

ligand_colors <- c("General" = "lawngreen",
                   "CAF-specific" = "royalblue",
                   "Endothelial-specific" = "gold")
target_colors <- c("p_emt" = "tomato")


vis_circos_obj <- prepare_circos_visualization(circos_links,
                                               ligand_colors = ligand_colors,
                                               target_colors = target_colors,
                                               celltype_order = NULL) 

Render the circos plot where all links have the same transparency. Here, only the widths of the blocks that indicate each target gene is proportional the ligand-target regulatory potential (~prior knowledge supporting the regulatory interaction).

make_circos_plot(vis_circos_obj, transparency = FALSE,  args.circos.text = list(cex = 0.5)) 

Render the circos plot where the degree of transparency determined by the regulatory potential value of a ligand-target interaction.

make_circos_plot(vis_circos_obj, transparency = TRUE,  args.circos.text = list(cex = 0.5)) 

To create a legend for the circos plot, we can use the ComplexHeatmap::Legend function and creating a gTree object from it with grid::grid.grabExpr. As the circos plot is drawn on base R graphics (i.e., it is not a ggplot object), we will get the plot using recordPlot().

par(bg = "transparent")

# Default celltype order
celltype_order <- unique(circos_links$ligand_type) %>% sort() %>% .[. != "General"] %>% c(., "General")

# Create legend
circos_legend <- ComplexHeatmap::Legend(
  labels = celltype_order,
  background = ligand_colors[celltype_order],
  type = "point",
  grid_height = unit(3, "mm"),
  grid_width = unit(3, "mm"),
  labels_gp = grid::gpar(fontsize = 8)
  )

circos_legend_grob <- grid::grid.grabExpr(ComplexHeatmap::draw(circos_legend))

make_circos_plot(vis_circos_obj, transparency = TRUE, args.circos.text = list(cex = 0.5))
p_circos_no_legend <- recordPlot()

We can combine the circos plot and the legend using cowplot::plot_grid.

cowplot::plot_grid(p_circos_no_legend, circos_legend_grob, rel_widths = c(1, 0.1))

We can save this plot to an svg file.

svg("ligand_target_circos.svg", width = 10, height = 10)
cowplot::plot_grid(p_circos_no_legend, circos_legend_grob, rel_widths = c(1, 0.1))
dev.off()

Visualize ligand-receptor interactions of the prioritized ligands in a circos plot

To create a ligand-receptor chord diagram, we can perform similar steps as above using the weighted ligand-receptor dataframe instead. However, as as prepare_circos_visualization accesses "target" and "target_type" columns, it is necessary to rename the columns accordingly even though the dataframe contains receptor and not target gene information.

lr_network_top_df <- get_weighted_ligand_receptor_links(best_upstream_ligands,
                                   expressed_receptors,
                                   lr_network,
                                   weighted_networks$lr_sig) %>%
  rename(ligand=from, target=to) %>% 
  mutate(target_type = "p_emt_receptor") %>% 
  inner_join(ligand_type_indication_df)

receptor_colors <- c("p_emt_receptor" = "darkred")

vis_circos_receptor_obj <- prepare_circos_visualization(lr_network_top_df,
                                                        ligand_colors = ligand_colors,
                                                        target_colors = receptor_colors) 

When drawing the plot, the argument link.visible = TRUE is also necessary for making all links visible, since no cutoff is used to filter out ligand-receptor interactions.

make_circos_plot(vis_circos_receptor_obj, transparency = FALSE,
                 link.visible = TRUE,  args.circos.text = list(cex = 0.8)) 

Just as above, if transparency = TRUE, the degree of transparency is determined by the prior interaction weight of the ligand-receptor interaction.

make_circos_plot(vis_circos_receptor_obj, transparency = TRUE,
                 link.visible = TRUE,  args.circos.text = list(cex = 0.8)) 

Adding an outer track to the circos plot (ligand-receptor-target circos plot)

In the paper of Bonnardel, T'Jonck et al. Stellate Cells, Hepatocytes, and Endothelial Cells Imprint the Kupffer Cell Identity on Monocytes Colonizing the Liver Macrophage Niche, we showed in Fig. 6B a ligand-receptor-target circos plot to visualize the main NicheNet predictions. This “ligand-receptor-target” circos plot was made by making first two separate circos plots: the ligand-target and ligand-receptor circos plot. Then these circos plots were overlayed in Inkscape (with the center of the two circles at the same location and the ligand-receptor circos plot bigger than the ligand-target one). To generate the combined circos plot as shown in Fig. 6B, we then manually removed all elements of the ligand-receptor circos plot except the outer receptor layer.

It is also possible to generate this kind of plot programmatically, given that you are able to group your ligands. Below, we demonstrate two approaches that uses either the circlize::draw.sector or circlize::highlight.sector function. The former is more complicated but gives you the flexibility to draw receptor arcs of different lengths, while the highlight.sector function is constrained to the widths of the targets in the inner track.

Using the draw.sector function

For demonstration purposes, we will use a subset of ligands and divide them into four groups that differ in the signaling pathways they interact with. Then, we assign targets and receptors to each ligand group based on their relative rankings (and summed weights in case the rankings are the same). The way we assigned targets and receptors to ligand groups here does not always lead to the most biologically meaningful results, as you will see in the final plot. Hence, it is best if you curate this list manually, e.g., through prior knowledge and by also looking at the ligand-target and ligand-receptor heatmap. Also keep in mind that there is no real correspondence between receptors and targets, as explained more here.

groups <- list(group1 = c("TGFB2", "ENG"),
               group2 = grep("BMP|GDF|INHBA", best_upstream_ligands, value = TRUE),
               group3 = grep("COL|MMP|TIMP", best_upstream_ligands, value = TRUE),
               group4 = "CXCL12")

# Create list of targets for each group of ligand
targets <- lapply(names(groups), function(i) {
  # Rank each target for each ligand
  active_ligand_target_links_df %>% group_by(ligand) %>% mutate(target_rank = dense_rank(desc(weight))) %>%
    filter(ligand %in% groups[[i]]) %>%
    # Make two metrics for each target -> summed weight and avg_rank
    group_by(target) %>%
    summarise(n = n(), summed_weight=sum(weight), avg_rank = mean(target_rank)) %>%
    # Only keep targets that are connected to at least half of ligands in the group
    filter(n > (length(groups[[i]])/2)) %>% mutate(type = i)
}) %>% do.call(rbind, .)

# Check if any targets are distinct per group (groups 2 and 4 have some)
lapply(paste0("group", 1:5), function(group_name) setdiff(targets %>% filter(type == group_name) %>% pull(target), targets %>% filter(type != group_name) %>% pull(target)))

# Assign target to a specific group, first based on ranking and second based on weight
targets_filtered <- targets %>% group_by(target) %>% filter(avg_rank == min(avg_rank)) %>%
  filter(summed_weight == max(summed_weight)) %>% ungroup()

# Do the same for receptors
receptor_colors <- c("#387D7A", "#9DA9A0", "#DD7373", "#725752") %>% setNames(names(groups))

receptors <- lapply(names(groups), function(i) {
  weighted_networks$lr_sig %>% filter(from %in% groups[[i]] & to %in% unique(lr_network_top_df$target)) %>% 
    group_by(from) %>% mutate(receptor_rank = dense_rank(desc(weight))) %>%
    group_by(to) %>%  summarise(summed_weight=sum(weight),  avg_weight=mean(weight),  avg_rank = mean(receptor_rank)) %>%
    mutate(type=i, color = receptor_colors[i]) %>% rename(receptor = to)
}) %>% do.call(rbind, .)

# Check distinct receptors per group
lapply(paste0("group", 1:5), function(group_name) setdiff(receptors %>% filter(type == group_name) %>% pull(receptor), receptors %>% filter(type != group_name) %>% pull(receptor)))

# Assign receptor to a specific group
receptors_filtered <- receptors %>% group_by(receptor) %>% filter(avg_rank == min(avg_rank)) %>%
  filter(summed_weight == max(summed_weight)) %>% ungroup()

Next, run the same functions as the plots before, after filtering the ligand assignment dataframe and the ligand-target links dataframe. We also add the group information to the targets

circos_links_subset <- get_ligand_target_links_oi(
  ligand_type_indication_df %>% filter(ligand %in% unlist(groups)),
  active_ligand_target_links_df %>% filter(ligand %in% unlist(groups)) %>%
    # Add group information
    inner_join(targets_filtered, by="target") %>%
    select(-target_type, n, summed_weight, avg_rank) %>%
    rename(target_type=type),
  cutoff = 0.66)

ligand_colors <- c("General" = "lawngreen",
                   "CAF-specific" = "royalblue",
                   "Endothelial-specific" = "gold")
target_colors <- rep("tomato", 4) %>% setNames(names(groups))

# Define specific gaps
vis_circos_obj_subset <- prepare_circos_visualization(circos_links_subset,
                             ligand_colors = ligand_colors,
                             target_colors = target_colors,
                             widths = list(width_same_cell_same_ligand_type = 0.6,
                                           width_different_cell = 4.5,
                                           width_ligand_target = 12,
                                           width_same_cell_same_target_type = 0.6))

Finally, we create the plot by adding an extra layer in preAllocateTracks, and we add a for loop at the end to draw the outer layer.

circos.par(gap.degree = vis_circos_obj_subset$gaps)

chordDiagram(vis_circos_obj_subset$links_circle,
             order=vis_circos_obj_subset$order,
             transparency=0,
             directional = 1,
             link.sort = TRUE,
             link.decreasing = FALSE,
             grid.col = vis_circos_obj_subset$ligand_colors,
             diffHeight = 0.005,
             direction.type = c("diffHeight", "arrows"),
             link.arr.type = "big.arrow",
             link.visible = vis_circos_obj_subset$links_circle$weight >=
               attr(vis_circos_obj_subset$links_circle, "cutoff_include_all_ligands"),
             annotationTrack = "grid",
    # Add extra track for outer layer
    preAllocateTracks = list(list(track.height = 0.025),
                             list(track.height = 0.25)))

# we go back to the first track and customize sector labels
circos.track(track.index = 2, panel.fun = function(x, y) {
    circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
        facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 0.8)
}, bg.border = NA) #


padding <- 0.05        # Gaps between receptor arcs
rou_adjustment <- 0.08 # Might need adjustment
for (group in unique(targets_filtered$type)){
  # Subset target and receptor
  targets_subset <- targets_filtered %>% filter(type == group)
  receptors_subset <- receptors_filtered %>% filter(type == group)

  # Get approximate position of outer ring
  pos <- circlize(c(0, 1), c(0, 1), track.index = 1)
  rou1 <- pos[1, "rou"]-rou_adjustment
  rou2 <- pos[2, "rou"]-rou_adjustment

  # Get range of angles from first to last target of the current group
  theta1 <- circlize:::get.sector.data((targets_subset %>% pull(target) %>% .[1]))["start.degree"]
  theta2 <- circlize:::get.sector.data((targets_subset %>% pull(target) %>% .[length(.)]))["end.degree"]

  # Scale the arc lengths according to the summed ligand-receptor weights
  receptors_subset_scaled <- receptors_subset %>%
    mutate(scaled_weight = summed_weight/sum(summed_weight)*(theta1-theta2))
  # For each receptor
  current_theta <- theta1
  for (i in 1:nrow(receptors_subset_scaled)){
    # Get end angle of the arc
    end_theta <- current_theta-(receptors_subset_scaled %>% slice(i) %>% pull(scaled_weight))
    d1 <- abs(end_theta-current_theta) # For gaps

    # Main function - we draw the arc here
    draw.sector(current_theta+(d1*-padding), end_theta-(d1*-padding),
                rou1, rou2,
                col = receptors_subset_scaled %>% slice(i) %>% pull(color),
                clock.wise = TRUE, border=NA)

    # Add text containing receptor name
    pos_text <- reverse.circlize((current_theta + end_theta)/2 + ifelse(current_theta < end_theta, 180, 0), (rou1 +  rou2)/2)
    # It's going to give a Note that point is out of plotting region
    suppressMessages(circos.text(pos_text[1,1], pos_text[1,2]+convert_y(7, "mm"),
                labels = receptors_subset_scaled %>% slice(i) %>% pull(receptor),
                niceFacing=TRUE, facing="clockwise", cex=0.6))

    current_theta <- end_theta
  }

}

circos.clear()

As mentioned previously, the resulting plot will require some further manual tweaking (e.g., ITGA5 should be in group 3).

Using the highlight.sector function

With this function, it is not possible to draw receptor arcs that end at different points from the target arcs. To demonstrate this, we will randomly assign the target genes into one of three groups (Receptors A, B, and C).

set.seed(10)
target_gene_groups <- data.frame(target_type = sample(c("Receptor A", "Receptor B", "Receptor C"),
                                                      length(unique(circos_links$target)), replace = TRUE),
                                 target = unique(circos_links$target))
head(target_gene_groups)

Again, define the colors of the receptors (outer layer) and targets (inner layer). Then, run prepare_circos_visualization to get the order of the targets and the gaps between them.

receptor_group_colors <- c("#387D7A", "#9DA9A0", "#704E2E") %>%
  setNames(unique(target_gene_groups$target_type))

target_colors <- rep("tomato", 3) %>% setNames(unique(target_gene_groups$target_type))

vis_circos_obj2 <- prepare_circos_visualization(circos_links %>% select(-target_type) %>%
                                                  # Combine target gene group information to circos_links
                                                  inner_join(target_gene_groups, by="target"),
                             ligand_colors = ligand_colors,
                             target_colors = target_colors,
                             widths = list(width_same_cell_same_ligand_type = 0.6,
                                           width_different_cell = 4.5,
                                           width_ligand_target = 12,
                                           width_same_cell_same_target_type = 0.6))

The general idea here is similar - adding an extra layer in preAllocateTracks and a for loop at the end to draw the outer layer - but the function allows for much cleaner code.

circos.par(gap.degree = vis_circos_obj2$gaps)
chordDiagram(vis_circos_obj2$links_circle,
             directional = 1,
             transparency=0,
             order=vis_circos_obj2$order,
             link.sort = TRUE,
             link.decreasing = FALSE,
             grid.col = vis_circos_obj2$ligand_colors,
             diffHeight = 0.005,
             direction.type = c("diffHeight", "arrows"),
             link.arr.type = "big.arrow",
             link.visible = vis_circos_obj2$links_circle$weight >=
               attr(vis_circos_obj2$links_circle, "cutoff_include_all_ligands"),
             annotationTrack = "grid",
    # Add extra track for outer layer
    preAllocateTracks = list(list(track.height = 0.025),
                             list(track.height = 0.2)))

# we go back to the first track and customize sector labels
circos.track(track.index = 2, panel.fun = function(x, y) {
    circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
        facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 0.8)
}, bg.border = NA) #

# Add outer layer
for (target_gene_group in unique(target_gene_groups$target_type)){
  highlight.sector(target_gene_groups %>% filter(target_type == target_gene_group) %>% pull(target),
                   track.index = 1,
                 col = receptor_group_colors[target_gene_group],
                 text = target_gene_group,
                 cex = 0.8, facing="bending.inside", niceFacing = TRUE, text.vjust = "5mm")
}

circos.clear()

References

Bonnardel et al., 2019, Immunity 51, 1–17, Stellate Cells, Hepatocytes, and Endothelial Cells Imprint the Kupffer Cell Identity on Monocytes Colonizing the Liver Macrophage Niche

saeyslab/nichenetr documentation built on April 27, 2024, 9:24 p.m.