偷偷问一下，关注了吗？
内容获取
1、购买打包合集（《KS科研分享与服务》付费内容打包集合），价格感人，可以加入微信VIP群（答疑交流群，甚至有小伙伴觉得群比代码更好），可以获取建号以来所有内容，群成员专享视频教程，提前更新，其他更多福利！
2、《KS科研分享与服务》公众号有QQ群，进入门槛是20元（完全是为了防止白嫖党，请理解），请考虑清楚。群里有免费推文的注释代码和示例数据（终身拥有），没有付费内容，群成员福利是购买单个付费内容半价！
需要者详情请联系作者（非需要者勿扰，处理太费时间）：
https://www.biorxiv.org/content/10.1101/2023.06.13.544751v1
一、准备工作
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 = NAbatches = 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.05frq_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 = TRUEgeneset_assessment = contrast_tbl$contrast %>%celltype_de, logFC_threshold, p_val_adj, p_val_threshold) %>%#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,= 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()= 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()= 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()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)] = 0senders_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
觉得我们分享有些用的，点个赞再走呗！
关注我们获取精彩内容：
关注不迷路：扫描下面二维码关注公众号！
B站视频号链接：https://space.bilibili.com/471040659?spm_id_from=333.1007.0.0
关注 KS科研分享与服务，
认清正版优质内容和服务！
优质内容持续输出，物超所值！
合作联系：ks_account@163.com
