NicheNET

官网教程：https://github.com/saeyslab/nichenetr

从 Seurat 对象开始执行 NicheNet 分析–逐步分析:https://github.com/saeyslab/nichenetr/blob/master/vignettes/seurat_steps.md

背景数据获取：https://zenodo.org/record/3260758#.ZCVH4XZBxPY

原文中的坑：geneset_oi就是一个基因列表，字符型

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

数据准备
data <-readRDS(“s.adult.rds”)
meta = data.frame(data@meta.data)

#以下三个文件在https://zenodo.org/record/3260758#.ZCVH4XZBxPY获取

ligand_target_matrix = readRDS(“ligand_target_matrix.rds”)

lr_network = readRDS(“lr_network.rds”)

weighted_networks =readRDS(“weighted_networks.rds”)
weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c(“from”,”to”))

开始分析
seuratObj = data

seuratObj@meta.data$group %>% table()

#设置分组类型，可为celltype

Idents(seuratObj) <- “group”

receiver
receiver = “ATII-3”
expressed_genes_receiver = get_expressed_genes(receiver, seuratObj, pct = 0.10)
background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]

sender
sender_celltypes = c(“Aberrant_Basaloid”,”ATI”, “Basal”, “Ciliated”, “Club”, “Goblet”,”Mesothelial”,”PNEC”,”Ionocyte”)
list_expressed_genes_sender = sender_celltypes %>% unique() %>% lapply(get_expressed_genes, seuratObj, 0.10) # lapply to get the expressed genes of every sender cell type separately here
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()

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()

#原文中的坑geneset_oi：就是一个基因列表

markers_sig <- read.delim(“E:\Ipf\1-data process\markers_sig(2_0.01).txt”, check.names = F)
markers_sig <-markers_sig[which(markers_sig$cluster %in% “18”),]

geneset_oi<-markers_sig$gene 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

ligand activity DotPlot（左2图）
best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand) %>% unique()
DotPlot(seuratObj, features = best_upstream_ligands %>% rev(), cols = “RdYlBu”) + RotatedAxis()+coord_flip()

ligand_pearson heatmap（左1图）
ligand_pearson_matrix = ligand_activities %>% select(pearson) %>% as.matrix() %>% magrittr::set_rownames(ligand_activities$test_ligand)
rownames(ligand_pearson_matrix) = rownames(ligand_pearson_matrix) %>% make.names()
colnames(ligand_pearson_matrix) = colnames(ligand_pearson_matrix) %>% make.names()
vis_ligand_pearson = ligand_pearson_matrix[order_ligands, ] %>% as.matrix(ncol = 1) %>% magrittr::set_colnames(“Pearson”)
p_ligand_pearson = vis_ligand_pearson %>% make_heatmap_ggplot(“Prioritized ligands”,”Ligand activity”, color = “darkorange”,legend_position = “top”, x_axis_position = “top”, legend_title = “Pearson correlation coefficient\ntarget gene prediction ability)”) + theme(legend.text = element_text(size = 9))
p_ligand_pearson

ligand activity heatmap（右1图）
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 200) %>% bind_rows() %>% drop_na()
active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.33)
order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev() %>% make.names()
order_targets = active_ligand_target_links_df$target %>% unique() %>% intersect(rownames(active_ligand_target_links)) %>% make.names()
rownames(active_ligand_target_links) = rownames(active_ligand_target_links) %>% make.names() # make.names() for heatmap visualization of genes like H2-T23
colnames(active_ligand_target_links) = colnames(active_ligand_target_links) %>% make.names() # make.names() for heatmap visualization of genes like H2-T23

vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()
p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot(“Prioritized ligands”,”Predicted target genes”, color = “purple”,legend_position = “top”, x_axis_position = “top”,legend_title = “Regulatory potential”) + theme(axis.text.x = element_text(face = “italic”)) + scale_fill_gradient2(low = “whitesmoke”, high = “purple”, breaks = c(0,0.0045,0.0090))
p_ligand_target_network