【视频教程】MultiNicheNet(2):适用于无重复样本的受配体差异分析流程演示

学术   2024-11-01 09:00   重庆  

偷偷问一下,关注了吗

内容获取


1、购买打包合集(《KS科研分享与服务》付费内容打包集合),价格感人,可以加入微信VIP群(答疑交流群,甚至有小伙伴觉得群比代码更好),可以获取建号以来所有内容,群成员专享视频教程,提前更新,其他更多福利!


2、《KS科研分享与服务》公众号有QQ群,进入门槛是20元(完全是为了防止白嫖党,请理解),请考虑清楚。群里有免费推文的注释代码和示例数据(终身拥有),没有付费内容,群成员福利是购买单个付费内容半价!


需要者详情请联系作者(非需要者勿扰,处理太费时间):


关于MultiNicheNet的分析教程,官网有及其详细的流程,以及各种情况的教程,只要肯花时间慢慢阅读,都没有问题,除了我们演示的情况,官网还有具有协变量数据、批次效应数据、添加蛋白组学数据等等教程,流程上大差不差。MultiNicheNet分析过程并不难,甚至可以说基本上不会遇见报错。我们之前演示了基本流程,适用于有重复样本(【视频教程】MultiNicheNet(1):多组单细胞通讯差异分析-基于基本流程)。这里我们演示一下可能大多数小伙伴数据情况,每个组并没有重复数据或者即使有重复,也是3个样本以下。总体分析流程和普通的没有任何区别,只是在某些步骤函数不太一样。因为没有重复,所以差异基因分析不能使用伪bulk,使用的是Findmarkers。

目前文章是在预印本:

https://www.biorxiv.org/content/10.1101/2023.06.13.544751v1

官网有详细的教程:
https://github.com/saeyslab/multinichenetr?tab=readme-ov-file

MultiNicheNet的分析分为三大块:”准备工作“,”分析工作“,”结果可视化‘,这样就会感觉很清晰:接下来看看过程:

一、准备工作

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

觉得我们分享有些用的,点个赞再走呗!

关注我们获取精彩内容:


关注不迷路:扫描下面二维码关注公众号!
B站视频号链接https://space.bilibili.com/471040659?spm_id_from=333.1007.0.0




关注 KS科研分享与服务,

认清正版优质内容和服务!

优质内容持续输出,物超所值!

合作联系:ks_account@163.com

新的板块-重要通知-双向选择

KS科研分享与服务
科研学习交流于分享,生信学习笔记,科研经历和生活!
 最新文章