NicheNet analysis

#Nichenet ligand-receptor interactions
#This section covers Figure 6D and single cell RNA sequencing data deposited in GSE152042

#Load packages

library(Seurat)

library(SeuratData)

library(SeuratWrappers)

library(cowplot)

library(patchwork)

library(ggplot2)

library(nichenetr)

library(tidyverse)

#Prepare data - "epithelial" and "Seurat" objects correspond to the resclustering analyses performed in https://doi.org/10.7554/eLife.62810 

epithelial <- RenameIdents(epithelial, '0' = "Basal 1", '1' = "Basal 2", '2' = "Inflammatory KC", '3' = "Inflammatory KC 2", '4' = "Differentiating KC", '5' = "Proliferating KC", '6' = "Inflammatory KC 3", '7' = "Inflammatory KC 4", '8' = "Proliferating KC 2", '9' = "Basal 3")

Stromal <- RenameIdents(Stromal, '0' = "Fb1", '1' = "MyoFb", '2' = "Fb2", '3' = "Pericyte", '4' = "Fb3", '5' = "Fb4", '6' = "Fb5")

DefaultAssay(Stromal) <- "RNA"

DefaultAssay(epithelial) <- "RNA"

#Rename conditions to "Healthy" and "Periodontitis"

epithelial@meta.data[['new_ids']] <- apply(epithelial@meta.data, 1, function(x) {ifelse(x[['stim']] == 'Healthy', 'Healthy', 'Periodontitis')})

Stromal@meta.data[['new_ids']] <- apply(Stromal@meta.data, 1, function(x) {ifelse(x[['stim']] == 'Healthy', 'Healthy', 'Periodontitis')})

#Subset healthy stromal subtypes

healthy_stromal <- subset(Stromal, idents = "Healthy")

DimPlot(healthy_stromal, reduction = "umap")

Idents(healthy_stromal) <- "seurat_clusters"

healthy_fibroblasts <- subset(healthy_stromal, idents = c("0", "2", "4", "5", "6"))

RenameIdents(healthy_fibroblasts, '0' = "Fb1", '2' = "Fb2", '4' = "Fb3", '5' = "Fb4", '6' = "Fb5")

#Subset basal epithelial cells

Basal1 <- subset(epithelial, idents = "Basal 1")

#The code below follows vignettes described in https://github.com/saeyslab/nichenetr: vignette("seurat_wrapper", package="nichenetr"); vignette("seurat_steps", package="nichenetr")

#Read in the NicheNet ligand-receptor prior network and ligand-target matrix
ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))
ligand_target_matrix[1:5,1:5]

lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
head(lr_network)

weighted_networks = readRDS(url("https://zenodo.org/record/3260758/files/weighted_networks.rds"))
weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c("from","to"))

head(weighted_networks$lr_sig)
head(weighted_networks$gr)

#read expression data of interacting cells

##seurat objects

DefaultAssay(healthy_fibroblasts) <- "RNA"

DefaultAssay(Basal1) <- "RNA"

#Define a "sender/niche" cell population and a "receiver/target" cell population and determine which genes are expressed in both populations. In this study, the receiver cell population is the epithelial basal cells (Basal 1), whereas the sender populations are the stromal subpopulations. We hypothesised that Fb2 is a niche population.

receiver = "Basal 1"

expressed_genes_receiver = get_expressed_genes(receiver, epithelial, pct = 0.10)

background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]

sender = "Fb2"

expressed_genes_sender = get_expressed_genes(sender, healthy_fibroblasts, pct = 0.10)

list_expressed_genes_sender = sender %>% unique() %>% lapply(get_expressed_genes, healthy_fibroblasts, 0.10)

expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()

#Define a gene set of interest genes in the “receiver/target” cell population that are potentially affected by ligands expressed by interacting cells (e.g. genes differentially expressed upon cell-cell interaction)

condition_oi = "Periodontitis"

condition_reference = "Healthy"

DE_table_receiver = FindMarkers(object = seurat_obj_receiver, ident.1 = condition_oi, ident.2 = condition_reference, min.pct = 0.10) %>% rownames_to_column("gene")

geneset_oi = DE_table_receiver %>% filter(p_val_adj <= 0.05 & abs(avg_logFC) >= 0.25) %>% pull(gene)

geneset_oi = geneset_oi %>% .[. %in% rownames(ligand_target_matrix)]

#define a set of potential ligands which are expressed by the sender niche

ligands = lr_network %>% pull(from) %>% unique()

receptors = lr_network %>% pull(to) %>% unique()

expressed_ligands = intersect(ligands,expressed_genes_sender)

expressed_receptors = intersect(receptors,expressed_genes_receiver)

potential_ligands = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>% pull(from) %>% unique()

#perform NicheNet ligand activity analysis: rank the potential ligands

ligand_activities = predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)

ligand_activities = ligand_activities %>% arrange(-pearson) %>% mutate(rank = rank(desc(pearson)))

ligand_activities

best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand) %>% unique()

#check which sender population expresses these top-ranked ligands

DotPlot(healthy_fibroblasts, features = best_upstream_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()