|
关于MultiNicheNet的分析教程，官网有及其详细的流程，以及各种情况的教程，只要肯花时间慢慢阅读，都没有问题，除了我们演示的情况，官网还有具有协变量数据、批次效应数据、添加蛋白组学数据等等教程，流程上大差不差。MultiNicheNet分析过程并不难，甚至可以说基本上不会遇见报错。我们之前演示了基本流程，适用于有重复样本（【视频教程】MultiNicheNet（1）：多组单细胞通讯差异分析-基于基本流程）。这里我们演示一下可能大多数小伙伴数据情况，每个组并没有重复数据或者即使有重复，也是3个样本以下。总体分析流程和普通的没有任何区别，只是在某些步骤函数不太一样。因为没有重复，所以差异基因分析不能使用伪bulk，使用的是Findmarkers。 https://www./content/10.1101/2023.06.13.544751v1 https://github.com/saeyslab/multinichenetr?tab=readme-ov-fileMultiNicheNet的分析分为三大块：”准备工作“，”分析工作“，”结果可视化'，这样就会感觉很清晰：接下来看看过程：library(SingleCellExperiment)library(dplyr)library(ggplot2)library(multinichenetr)library(Seurat)library(tidyverse)library(nichenetr) #load single cell data---------------------------------------------------------------------------------------load("D:/KS项目/公众号文章/nichenet单细胞通讯分析及可视化/scRNA_nichenet.RData")table(sce$celltype, sce$group)# HC SD# Endo 434 1201# Fib 878 1083# Musc 526 899# Ker 792 1810# APCs 217 825# Tcell 343 159# Mast 67 298# LY 33 134table(sce$orig.ident, sce$celltype)# Endo Fib Musc Ker APCs Tcell Mast LY Mela# HC_scRNA 434 878 526 792 217 343 67 33 47# SD_scRNA 1201 1083 899 1810 825 159 298 134 89
#convert to SCE object---------------------------------------------------------------------------------------sce = Seurat::as.SingleCellExperiment(sce, assay = "RNA")sce = alias_to_symbol_SCE(sce, "human") %>% makenames_SCE()
##set up MultiNicheNet---------------------------------------------------------------------------------------#lrlr_network_all = readRDS('D:/KS项目/公众号文章/Multinechenetr多组受配体分析/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('D:/KS项目/公众号文章/Multinechenetr多组受配体分析/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()]
#Prepare settings of the MultiNicheNet cell-cell communication analysis---------------------------------------------------------------------------------------#定义metadata中关于group, sample and cell type IDs的列sample_id = "orig.ident"group_id = "group"celltype_id = "celltype"covariates = NA batches = NA
#set sce name: make.names---------------------------------------------------------------------------------------SummarizedExperiment::colData(sce)$orig.ident = SummarizedExperiment::colData(sce)$orig.ident %>% make.names()SummarizedExperiment::colData(sce)$group = SummarizedExperiment::colData(sce)$group %>% make.names()SummarizedExperiment::colData(sce)$celltype = SummarizedExperiment::colData(sce)$celltype %>% make.names()
#set contrasts of interest，and sender、receiver cells---------------------------------------------------------------------------------------contrasts_oi = c("'SD-HC','HC-SD'")contrast_tbl = tibble(contrast = c("SD-HC","HC-SD"),group = c("SD","HC"))
senders_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique()
#step1-celltype filtering---------------------------------------------------------------------------------------min_cells = 10abundance_info = get_abundance_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, batches = batches)
#step2:Gene filtering---------------------------------------------------------------------------------------min_sample_prop = 1fraction_cutoff = 0.05
frq_list = get_frac_exprs_sampleAgnostic(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)#Now only keep genes that are expressed by at least one cell type:genes_oi = frq_list$expressed_df %>% filter(expressed == TRUE) %>% pull(gene) %>% unique() sce = sce[genes_oi, ]
#step3:Expression calculation---------------------------------------------------------------------------------------abundance_expression_info = process_abundance_expression_info( sce = sce, sample_id = group_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)
#step4:ifferential expression (DE) analysis--------------------------------------------------------------------------------------#DE: FindMarkers approach.DE_info = get_DE_info_sampleAgnostic( sce = sce, group_id = group_id, celltype_id = celltype_id, contrasts_oi = contrasts_oi, expressed_df = frq_list$expressed_df, min_cells = min_cells, contrast_tbl = contrast_tbl)
celltype_de = DE_info$celltype_de_findmarkers
#Combine DE information for ligand-senders and receptors-receiverssender_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)
#step5:Ligand activity prediction--------------------------------------------------------------------------------------#inspect the geneset_oi-vs-background ratioslogFC_threshold = 0.3 # lower here for FindMarkers than for Pseudobulk-EdgeRp_val_threshold = 0.05p_val_adj = TRUE geneset_assessment = contrast_tbl$contrast %>% lapply(process_geneset_data, celltype_de, logFC_threshold, p_val_adj, p_val_threshold) %>% bind_rows()
#Perform the ligand activity analysis and ligand-target inferencetop_n_target = 250ligand_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 = p_val_adj, top_n_target = top_n_target, verbose = F, n.cores = 2 )))
#step6:Prioritization--------------------------------------------------------------------------------------ligand_activity_down = FALSEsender_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(group_id, batches)] %>% tibble::as_tibble() %>% dplyr::distinct() colnames(grouping_tbl) = c("group",batches) grouping_tbl = grouping_tbl %>% mutate(sample = group) grouping_tbl = grouping_tbl %>% tibble::as_tibble()} else { grouping_tbl = metadata_combined[,c(group_id)] %>% tibble::as_tibble() %>% dplyr::distinct() colnames(grouping_tbl) = c("group") grouping_tbl = grouping_tbl %>% mutate(sample = group) %>% select(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 = "no_frac_LR_expr", # fraction_cutoff = fraction_cutoff, abundance_data_receiver = abundance_expression_info$abundance_data_receiver, abundance_data_sender = abundance_expression_info$abundance_data_sender, ligand_activity_down = ligand_activity_down))
#step7:Save all the output of MultiNicheNet--------------------------------------------------------------------------------------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, lr_target_prior_cor = tibble())
saveRDS(multinichenet_output, paste0(path, "multinichenet_output.rds"))
prioritized_tbl_oi_all = get_top_n_lr_pairs(prioritization_tables, top_n = 50, rank_per_group = F)prioritized_tbl_oi = 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$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)
circos_list = make_circos_group_comparison(prioritized_tbl_oi, colors_sender, colors_receiver)
#---------------------------------------------------------------------------------group_oi = "SD"prioritized_tbl_oi_M_30 = get_top_n_lr_pairs( prioritization_tables, top_n = 30, groups_oi = group_oi)plot_oi = make_sample_lr_prod_activity_plots( prioritization_tables, prioritized_tbl_oi_M_30)plot_oi
 觉得我们分享有些用的，点个赞再走呗！
|
