【视频教程】MultiNicheNet（1）：多组单细胞通讯差异分析-基于基本流程

# install.packages("devtools")# devtools::install_github("saeyslab/nichenetr")devtools::install_github("saeyslab/multinichenetr")library(multinichenetr)

#加载数据库library(SingleCellExperiment)library(dplyr)library(ggplot2)library(nichenetr)library(multinichenetr)library(Seurat)
#lrlr_network_all = readRDS('./multinechenetr-human-network/lr_network_human_allInfo_30112033.rds') %>%   mutate(ligand = convert_alias_to_symbols(ligand, organism = "human"),          receptor = convert_alias_to_symbols(receptor, organism = "human"))
lr_network_all = lr_network_all  %>% mutate(ligand = make.names(ligand), receptor = make.names(receptor)) #对lr_network_all数据框中的ligand列和receptor列的值进行处理，使得这些值变为有效的R变量名。lr_network = lr_network_all %>% distinct(ligand, receptor)#去除重复的行，只保留唯一的组合
#weight matrixligand_target_matrix = readRDS('./multinechenetr-human-network/ligand_target_matrix_nsga2r_final.rds')# target genes in rows, ligands in columnscolnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = "human") %>% make.names()rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = "human") %>% make.names()
lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix))ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()]

#load datasce_subset_misc <- readRDS("D:/KS项目/公众号文章/Multinechenetr多组受配体分析/sce_subset_misc.rds")sce = alias_to_symbol_SCE(sce_subset_misc, "human") %>% makenames_SCE()# sce1 <- as.Seurat(sce_subset_misc)# unique(sce1$Annotation_v2.0)# unique(sce1$MIS.C.AgeTier)# unique(sce1$ShortID)
#！！！multinichenetr建议每个样本中每种celltype有足够的数量table(sce1$ShortID, sce1$MIS.C.AgeTier)#      A    M    S# A1 1335    0    0# A2  863    0    0# A3  816    0    0# A4 1312    0    0# M1    0 1052    0# M2    0  229    0# M3    0  704    0# M4    0  927    0# M5    0  248    0# M6    0  368    0# M7    0 1756    0# S1    0    0  812# S2    0    0  497# S3    0    0  731# S4    0    0 1174# S5    0    0  519

table(sce1$ShortID, sce1$Annotation_v2.0)#           L_NK_CD56._CD16.   L_T_TIM3._CD38._HLADR.   M_Monocyte_CD16# A1             1302                     27               6# A2              712                    137              14# A3              462                    291              63# A4             1273                     36               3# M1              158                    873              21# M2               53                    102              74# M3              273                    426               5# M4              573                    130             224# M5               96                     56              96# M6               47                    320               1# M7              367                   1193             196# S1              620                     76             116# S2              355                     38             104# S3              352                     50             329# S4             1082                     21              71# S5              329                     23             167

# sample_id = "ShortID"# group_id = "MIS.C.AgeTier"# celltype_id = "Annotation_v2.0"# rm(sce1,sce_subset_misc)

#设置分组、celltype注释、样本idsample_id = "ShortID"group_id = "MIS.C.AgeTier"# celltype_id = "Annotation_v2.0"#有协变量或者批次效应的数据，可设置，本数据不包含，设置NA即可covariates = NAbatches = NA
#SummarizedExperiment::colData(sce)$ShortID = SummarizedExperiment::colData(sce)$ShortID %>% make.names()#SummarizedExperiment::colData(sce)$MIS.C.AgeTier = SummarizedExperiment::colData(sce)$MIS.C.AgeTier %>% make.names()#SummarizedExperiment::colData(sce)$Annotation_v2.0 = SummarizedExperiment::colData(sce)$Annotation_v2.0 %>% make.names()

#设置比较组contrasts_oi = c("'M-(S+A)/2','S-(M+A)/2','A-(S+M)/2'")contrast_tbl = tibble(contrast = c("M-(S+A)/2","S-(M+A)/2", "A-(S+M)/2"),                      group = c("M","S","A"))

#设置感兴趣的sender、receiver cellslibrary(dplyr)sce$celltype <- ""sce$celltype <- case_when(sce$Annotation_v2.0=="L_T_TIM3._CD38._HLADR." ~ "HLADR",                          sce$Annotation_v2.0=="L_NK_CD56._CD16." ~ "NK",                          sce$Annotation_v2.0=="M_Monocyte_CD16" ~ "Monocyte")
celltype_id = "celltype"SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()# [1] "HLADR"    "NK"       "Monocyte"
# senders_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()# receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()senders_oi = c("HLADR","NK","Monocyte")receivers_oi = c("HLADR","NK","Monocyte")sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% c(senders_oi, receivers_oi)]

#过滤细胞abundance_info = get_abundance_info(  sce = sce,  sample_id = sample_id,   group_id = group_id,   celltype_id = celltype_id,  min_cells = 10,  senders_oi = senders_oi,   receivers_oi = receivers_oi,  batches = batches)

#过滤基因min_cells=10fraction_cutoff = 0.05min_sample_prop = 0.5
frq_list = get_frac_exprs(  sce = sce,   sample_id = sample_id, celltype_id =  celltype_id, group_id = group_id,   batches = batches,   min_cells = min_cells,   fraction_cutoff = fraction_cutoff,   min_sample_prop = min_sample_prop)#保留至少在一种celltype中表达的基因genes_oi = frq_list$expressed_df %>% filter(expressed == TRUE) %>% pull(gene) %>% unique() sce = sce[genes_oi, ]

#伪bulk表达量计算abundance_expression_info = process_abundance_expression_info(  sce = sce,   sample_id = sample_id, group_id = group_id, celltype_id = celltype_id,   min_cells = min_cells,   senders_oi = senders_oi,   receivers_oi = receivers_oi,   lr_network = lr_network,   batches = batches,   frq_list = frq_list,   abundance_info = abundance_info)

#差异分析#差异分析DE_info = get_DE_info(  sce = sce,   sample_id = sample_id, group_id = group_id, celltype_id = celltype_id,   batches = batches, covariates = covariates,   contrasts_oi = contrasts_oi,   min_cells = min_cells,   expressed_df = frq_list$expressed_df)#认为在每种细胞类型中表达的基因
celltype_de = DE_info$celltype_de$de_output_tidy
#将差异结果于ligands和receptor结合，看看他们表达的差异性sender_receiver_de = multinichenetr::combine_sender_receiver_de(  sender_de = celltype_de,  receiver_de = celltype_de,  senders_oi = senders_oi,  receivers_oi = receivers_oi,  lr_network = lr_network)

#ligand activity计算、配体靶基因预测logFC_threshold = 0.50p_val_threshold = 0.05p_val_adj = FALSE geneset_assessment = contrast_tbl$contrast %>%   lapply(process_geneset_data,          celltype_de, logFC_threshold, p_val_adj, p_val_threshold) %>% bind_rows() geneset_assessment

ligand_activities_targets_DEgenes = suppressMessages(suppressWarnings(  get_ligand_activities_targets_DEgenes(    receiver_de = celltype_de,    receivers_oi = intersect(receivers_oi, celltype_de$cluster_id %>% unique()),    ligand_target_matrix = ligand_target_matrix,    logFC_threshold = logFC_threshold,    p_val_threshold = p_val_threshold,    p_val_adj = F,    top_n_target = 250,    verbose = T,     n.cores = 2  )))

#受配体对优先级排序，得到最终结果sender_receiver_tbl = sender_receiver_de %>% distinct(sender, receiver)metadata_combined = SummarizedExperiment::colData(sce) %>% tibble::as_tibble()
if(!is.na(batches)){  grouping_tbl = metadata_combined[,c(sample_id, group_id, batches)] %>%     tibble::as_tibble() %>% distinct()  colnames(grouping_tbl) = c("sample","group",batches)} else {  grouping_tbl = metadata_combined[,c(sample_id, group_id)] %>% tibble::as_tibble() %>% distinct()  colnames(grouping_tbl) = c("sample","group")}

prioritization_tables = suppressMessages(multinichenetr::generate_prioritization_tables(  sender_receiver_info = abundance_expression_info$sender_receiver_info,  sender_receiver_de = sender_receiver_de,  ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes,  contrast_tbl = contrast_tbl,  sender_receiver_tbl = sender_receiver_tbl,  grouping_tbl = grouping_tbl,  scenario = "regular", # all prioritization criteria will be weighted equally  fraction_cutoff = 0.05,   abundance_data_receiver = abundance_expression_info$abundance_data_receiver,  abundance_data_sender = abundance_expression_info$abundance_data_sender,  ligand_activity_down = F))

#保存结果：path = "./"
multinichenet_output = list(  celltype_info = abundance_expression_info$celltype_info,  celltype_de = celltype_de,  sender_receiver_info = abundance_expression_info$sender_receiver_info,  sender_receiver_de =  sender_receiver_de,  ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes,  prioritization_tables = prioritization_tables,  grouping_tbl = grouping_tbl) 

saveRDS(multinichenet_output, paste0(path, "multinichenet_output.rds"))

#弦图#可视化每个组中，感兴趣sender与receiver最优先级的受配体对#首先提取数据，可视化排序score top多少的受配体对prioritized_tbl_oi_all = get_top_n_lr_pairs(  multinichenet_output$prioritization_tables,   top_n = 30,   rank_per_group = T)

prioritized_tbl_oi =   multinichenet_output$prioritization_tables$group_prioritization_tbl %>%  filter(id %in% prioritized_tbl_oi_all$id) %>%  distinct(id, sender, receiver, ligand, receptor, group) %>%   left_join(prioritized_tbl_oi_all)
# prioritized_tbl_oi <- na.omit(prioritized_tbl_oi)prioritized_tbl_oi$prioritization_score[is.na(prioritized_tbl_oi$prioritization_score)] = 0

senders_receivers = union(prioritized_tbl_oi$sender %>% unique(), prioritized_tbl_oi$receiver %>% unique()) %>% sort()
colors_sender = RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral') %>% magrittr::set_names(senders_receivers)colors_receiver = RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral') %>% magrittr::set_names(senders_receivers)

par(mfrow = c(2,2), xpd=TRUE)circos_list = make_circos_group_comparison(prioritized_tbl_oi, colors_sender, colors_receiver)

#整合气泡图#可视化每组受配体详细信息prioritized_tbl_oi_M_30 = prioritized_tbl_oi_all %>%  filter(group == "M")plot_oi = make_sample_lr_prod_activity_plots_Omnipath(  multinichenet_output$prioritization_tables,   prioritized_tbl_oi_M_30 %>% inner_join(lr_network_all))plot_oi

#可视化target基因#可以设定感兴趣的组，感兴趣的receiver cell
group_oi = "M"receiver_oi = "Monocyte"prioritized_tbl_oi_M_10 = get_top_n_lr_pairs(  multinichenet_output$prioritization_tables,   10,  groups_oi = group_oi, receivers_oi = receiver_oi)

combined_plot = make_ligand_activity_target_plot(  group_oi,   receiver_oi,   prioritized_tbl_oi_M_10,  multinichenet_output$prioritization_tables,   multinichenet_output$ligand_activities_targets_DEgenes,   contrast_tbl,   multinichenet_output$grouping_tbl,   multinichenet_output$celltype_info,   ligand_target_matrix,   plot_legend = FALSE)combined_plot