偷偷问一下,关注了吗?
内容获取
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 134
table(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---------------------------------------------------------------------------------------
#lr
lr_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 matrix
ligand_target_matrix = readRDS('D:/KS项目/公众号文章/Multinechenetr多组受配体分析/multinechenetr-human-network/ligand_target_matrix_nsga2r_final.rds')# target genes in rows, ligands in columns
colnames(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 = 10
abundance_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 = 1
fraction_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-receivers
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)
#step5:Ligand activity prediction--------------------------------------------------------------------------------------
#inspect the geneset_oi-vs-background ratios
logFC_threshold = 0.3 # lower here for FindMarkers than for Pseudobulk-EdgeR
p_val_threshold = 0.05
p_val_adj = TRUE
geneset_assessment = contrast_tbl$contrast %>%
celltype_de, logFC_threshold, p_val_adj, p_val_threshold) %>%
#Perform the ligand activity analysis and ligand-target inference
top_n_target = 250
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 = p_val_adj,
top_n_target = top_n_target,
verbose = F,
2 =
)
))
#step6:Prioritization--------------------------------------------------------------------------------------
ligand_activity_down = FALSE
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(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)] = 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