【视频演示】Nichenetr V2（三）：nichenet分析之基于ligand-receptor的ligands排序

#load packages and datalibrary(nichenetr) library(Seurat)library(SeuratObject)library(tidyverse)

#load datasetwd("D:\\KS项目\\公众号文章\\nichenet单细胞通讯分析及可视化")load("D:/KS项目/公众号文章/nichenet单细胞通讯分析及可视化/scRNA_nichenet.RData")lr_network <- readRDS('./nichent-human-network/lr_network_human_21122021.rds')ligand_target_matrix <- readRDS('./nichent-human-network/ligand_target_matrix_nsga2r_final.rds')# target genes in rows, ligands in columnsweighted_networks <- readRDS('./nichent-human-network/weighted_networks_nsga2r_final.rds')lr_network <- lr_network %>% distinct(from, to)

#Perform the NicheNet analysis## 1. Define set of potential ligandsreceiver = "Tcell"expressed_genes_receiver <- get_expressed_genes(receiver, sce, pct = 0.1)
sender_celltypes <- c("Endo","Fib","Musc","Ker","APCs","Tcell","Mast","LY","Mela")list_expressed_genes_sender <- sender_celltypes %>% unique() %>% lapply(get_expressed_genes, sce, 0.1)expressed_genes_sender <- list_expressed_genes_sender %>% unlist() %>% unique()

all_ligands <- unique(lr_network$from)all_receptors <- unique(lr_network$to)
expressed_ligands <- intersect(all_ligands, expressed_genes_sender)expressed_receptors <- intersect(all_receptors, expressed_genes_receiver)
potential_ligands <- lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>%  pull(from) %>% unique()

# 2. Define the gene set of interestseurat_obj_receiver <- subset(sce, idents = receiver)DE_table_receiver <-  FindMarkers(object = seurat_obj_receiver,                                  ident.1 = "SD",                                   ident.2 = "HC",                                  group.by = "group",                                  min.pct = 0.1) %>% rownames_to_column("gene")
geneset_oi <- DE_table_receiver %>% filter(p_val_adj <= 0.05 & abs(avg_log2FC) >= 0.25) %>% pull(gene)geneset_oi <- geneset_oi %>% .[. %in% rownames(ligand_target_matrix)]
# 3. Define background genesbackground_expressed_genes <- expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
# 4. Perform NicheNet ligand activity analysisligand_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(-aupr_corrected) %>% mutate(rank = rank(desc(aupr_corrected)))

lr_network_filtered <-  lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors)# 只计算实验组, with genes that are in the ligand-receptor networkDE_table <- FindAllMarkers(subset(sce, subset = group == "SD"),                           min.pct = 0, logfc.threshold = 0, return.thresh = 1,                           features = unique(unlist(lr_network_filtered))) 

# Average expression information - only for SD conditionexpression_info <- get_exprs_avg(sce, "celltype", condition_colname = "group",                                  condition_oi = "SD",                                 features = unique(unlist(lr_network_filtered)))
# Calculate condition specificity - only for datasets with two conditions!condition_markers <- FindMarkers(object = sce, ident.1 = "SD", ident.2 = "HC",                                 group.by = "group", min.pct = 0, logfc.threshold = 0,                                 features = unique(unlist(lr_network_filtered))) %>% rownames_to_column("gene")

# Combine DE of senders and receivers -> used for prioritizationprocessed_DE_table <- process_table_to_ic(DE_table, table_type = "celltype_DE", lr_network_filtered,                                          senders_oi = sender_celltypes, receivers_oi = receiver)
processed_expr_table <- process_table_to_ic(expression_info, table_type = "expression", lr_network_filtered)
processed_condition_markers <- process_table_to_ic(condition_markers, table_type = "group_DE", lr_network_filtered)

# 不跑，一步法包装函数# info_tables <- generate_info_tables(seuratObj,#                                     celltype_colname = "celltype",#                                     senders_oi = sender_celltypes,#                                     receivers_oi = receiver,#                                     lr_network = lr_network_filtered,#                                     condition_colname = "aggregate",#                                     condition_oi = condition_oi,#                                     condition_reference = condition_reference,#                                     scenario = "case_control")

prioritizing_weights = c("de_ligand" = 1,                         "de_receptor" = 1,                         "activity_scaled" = 1,                         "exprs_ligand" = 1,                         "exprs_receptor" = 1,                         "ligand_condition_specificity" = 1,                         "receptor_condition_specificity" = 1)
prior_table <- generate_prioritization_tables(processed_expr_table,                                              processed_DE_table,                                              ligand_activities,                                              processed_condition_markers,                                              prioritizing_weights)
#选取用去看排序的列即可prior_table <- prior_table %>% select(c('sender', 'receiver', 'ligand', 'receptor', 'scaled_p_val_adapted_ligand', 'scaled_p_val_adapted_receptor', 'scaled_avg_exprs_ligand', 'scaled_avg_exprs_receptor', 'scaled_p_val_adapted_ligand_group', 'scaled_p_val_adapted_receptor_group', 'scaled_activity'))