作者仗谆，追風(fēng)少年i
系統(tǒng)總結(jié)nichenet多條件通訊差異分析
nichenet的原理大家可以參考文章10X單細(xì)胞（10X空間轉(zhuǎn)錄組）通訊分析之NicheNet，10X單細(xì)胞（10X空間轉(zhuǎn)錄組）空間相關(guān)性分析和cellphoneDB與NicheNet聯(lián)合進(jìn)行細(xì)胞通訊分析回铛，單細(xì)胞分析之細(xì)胞交互-5：NicheNet多組間互作比較空民。
Differential NicheNet analysis between niches of interest
關(guān)于niche摩渺，在空間轉(zhuǎn)錄組上很常見偶洋，就是生態(tài)位饲趋、微環(huán)境幽邓。
The goal of Differential NicheNet is to predict ligand-receptors pairs that are both differentially expressed and active between different niches of interest.
Load in packages
library(nichenetr)
library(RColorBrewer)
library(tidyverse)
library(Seurat)
示例數(shù)據(jù)
seurat_obj = readRDS(url("https://zenodo.org/record/5840787/files/seurat_obj_subset_integrated_zonation.rds"))
DimPlot(seurat_obj, group.by = "celltype", label = TRUE)
seurat_obj = SetIdent(seurat_obj, value = "celltype")
Read in the NicheNet ligand-receptor network and ligand-target matrix
ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))
ligand_target_matrix[1:5,1:5] # target genes in rows, ligands in columns
## CXCL1 CXCL2 CXCL3 CXCL5 PPBP
## A1BG 3.534343e-04 4.041324e-04 3.729920e-04 3.080640e-04 2.628388e-04
## A1BG-AS1 1.650894e-04 1.509213e-04 1.583594e-04 1.317253e-04 1.231819e-04
## A1CF 5.787175e-04 4.596295e-04 3.895907e-04 3.293275e-04 3.211944e-04
## A2M 6.027058e-04 5.996617e-04 5.164365e-04 4.517236e-04 4.590521e-04
## A2M-AS1 8.898724e-05 8.243341e-05 7.484018e-05 4.912514e-05 5.120439e-05
lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
lr_network = lr_network %>% mutate(bonafide = ! database %in% c("ppi_prediction","ppi_prediction_go"))
lr_network = lr_network %>% dplyr::rename(ligand = from, receptor = to) %>% distinct(ligand, receptor, bonafide)
head(lr_network)
## # A tibble: 6 x 3
## ligand receptor bonafide
## <chr> <chr> <lgl>
## 1 CXCL1 CXCR2 TRUE
## 2 CXCL2 CXCR2 TRUE
## 3 CXCL3 CXCR2 TRUE
## 4 CXCL5 CXCR2 TRUE
## 5 PPBP CXCR2 TRUE
## 6 CXCL6 CXCR2 TRUE
人鼠基因轉(zhuǎn)換
if(organism == "mouse"){
lr_network = lr_network %>% mutate(ligand = convert_human_to_mouse_symbols(ligand), receptor = convert_human_to_mouse_symbols(receptor)) %>% drop_na()
colnames(ligand_target_matrix) = ligand_target_matrix %>% colnames() %>% convert_human_to_mouse_symbols()
rownames(ligand_target_matrix) = ligand_target_matrix %>% rownames() %>% convert_human_to_mouse_symbols()
ligand_target_matrix = ligand_target_matrix %>% .[!is.na(rownames(ligand_target_matrix)), !is.na(colnames(ligand_target_matrix))]
}
1炮温、Define the niches/microenvironments of interest
Each niche should have at least one “sender/niche” cell population and one “receiver/target” cell population (present in your expression data)
In this case study, we are interested to find differences in cell-cell interactions to hepatic macrophages in three different niches: 1) the Kupffer cell niche, 2) the bile-duct or lipid-associated macrophage niche, and 3) the capsule macrophage niche.
Based on imaging and spatial transcriptomics, the composition of each niche was defined as follows:(借助空間轉(zhuǎn)錄組定義生態(tài)位)
The receiver cell population in the Kupffer cell niche is the “KCs” cell type, the sender cell types are: “LSECs_portal”,“Hepatocytes_portal”, and “Stellate cells_portal”. The receiver cell population in the lipid-associated macrophage (MoMac2) niche is the “MoMac2” cell type, the sender cell types are: “Cholangiocytes”, and “Fibroblast 2”. The receiver cell population in the capsule macrophage (MoMac1) niche is the “MoMac1” cell type, the sender cell types are: “Capsule fibroblasts”, and “Mesothelial cells”.
niches = list(
"KC_niche" = list(
"sender" = c("LSECs_portal","Hepatocytes_portal","Stellate cells_portal"),
"receiver" = c("KCs")),
"MoMac2_niche" = list(
"sender" = c("Cholangiocytes","Fibroblast 2"),
"receiver" = c("MoMac2")),
"MoMac1_niche" = list(
"sender" = c("Capsule fibroblasts","Mesothelial cells"),
"receiver" = c("MoMac1"))
)
2、Calculate differential expression between the niches
determine DE between the different niches for both senders and receivers to define the DE of L-R pairs.
Calculate DE
The method to calculate the differential expression is here the standard Seurat Wilcoxon test, but this can be replaced if wanted by the user (only requirement: output tables DE_sender_processed and DE_receiver_processed should be in the same format as shown here).
DE will be calculated for each pairwise sender (or receiver) cell type comparision between the niches (so across niches, not within niche).
assay_oi = "SCT" # other possibilities: RNA,...
seurat_obj = PrepSCTFindMarkers(seurat_obj, assay = "SCT", verbose = FALSE)
DE_sender = calculate_niche_de(seurat_obj = seurat_obj %>% subset(features = lr_network$ligand %>% intersect(rownames(seurat_obj))), niches = niches, type = "sender", assay_oi = assay_oi) # only ligands important for sender cell types
## [1] "Calculate Sender DE between: LSECs_portal and Cholangiocytes"
## [2] "Calculate Sender DE between: LSECs_portal and Fibroblast 2"
## [3] "Calculate Sender DE between: LSECs_portal and Capsule fibroblasts"
## [4] "Calculate Sender DE between: LSECs_portal and Mesothelial cells"
## [1] "Calculate Sender DE between: Hepatocytes_portal and Cholangiocytes"
## [2] "Calculate Sender DE between: Hepatocytes_portal and Fibroblast 2"
## [3] "Calculate Sender DE between: Hepatocytes_portal and Capsule fibroblasts"
## [4] "Calculate Sender DE between: Hepatocytes_portal and Mesothelial cells"
## [1] "Calculate Sender DE between: Stellate cells_portal and Cholangiocytes"
## [2] "Calculate Sender DE between: Stellate cells_portal and Fibroblast 2"
## [3] "Calculate Sender DE between: Stellate cells_portal and Capsule fibroblasts"
## [4] "Calculate Sender DE between: Stellate cells_portal and Mesothelial cells"
## [1] "Calculate Sender DE between: Cholangiocytes and LSECs_portal"
## [2] "Calculate Sender DE between: Cholangiocytes and Hepatocytes_portal"
## [3] "Calculate Sender DE between: Cholangiocytes and Stellate cells_portal"
## [4] "Calculate Sender DE between: Cholangiocytes and Capsule fibroblasts"
## [5] "Calculate Sender DE between: Cholangiocytes and Mesothelial cells"
## [1] "Calculate Sender DE between: Fibroblast 2 and LSECs_portal"
## [2] "Calculate Sender DE between: Fibroblast 2 and Hepatocytes_portal"
## [3] "Calculate Sender DE between: Fibroblast 2 and Stellate cells_portal"
## [4] "Calculate Sender DE between: Fibroblast 2 and Capsule fibroblasts"
## [5] "Calculate Sender DE between: Fibroblast 2 and Mesothelial cells"
## [1] "Calculate Sender DE between: Capsule fibroblasts and LSECs_portal"
## [2] "Calculate Sender DE between: Capsule fibroblasts and Hepatocytes_portal"
## [3] "Calculate Sender DE between: Capsule fibroblasts and Stellate cells_portal"
## [4] "Calculate Sender DE between: Capsule fibroblasts and Cholangiocytes"
## [5] "Calculate Sender DE between: Capsule fibroblasts and Fibroblast 2"
## [1] "Calculate Sender DE between: Mesothelial cells and LSECs_portal"
## [2] "Calculate Sender DE between: Mesothelial cells and Hepatocytes_portal"
## [3] "Calculate Sender DE between: Mesothelial cells and Stellate cells_portal"
## [4] "Calculate Sender DE between: Mesothelial cells and Cholangiocytes"
## [5] "Calculate Sender DE between: Mesothelial cells and Fibroblast 2"
DE_receiver = calculate_niche_de(seurat_obj = seurat_obj %>% subset(features = lr_network$receptor %>% unique()), niches = niches, type = "receiver", assay_oi = assay_oi) # only receptors now, later on: DE analysis to find targets
## # A tibble: 3 x 2
## receiver receiver_other_niche
## <chr> <chr>
## 1 KCs MoMac2
## 2 KCs MoMac1
## 3 MoMac2 MoMac1
## [1] "Calculate receiver DE between: KCs and MoMac2" "Calculate receiver DE between: KCs and MoMac1"
## [1] "Calculate receiver DE between: MoMac2 and KCs" "Calculate receiver DE between: MoMac2 and MoMac1"
## [1] "Calculate receiver DE between: MoMac1 and KCs" "Calculate receiver DE between: MoMac1 and MoMac2"
DE_sender = DE_sender %>% mutate(avg_log2FC = ifelse(avg_log2FC == Inf, max(avg_log2FC[is.finite(avg_log2FC)]), ifelse(avg_log2FC == -Inf, min(avg_log2FC[is.finite(avg_log2FC)]), avg_log2FC)))
DE_receiver = DE_receiver %>% mutate(avg_log2FC = ifelse(avg_log2FC == Inf, max(avg_log2FC[is.finite(avg_log2FC)]), ifelse(avg_log2FC == -Inf, min(avg_log2FC[is.finite(avg_log2FC)]), avg_log2FC)))
注意這里的差異分析牵舵，任意兩對(duì)的細(xì)胞類型都納入分析柒啤。
Process DE results
expression_pct = 0.10 ###這個(gè)值最好大一點(diǎn)
DE_sender_processed = process_niche_de(DE_table = DE_sender, niches = niches, expression_pct = expression_pct, type = "sender")
DE_receiver_processed = process_niche_de(DE_table = DE_receiver, niches = niches, expression_pct = expression_pct, type = "receiver")
Combine sender-receiver DE based on L-R pairs
specificity_score_LR_pairs = "min_lfc"
DE_sender_receiver = combine_sender_receiver_de(DE_sender_processed, DE_receiver_processed, lr_network, specificity_score = specificity_score_LR_pairs)
3、Optional: Calculate differential expression between the different spatial regions
為了改進(jìn)細(xì)胞-細(xì)胞相互作用預(yù)測(cè)畸颅，如果可能且適用担巩，可以考慮空間信息∶怀矗空間信息可以來自顯微鏡數(shù)據(jù)涛癌，也可以來自空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)，例如 Visium送火。
有幾種方法可以將空間信息合并到差分 NicheNet pipeline中拳话。首先，如果細(xì)胞類型位于相同的空間位置种吸，則只能將它們視為屬于同一生態(tài)位弃衍。另一種方法是在優(yōu)先框架中包括一種細(xì)胞類型內(nèi)配體-受體對(duì)的空間差異表達(dá)。
例如：有一個(gè)細(xì)胞類型 X坚俗，位于區(qū)域 A 和 B镜盯，想研究區(qū)域 A 的細(xì)胞間通信。首先在生態(tài)位定義中僅添加區(qū)域 A 的 celltypeX猖败，然后計(jì)算 celltypeX-regionA 之間的 DE和 celltypeX-regionB 為 regionA 特異性配體提供更高的優(yōu)先權(quán)形耗。
在本案例研究中，感興趣的區(qū)域是肝臟的門靜脈周圍區(qū)域辙浑，因?yàn)樾∈笾械?KCs 主要位于門靜脈周圍區(qū)域激涤。因此，與中心周圍區(qū)域相比判呕，分析將賦予在 KCs 生態(tài)位細(xì)胞中更高表達(dá)的配體的權(quán)重倦踢。
include_spatial_info_sender = TRUE # if not spatial info to include: put this to false
include_spatial_info_receiver = FALSE # if spatial info to include: put this to true
spatial_info = tibble(celltype_region_oi = c("LSECs_portal","Hepatocytes_portal","Stellate cells_portal"),
celltype_other_region = c("LSECs_central","Hepatocytes_central","Stellate cells_central")
) %>%
mutate(niche = "KC_niche", celltype_type = "sender")
specificity_score_spatial = "lfc"
# this is how this should be defined if you don't have spatial info
# mock spatial info
if(include_spatial_info_sender == FALSE & include_spatial_info_receiver == FALSE){
spatial_info = tibble(celltype_region_oi = NA, celltype_other_region = NA) %>% mutate(niche = niches %>% names() %>% head(1), celltype_type = "sender")
}
if(include_spatial_info_sender == TRUE){
sender_spatial_DE = calculate_spatial_DE(seurat_obj = seurat_obj %>% subset(features = lr_network$ligand %>% unique()), spatial_info = spatial_info %>% filter(celltype_type == "sender"), assay_oi = assay_oi)
sender_spatial_DE_processed = process_spatial_de(DE_table = sender_spatial_DE, type = "sender", lr_network = lr_network, expression_pct = expression_pct, specificity_score = specificity_score_spatial)
# add a neutral spatial score for sender celltypes in which the spatial is not known / not of importance
sender_spatial_DE_others = get_non_spatial_de(niches = niches, spatial_info = spatial_info, type = "sender", lr_network = lr_network)
sender_spatial_DE_processed = sender_spatial_DE_processed %>% bind_rows(sender_spatial_DE_others)
sender_spatial_DE_processed = sender_spatial_DE_processed %>% mutate(scaled_ligand_score_spatial = scale_quantile_adapted(ligand_score_spatial))
} else {
# # add a neutral spatial score for all sender celltypes (for none of them, spatial is relevant in this case)
sender_spatial_DE_processed = get_non_spatial_de(niches = niches, spatial_info = spatial_info, type = "sender", lr_network = lr_network)
sender_spatial_DE_processed = sender_spatial_DE_processed %>% mutate(scaled_ligand_score_spatial = scale_quantile_adapted(ligand_score_spatial))
}
## [1] "Calculate Spatial DE between: LSECs_portal and LSECs_central"
## [1] "Calculate Spatial DE between: Hepatocytes_portal and Hepatocytes_central"
## [1] "Calculate Spatial DE between: Stellate cells_portal and Stellate cells_central"
if(include_spatial_info_receiver == TRUE){
receiver_spatial_DE = calculate_spatial_DE(seurat_obj = seurat_obj %>% subset(features = lr_network$receptor %>% unique()), spatial_info = spatial_info %>% filter(celltype_type == "receiver"), assay_oi = assay_oi)
receiver_spatial_DE_processed = process_spatial_de(DE_table = receiver_spatial_DE, type = "receiver", lr_network = lr_network, expression_pct = expression_pct, specificity_score = specificity_score_spatial)
# add a neutral spatial score for receiver celltypes in which the spatial is not known / not of importance
receiver_spatial_DE_others = get_non_spatial_de(niches = niches, spatial_info = spatial_info, type = "receiver", lr_network = lr_network)
receiver_spatial_DE_processed = receiver_spatial_DE_processed %>% bind_rows(receiver_spatial_DE_others)
receiver_spatial_DE_processed = receiver_spatial_DE_processed %>% mutate(scaled_receptor_score_spatial = scale_quantile_adapted(receptor_score_spatial))
} else {
# # add a neutral spatial score for all receiver celltypes (for none of them, spatial is relevant in this case)
receiver_spatial_DE_processed = get_non_spatial_de(niches = niches, spatial_info = spatial_info, type = "receiver", lr_network = lr_network)
receiver_spatial_DE_processed = receiver_spatial_DE_processed %>% mutate(scaled_receptor_score_spatial = scale_quantile_adapted(receptor_score_spatial))
}
4. Calculate ligand activities and infer active ligand-target links
在這一步中，將預(yù)測(cè)每個(gè)配體在不同生態(tài)位中的每種受體細(xì)胞類型的配體活性侠草。 這類似于在正常 NicheNet 管道中進(jìn)行的配體活性分析辱挥。
為了計(jì)算配體活性，首先需要為每個(gè)生態(tài)位定義一個(gè)感興趣的geneset边涕。 在本案例研究中晤碘，與膠囊和膽管巨噬細(xì)胞相比褂微，Kupffer 細(xì)胞生態(tài)位感興趣的基因組是 Kupffer 細(xì)胞中上調(diào)的基因。 與膠囊巨噬細(xì)胞和枯否細(xì)胞相比园爷，膽管巨噬細(xì)胞生態(tài)位的基因組是膽管巨噬細(xì)胞中上調(diào)的基因宠蚂。 同樣對(duì)于感興趣的膠囊巨噬細(xì)胞基因組。
也可以認(rèn)為定義感興趣的geneset
lfc_cutoff = 0.15 # recommended for 10x as min_lfc cutoff.
specificity_score_targets = "min_lfc"
DE_receiver_targets = calculate_niche_de_targets(seurat_obj = seurat_obj, niches = niches, lfc_cutoff = lfc_cutoff, expression_pct = expression_pct, assay_oi = assay_oi)
## [1] "Calculate receiver DE between: KCs and MoMac2" "Calculate receiver DE between: KCs and MoMac1"
## [1] "Calculate receiver DE between: MoMac2 and KCs" "Calculate receiver DE between: MoMac2 and MoMac1"
## [1] "Calculate receiver DE between: MoMac1 and KCs" "Calculate receiver DE between: MoMac1 and MoMac2"
DE_receiver_processed_targets = process_receiver_target_de(DE_receiver_targets = DE_receiver_targets, niches = niches, expression_pct = expression_pct, specificity_score = specificity_score_targets)
background = DE_receiver_processed_targets %>% pull(target) %>% unique()
geneset_KC = DE_receiver_processed_targets %>% filter(receiver == niches$KC_niche$receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1) %>% pull(target) %>% unique()
geneset_MoMac2 = DE_receiver_processed_targets %>% filter(receiver == niches$MoMac2_niche$receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1) %>% pull(target) %>% unique()
geneset_MoMac1 = DE_receiver_processed_targets %>% filter(receiver == niches$MoMac1_niche$receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1) %>% pull(target) %>% unique()
# Good idea to check which genes will be left out of the ligand activity analysis (=when not present in the rownames of the ligand-target matrix).
# If many genes are left out, this might point to some issue in the gene naming (eg gene aliases and old gene symbols, bad human-mouse mapping)
geneset_KC %>% setdiff(rownames(ligand_target_matrix))
## [1] "Fcna" "Wfdc17" "AW112010" "mt-Co1" "mt-Nd2" "C4b" "Adgre4" "mt-Co3"
## [9] "Pira2" "mt-Co2" "mt-Nd4" "mt-Atp6" "mt-Nd1" "mt-Nd3" "Ear2" "2900097C17Rik"
## [17] "Iigp1" "Trim30a" "B430306N03Rik" "mt-Cytb" "Pilrb2" "Anapc15" "Arf2" "Gbp8"
## [25] "AC149090.1" "Cd209f" "Xlr" "Ifitm6"
geneset_MoMac2 %>% setdiff(rownames(ligand_target_matrix))
## [1] "Chil3" "Lyz1" "Ccl9" "Tmsb10" "Ly6c2" "Gm21188" "Gm10076" "Ms4a6c"
## [9] "Calm3" "Atp5e" "Ftl1-ps1" "S100a11" "Clec4a3" "Snrpe" "Cox6c" "Ly6i"
## [17] "1810058I24Rik" "Rpl34" "Aph1c" "Atp5o.1"
geneset_MoMac1 %>% setdiff(rownames(ligand_target_matrix))
## [1] "H2-Ab1" "Malat1" "H2-Aa" "Hspa1b" "Gm26522" "Ly6a" "H2-D1" "Klra2" "Bcl2a1d" "Kcnq1ot1"
length(geneset_KC)
## [1] 443
length(geneset_MoMac2)
## [1] 339
length(geneset_MoMac1)
## [1] 84
It is always useful to check the number of genes in the geneset before doing the ligand activity analysis. We recommend having between 20 and 1000 genes in the geneset of interest, and a background of at least 5000 genes for a proper ligand activity analysis. If you retrieve too many DE genes, it is recommended to use a higher lfc_cutoff threshold. We recommend using a cutoff of 0.15 if you have > 2 receiver cells/niches to compare and use the min_lfc as specificity score. If you have only 2 receivers/niche, we recommend using a higher threshold (such as using 0.25). If you have single-cell data like Smart-seq2 with high sequencing depth, we recommend to also use higher threshold.
top_n_target = 250
niche_geneset_list = list(
"KC_niche" = list(
"receiver" = "KCs",
"geneset" = geneset_KC,
"background" = background),
"MoMac1_niche" = list(
"receiver" = "MoMac1",
"geneset" = geneset_MoMac1 ,
"background" = background),
"MoMac2_niche" = list(
"receiver" = "MoMac2",
"geneset" = geneset_MoMac2 ,
"background" = background)
)
ligand_activities_targets = get_ligand_activities_targets(niche_geneset_list = niche_geneset_list, ligand_target_matrix = ligand_target_matrix, top_n_target = top_n_target)
## [1] "Calculate Ligand activities for: KCs"
## [1] "Calculate Ligand activities for: MoMac1"
## [1] "Calculate Ligand activities for: MoMac2"
5. Calculate (scaled) expression of ligands, receptors and targets across cell types of interest (log expression values and expression fractions)
在這一步中童社，將計(jì)算所有感興趣的細(xì)胞類型的配體求厕、受體和靶基因的平均（縮放）表達(dá)和表達(dá)分?jǐn)?shù)。 現(xiàn)在扰楼，這通過 Seurat 的 DotPlot 函數(shù)進(jìn)行了分析呀癣，但當(dāng)然也可以通過其他方式完成。
features_oi = union(lr_network$ligand, lr_network$receptor) %>% union(ligand_activities_targets$target) %>% setdiff(NA)
dotplot = suppressWarnings(Seurat::DotPlot(seurat_obj %>% subset(idents = niches %>% unlist() %>% unique()), features = features_oi, assay = assay_oi))
exprs_tbl = dotplot$data %>% as_tibble()
exprs_tbl = exprs_tbl %>% rename(celltype = id, gene = features.plot, expression = avg.exp, expression_scaled = avg.exp.scaled, fraction = pct.exp) %>%
mutate(fraction = fraction/100) %>% as_tibble() %>% select(celltype, gene, expression, expression_scaled, fraction) %>% distinct() %>% arrange(gene) %>% mutate(gene = as.character(gene))
exprs_tbl_ligand = exprs_tbl %>% filter(gene %in% lr_network$ligand) %>% rename(sender = celltype, ligand = gene, ligand_expression = expression, ligand_expression_scaled = expression_scaled, ligand_fraction = fraction)
exprs_tbl_receptor = exprs_tbl %>% filter(gene %in% lr_network$receptor) %>% rename(receiver = celltype, receptor = gene, receptor_expression = expression, receptor_expression_scaled = expression_scaled, receptor_fraction = fraction)
exprs_tbl_target = exprs_tbl %>% filter(gene %in% ligand_activities_targets$target) %>% rename(receiver = celltype, target = gene, target_expression = expression, target_expression_scaled = expression_scaled, target_fraction = fraction)
exprs_tbl_ligand = exprs_tbl_ligand %>% mutate(scaled_ligand_expression_scaled = scale_quantile_adapted(ligand_expression_scaled)) %>% mutate(ligand_fraction_adapted = ligand_fraction) %>% mutate_cond(ligand_fraction >= expression_pct, ligand_fraction_adapted = expression_pct) %>% mutate(scaled_ligand_fraction_adapted = scale_quantile_adapted(ligand_fraction_adapted))
exprs_tbl_receptor = exprs_tbl_receptor %>% mutate(scaled_receptor_expression_scaled = scale_quantile_adapted(receptor_expression_scaled)) %>% mutate(receptor_fraction_adapted = receptor_fraction) %>% mutate_cond(receptor_fraction >= expression_pct, receptor_fraction_adapted = expression_pct) %>% mutate(scaled_receptor_fraction_adapted = scale_quantile_adapted(receptor_fraction_adapted))
6. Expression fraction and receptor
在這一步中弦赖，將根據(jù)受體的表達(dá)強(qiáng)度對(duì)配體-受體相互作用進(jìn)行評(píng)分项栏，這樣就可以對(duì)特定細(xì)胞類型中特定配體的最強(qiáng)表達(dá)受體給予更高的分?jǐn)?shù)。 這不會(huì)影響以后單個(gè)配體的等級(jí)蹬竖，但將有助于確定每個(gè)配體最重要的受體的優(yōu)先級(jí).
exprs_sender_receiver = lr_network %>%
inner_join(exprs_tbl_ligand, by = c("ligand")) %>%
inner_join(exprs_tbl_receptor, by = c("receptor")) %>% inner_join(DE_sender_receiver %>% distinct(niche, sender, receiver))
ligand_scaled_receptor_expression_fraction_df = exprs_sender_receiver %>% group_by(ligand, receiver) %>% mutate(rank_receptor_expression = dense_rank(receptor_expression), rank_receptor_fraction = dense_rank(receptor_fraction)) %>% mutate(ligand_scaled_receptor_expression_fraction = 0.5*( (rank_receptor_fraction / max(rank_receptor_fraction)) + ((rank_receptor_expression / max(rank_receptor_expression))) ) ) %>% distinct(ligand, receptor, receiver, ligand_scaled_receptor_expression_fraction, bonafide) %>% distinct() %>% ungroup()
7. Prioritization of ligand-receptor and ligand-target links
在這一步中忘嫉，將結(jié)合上述所有計(jì)算信息來優(yōu)先考慮配體-受體-目標(biāo)link。 在 0 和 1 之間縮放每個(gè)感興趣的屬性案腺，最終的優(yōu)先級(jí)分?jǐn)?shù)是所有感興趣屬性的縮放分?jǐn)?shù)的加權(quán)和庆冕。
We provide the user the option to consider the following properties for prioritization (of which the weights are defined in prioritizing_weights) :
Ligand DE score: niche-specific expression of the ligand: by default, this the minimum logFC between the sender of interest and all the senders of the other niche(s). The higher the min logFC, the higher the niche-specificity of the ligand. Therefore we recommend to give this factor a very high weight. prioritizing_weights argument: "scaled_ligand_score". Recommended weight: 5 (at least 1, max 5).
Scaled ligand expression: scaled expression of a ligand in one sender compared to the other cell types in the dataset. This might be useful to rescue potentially interesting ligands that have a high scaled expression value, but a relatively small min logFC compared to the other niche. One reason why this logFC might be small occurs when (some) genes are not picked up efficiently by the used sequencing technology (or other reasons for low RNA expression of ligands). For example, we have observed that many ligands from the Tgf-beta/BMP family are not picked up efficiently with single-nuclei RNA sequencing compared to single-cell sequencing. prioritizing_weights argument: "scaled_ligand_expression_scaled". Recommended weight: 1 (unless technical reason for lower gene detection such as while using Nuc-seq: then recommended to use a higher weight: 2).
Ligand expression fraction: Ligands that are expressed in a smaller fraction of cells of a cell type than defined by exprs_cutoff(default: 0.10) will get a lower ranking, proportional to their fraction (eg ligand expressed in 9% of cells will be ranked higher than ligand expressed in 0.5% of cells). We opted for this weighting based on fraction, instead of removing ligands that are not expressed in more cells than this cutoff, because some interesting ligands could be removed that way. Fraction of expression is not taken into account for the prioritization if it is already higher than the cutoff. prioritizing_weights argument: "ligand_fraction". Recommended weight: 1.
Ligand spatial DE score: spatial expression specificity of the ligand. If the niche of interest is at a specific tissue location, but some of the sender cell types of that niche are also present in other locations, it can be very informative to further prioritize ligands of that sender by looking how they are DE between the spatial location of interest compared to the other locations. prioritizing_weights argument: "scaled_ligand_score_spatial". Recommended weight: 2 (or 0 if not applicable).
Receptor DE score: niche-specific expression of the receptor: by default, this the minimum logFC between the receiver of interest and all the receiver of the other niche(s). The higher the min logFC, the higher the niche-specificity of the receptor. Based on our experience, we don’t suggest to give this as high importance as the ligand DE, but this might depend on the specific case study. prioritizing_weights argument: "scaled_receptor_score". Recommended weight: 0.5 (at least 0.5, and lower than "scaled_ligand_score").
Scaled receptor expression: scaled expression of a receptor in one receiver compared to the other cell types in the dataset. This might be useful to rescue potentially interesting receptors that have a high scaled expression value, but a relatively small min logFC compared to the other niche. One reason why this logFC might be small occurs when (some) genes are not picked up efficiently by the used sequencing technology. prioritizing_weights argument: "scaled_receptor_expression_scaled". Recommended weight: 0.5.
Receptor expression fraction: Receptors that are expressed in a smaller fraction of cells of a cell type than defined by exprs_cutoff(default: 0.10) will get a lower ranking, proportional to their fraction (eg receptor expressed in 9% of cells will be ranked higher than receptor expressed in 0.5% of cells). We opted for this weighting based on fraction, instead of removing receptors that are not expressed in more cells than this cutoff, because some interesting receptors could be removed that way. Fraction of expression is not taken into account for the prioritization if it is already higher than the cutoff. prioritizing_weights argument: "receptor_fraction". Recommended weight: 1.
Receptor expression strength: this factor let us give higher weights to the most highly expressed receptor of a ligand in the receiver. This let us rank higher one member of a receptor family if it higher expressed than the other members. prioritizing_weights argument: "ligand_scaled_receptor_expression_fraction". Recommended value: 1 (minimum: 0.5).
Receptor spatial DE score: spatial expression specificity of the receptor. If the niche of interest is at a specific tissue location, but the receiver cell type of that niche is also present in other locations, it can be very informative to further prioritize receptors of that receiver by looking how they are DE between the spatial location of interest compared to the other locations. prioritizing_weights argument: "scaled_receptor_score_spatial". Recommended weight: 1 (or 0 if not applicable).
Absolute ligand activity: to further prioritize ligand-receptor pairs based on their predicted effect of the ligand-receptor interaction on the gene expression in the receiver cell type - absolute ligand activity accords to ‘a(chǎn)bsolute’ enrichment of target genes of a ligand within the affected receiver genes. prioritizing_weights argument: "scaled_activity". Recommended weight: 0, unless absolute enrichment of target genes is of specific interest.
Normalized ligand activity: to further prioritize ligand-receptor pairs based on their predicted effect of the ligand-receptor interaction on the gene expression in the receiver cell type - normalization of activity is done because we found that some datasets/conditions/niches have higher baseline activity values than others - normalized ligand activity accords to ‘relative’ enrichment of target genes of a ligand within the affected receiver genes. prioritizing_weights argument: "scaled_activity_normalized". Recommended weight: at least 1.
Prior knowledge quality of the L-R interaction: the NicheNet LR network consists of two types of interactions: L-R pairs documented in curated databases, and L-R pairs predicted based on gene annotation and PPIs. The former are categorized as ‘bona fide’ interactions. To rank bona fide interactions higher, but not exlude potentially interesting non-bona-fide ones, we give bona fide interactions a score of 1, and non-bona-fide interactions a score fof 0.5. prioritizing_weights argument: "bona_fide" Recommend weight: at least 1.
prioritizing_weights = c("scaled_ligand_score" = 5,
"scaled_ligand_expression_scaled" = 1,
"ligand_fraction" = 1,
"scaled_ligand_score_spatial" = 2,
"scaled_receptor_score" = 0.5,
"scaled_receptor_expression_scaled" = 0.5,
"receptor_fraction" = 1,
"ligand_scaled_receptor_expression_fraction" = 1,
"scaled_receptor_score_spatial" = 0,
"scaled_activity" = 0,
"scaled_activity_normalized" = 1,
"bona_fide" = 1)
output = list(DE_sender_receiver = DE_sender_receiver, ligand_scaled_receptor_expression_fraction_df = ligand_scaled_receptor_expression_fraction_df, sender_spatial_DE_processed = sender_spatial_DE_processed, receiver_spatial_DE_processed = receiver_spatial_DE_processed,
ligand_activities_targets = ligand_activities_targets, DE_receiver_processed_targets = DE_receiver_processed_targets, exprs_tbl_ligand = exprs_tbl_ligand, exprs_tbl_receptor = exprs_tbl_receptor, exprs_tbl_target = exprs_tbl_target)
prioritization_tables = get_prioritization_tables(output, prioritizing_weights)
prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(receiver == niches[[1]]$receiver) %>% head(10)
## # A tibble: 10 x 37
## niche receiver sender ligand_receptor ligand receptor bonafide ligand_score ligand_signific~ ligand_present ligand_expressi~
## <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl>
## 1 KC_niche KCs Hepatocytes~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE 3.18 1 1 14.7
## 2 KC_niche KCs Hepatocytes~ Apoa1--Msr1 Apoa1 Msr1 FALSE 3.18 1 1 14.7
## 3 KC_niche KCs Hepatocytes~ Apoa1--Abca1 Apoa1 Abca1 FALSE 3.18 1 1 14.7
## 4 KC_niche KCs Hepatocytes~ Apoa1--Scarb1 Apoa1 Scarb1 FALSE 3.18 1 1 14.7
## 5 KC_niche KCs Hepatocytes~ Apoa1--Derl1 Apoa1 Derl1 FALSE 3.18 1 1 14.7
## 6 KC_niche KCs Hepatocytes~ Serpina1a--Lrp1 Serpina~ Lrp1 TRUE 2.64 1 1 6.97
## 7 KC_niche KCs Hepatocytes~ Apoa1--Atp5b Apoa1 Atp5b FALSE 3.18 1 1 14.7
## 8 KC_niche KCs Hepatocytes~ Trf--Tfrc Trf Tfrc TRUE 1.61 1 1 6.19
## 9 KC_niche KCs Hepatocytes~ Apoa1--Cd36 Apoa1 Cd36 FALSE 3.18 1 1 14.7
## 10 KC_niche KCs LSECs_portal Cxcl10--Fpr1 Cxcl10 Fpr1 FALSE 1.66 1 1 2.35
## # ... with 26 more variables: ligand_expression_scaled <dbl>, ligand_fraction <dbl>, ligand_score_spatial <dbl>, receptor_score <dbl>,
## # receptor_significant <dbl>, receptor_present <dbl>, receptor_expression <dbl>, receptor_expression_scaled <dbl>,
## # receptor_fraction <dbl>, receptor_score_spatial <dbl>, ligand_scaled_receptor_expression_fraction <dbl>,
## # avg_score_ligand_receptor <dbl>, activity <dbl>, activity_normalized <dbl>, scaled_ligand_score <dbl>,
## # scaled_ligand_expression_scaled <dbl>, scaled_receptor_score <dbl>, scaled_receptor_expression_scaled <dbl>,
## # scaled_avg_score_ligand_receptor <dbl>, scaled_ligand_score_spatial <dbl>, scaled_receptor_score_spatial <dbl>,
## # scaled_ligand_fraction_adapted <dbl>, scaled_receptor_fraction_adapted <dbl>, scaled_activity <dbl>, ...
prioritization_tables$prioritization_tbl_ligand_target %>% filter(receiver == niches[[1]]$receiver) %>% head(10)
## # A tibble: 10 x 20
## niche receiver sender ligand_receptor ligand receptor bonafide target target_score target_signific~ target_present target_expressi~
## <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Abca1 0.197 1 1 0.979
## 2 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Actb 0.279 1 1 21.6
## 3 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Ehd1 0.272 1 1 0.402
## 4 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Hmox1 1.16 1 1 5.23
## 5 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Sgk1 0.265 1 1 0.629
## 6 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Tcf7l2 0.811 1 1 1.32
## 7 KC_niche KCs Hepato~ Apoa1--Lrp1 Apoa1 Lrp1 FALSE Tsc22~ 0.263 1 1 0.635
## 8 KC_niche KCs Hepato~ Apoa1--Msr1 Apoa1 Msr1 FALSE Abca1 0.197 1 1 0.979
## 9 KC_niche KCs Hepato~ Apoa1--Msr1 Apoa1 Msr1 FALSE Actb 0.279 1 1 21.6
## 10 KC_niche KCs Hepato~ Apoa1--Msr1 Apoa1 Msr1 FALSE Ehd1 0.272 1 1 0.402
## # ... with 8 more variables: target_expression_scaled <dbl>, target_fraction <dbl>, ligand_target_weight <dbl>, activity <dbl>,
## # activity_normalized <dbl>, scaled_activity <dbl>, scaled_activity_normalized <dbl>, prioritization_score <dbl>
prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(receiver == niches[[2]]$receiver) %>% head(10)
## # A tibble: 10 x 37
## niche receiver sender ligand_receptor ligand receptor bonafide ligand_score ligand_significa~ ligand_present ligand_expressi~
## <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl>
## 1 MoMac2_niche MoMac2 Cholangi~ Spp1--Itga4 Spp1 Itga4 TRUE 6.09 1 1 72.4
## 2 MoMac2_niche MoMac2 Cholangi~ Spp1--Cd44 Spp1 Cd44 TRUE 6.09 1 1 72.4
## 3 MoMac2_niche MoMac2 Cholangi~ Spp1--Itgb5 Spp1 Itgb5 TRUE 6.09 1 1 72.4
## 4 MoMac2_niche MoMac2 Cholangi~ Spp1--Itgav Spp1 Itgav TRUE 6.09 1 1 72.4
## 5 MoMac2_niche MoMac2 Cholangi~ Spp1--Itgb1 Spp1 Itgb1 TRUE 6.09 1 1 72.4
## 6 MoMac2_niche MoMac2 Cholangi~ Spp1--Itga9 Spp1 Itga9 TRUE 6.09 1 1 72.4
## 7 MoMac2_niche MoMac2 Cholangi~ Spp1--Ncstn Spp1 Ncstn FALSE 6.09 1 1 72.4
## 8 MoMac2_niche MoMac2 Cholangi~ Spp1--Itga5 Spp1 Itga5 FALSE 6.09 1 1 72.4
## 9 MoMac2_niche MoMac2 Fibrobla~ Lama2--Rpsa Lama2 Rpsa TRUE 1.51 1 1 3.19
## 10 MoMac2_niche MoMac2 Cholangi~ Cyr61--Itgb2 Cyr61 Itgb2 TRUE 0.812 1 1 3.11
## # ... with 26 more variables: ligand_expression_scaled <dbl>, ligand_fraction <dbl>, ligand_score_spatial <dbl>, receptor_score <dbl>,
## # receptor_significant <dbl>, receptor_present <dbl>, receptor_expression <dbl>, receptor_expression_scaled <dbl>,
## # receptor_fraction <dbl>, receptor_score_spatial <dbl>, ligand_scaled_receptor_expression_fraction <dbl>,
## # avg_score_ligand_receptor <dbl>, activity <dbl>, activity_normalized <dbl>, scaled_ligand_score <dbl>,
## # scaled_ligand_expression_scaled <dbl>, scaled_receptor_score <dbl>, scaled_receptor_expression_scaled <dbl>,
## # scaled_avg_score_ligand_receptor <dbl>, scaled_ligand_score_spatial <dbl>, scaled_receptor_score_spatial <dbl>,
## # scaled_ligand_fraction_adapted <dbl>, scaled_receptor_fraction_adapted <dbl>, scaled_activity <dbl>, ...
prioritization_tables$prioritization_tbl_ligand_target %>% filter(receiver == niches[[2]]$receiver) %>% head(10)
## # A tibble: 10 x 20
## niche receiver sender ligand_receptor ligand receptor bonafide target target_score target_signific~ target_present target_expressi~
## <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Ahnak 1.05 1 1 1.36
## 2 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Cdkn1a 0.609 1 1 0.801
## 3 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Cxcr4 0.374 1 1 0.717
## 4 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Dhrs3 0.371 1 1 0.743
## 5 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Fn1 0.360 1 1 0.456
## 6 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Gadd4~ 0.180 1 1 0.474
## 7 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Gapdh 0.656 1 1 3.27
## 8 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Gdf15 0.479 1 1 0.521
## 9 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Gsn 0.221 1 1 0.647
## 10 MoMac2~ MoMac2 Cholang~ Spp1--Itga4 Spp1 Itga4 TRUE Plec 0.154 1 1 0.164
## # ... with 8 more variables: target_expression_scaled <dbl>, target_fraction <dbl>, ligand_target_weight <dbl>, activity <dbl>,
## # activity_normalized <dbl>, scaled_activity <dbl>, scaled_activity_normalized <dbl>, prioritization_score <dbl>
prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(receiver == niches[[3]]$receiver) %>% head(10)
## # A tibble: 10 x 37
## niche receiver sender ligand_receptor ligand receptor bonafide ligand_score ligand_signific~ ligand_present ligand_expressi~
## <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl>
## 1 MoMac1_niche MoMac1 Mesotheli~ C3--C3ar1 C3 C3ar1 TRUE 3.52 1 1 22.6
## 2 MoMac1_niche MoMac1 Capsule f~ C3--C3ar1 C3 C3ar1 TRUE 3.42 1 1 20.9
## 3 MoMac1_niche MoMac1 Mesotheli~ C3--Itgb2 C3 Itgb2 TRUE 3.52 1 1 22.6
## 4 MoMac1_niche MoMac1 Mesotheli~ C3--Itgax C3 Itgax TRUE 3.52 1 1 22.6
## 5 MoMac1_niche MoMac1 Mesotheli~ C3--Lrp1 C3 Lrp1 TRUE 3.52 1 1 22.6
## 6 MoMac1_niche MoMac1 Capsule f~ C3--Itgb2 C3 Itgb2 TRUE 3.42 1 1 20.9
## 7 MoMac1_niche MoMac1 Capsule f~ C3--Itgax C3 Itgax TRUE 3.42 1 1 20.9
## 8 MoMac1_niche MoMac1 Capsule f~ C3--Lrp1 C3 Lrp1 TRUE 3.42 1 1 20.9
## 9 MoMac1_niche MoMac1 Capsule f~ Rarres2--Cmklr1 Rarre~ Cmklr1 TRUE 2.50 1 1 15.8
## 10 MoMac1_niche MoMac1 Mesotheli~ C3--Ccr5 C3 Ccr5 FALSE 3.52 1 1 22.6
## # ... with 26 more variables: ligand_expression_scaled <dbl>, ligand_fraction <dbl>, ligand_score_spatial <dbl>, receptor_score <dbl>,
## # receptor_significant <dbl>, receptor_present <dbl>, receptor_expression <dbl>, receptor_expression_scaled <dbl>,
## # receptor_fraction <dbl>, receptor_score_spatial <dbl>, ligand_scaled_receptor_expression_fraction <dbl>,
## # avg_score_ligand_receptor <dbl>, activity <dbl>, activity_normalized <dbl>, scaled_ligand_score <dbl>,
## # scaled_ligand_expression_scaled <dbl>, scaled_receptor_score <dbl>, scaled_receptor_expression_scaled <dbl>,
## # scaled_avg_score_ligand_receptor <dbl>, scaled_ligand_score_spatial <dbl>, scaled_receptor_score_spatial <dbl>,
## # scaled_ligand_fraction_adapted <dbl>, scaled_receptor_fraction_adapted <dbl>, scaled_activity <dbl>, ...
prioritization_tables$prioritization_tbl_ligand_target %>% filter(receiver == niches[[3]]$receiver) %>% head(10)
## # A tibble: 10 x 20
## niche receiver sender ligand_receptor ligand receptor bonafide target target_score target_signific~ target_present target_expressi~
## <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Btg2 0.615 1 1 1.51
## 2 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Ccnd2 0.505 1 1 0.490
## 3 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Cdk6 0.221 1 1 0.320
## 4 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Ier5 0.396 1 1 1.16
## 5 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Il1b 0.956 1 1 3.74
## 6 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Jun 0.765 1 1 1.93
## 7 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Pdgfb 0.243 1 1 0.510
## 8 MoMac1~ MoMac1 Mesothe~ C3--C3ar1 C3 C3ar1 TRUE Ubc 0.306 1 1 2.16
## 9 MoMac1~ MoMac1 Capsule~ C3--C3ar1 C3 C3ar1 TRUE Btg2 0.615 1 1 1.51
## 10 MoMac1~ MoMac1 Capsule~ C3--C3ar1 C3 C3ar1 TRUE Ccnd2 0.505 1 1 0.490
## # ... with 8 more variables: target_expression_scaled <dbl>, target_fraction <dbl>, ligand_target_weight <dbl>, activity <dbl>,
## # activity_normalized <dbl>, scaled_activity <dbl>, scaled_activity_normalized <dbl>, prioritization_score <dbl>
prioritization_tables$prioritization_tbl_ligand_receptor = prioritization_tables$prioritization_tbl_ligand_receptor %>% mutate(receiver = factor(receiver, levels = c("KCs","MoMac1","MoMac2")), niche = factor(niche, levels = c("KC_niche","MoMac1_niche","MoMac2_niche")))
prioritization_tables$prioritization_tbl_ligand_target = prioritization_tables$prioritization_tbl_ligand_target %>% mutate(receiver = factor(receiver, levels = c("KCs","MoMac1","MoMac2")), niche = factor(niche, levels = c("KC_niche","MoMac1_niche","MoMac2_niche")))
8、可視化
Differential expression of ligand and expression
Before visualization, we need to define the most important ligand-receptor pairs per niche. We will do this by first determining for which niche the highest score is found for each ligand/ligand-receptor pair. And then getting the top 50 ligands per niche.
top_ligand_niche_df = prioritization_tables$prioritization_tbl_ligand_receptor %>% select(niche, sender, receiver, ligand, receptor, prioritization_score) %>% group_by(ligand) %>% top_n(1, prioritization_score) %>% ungroup() %>% select(ligand, receptor, niche) %>% rename(top_niche = niche)
top_ligand_receptor_niche_df = prioritization_tables$prioritization_tbl_ligand_receptor %>% select(niche, sender, receiver, ligand, receptor, prioritization_score) %>% group_by(ligand, receptor) %>% top_n(1, prioritization_score) %>% ungroup() %>% select(ligand, receptor, niche) %>% rename(top_niche = niche)
ligand_prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% select(niche, sender, receiver, ligand, prioritization_score) %>% group_by(ligand, niche) %>% top_n(1, prioritization_score) %>% ungroup() %>% distinct() %>% inner_join(top_ligand_niche_df) %>% filter(niche == top_niche) %>% group_by(niche) %>% top_n(50, prioritization_score) %>% ungroup() # get the top50 ligands per niche
receiver_oi = "KCs"
filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% pull(ligand) %>% unique()
prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand, receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup()
lfc_plot = make_ligand_receptor_lfc_plot(receiver_oi, prioritized_tbl_oi, prioritization_tables$prioritization_tbl_ligand_receptor, plot_legend = FALSE, heights = NULL, widths = NULL)
lfc_plot
Show the spatialDE as additional information
lfc_plot_spatial = make_ligand_receptor_lfc_spatial_plot(receiver_oi, prioritized_tbl_oi, prioritization_tables$prioritization_tbl_ligand_receptor, ligand_spatial = include_spatial_info_sender, receptor_spatial = include_spatial_info_receiver, plot_legend = FALSE, heights = NULL, widths = NULL)
lfc_plot_spatial
從這個(gè)圖中劈榨，您可以看到一些 KC niche配體访递，如 Dll4（由 LSEC）和 Il34（由星狀細(xì)胞）在門靜脈周圍 LSEC/星狀細(xì)胞中的表達(dá)高于中央周圍的。 這可能是有趣的信息同辣，因?yàn)橹?KC 主要位于門靜脈拷姿。 然而，其他配體旱函，如 Gdf2（星狀細(xì)胞）不優(yōu)先由門靜脈周圍星狀細(xì)胞表達(dá)响巢，但這并不意味著它們不有趣。 如下圖所示棒妨，該配體具有最高的配體活性之一踪古，這意味著其靶基因在 KC 特異性基因中存在很強(qiáng)的富集。
Ligand expression, activity and target genes
exprs_activity_target_plot = make_ligand_activity_target_exprs_plot(receiver_oi, prioritized_tbl_oi, prioritization_tables$prioritization_tbl_ligand_receptor, prioritization_tables$prioritization_tbl_ligand_target, output$exprs_tbl_ligand, output$exprs_tbl_target, lfc_cutoff, ligand_target_matrix, plot_legend = FALSE, heights = NULL, widths = NULL)
exprs_activity_target_plot$combined_plot
On this plot, we can see that some strongly DE ligand-receptor pairs in the KC niche, have also high scaled ligand activity on KCs - making them strong predictions for further validation. An example of this is Gdf2 and Bmp10, who bind the receptor Acvrl1 (ALK1). The role of Gdf2/Bmp10-Acvrl1 in KC development was experimentally validated in the Guilliams et al paper.
important: ligand-receptor pairs with both high differential expression and ligand activity (=target gene enrichment) are very interesting predictions as key regulators of your intercellular communication process of interest !
filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% top_n(20, prioritization_score) %>% pull(ligand) %>% unique()
prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand, receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup()
exprs_activity_target_plot = make_ligand_activity_target_exprs_plot(receiver_oi, prioritized_tbl_oi, prioritization_tables$prioritization_tbl_ligand_receptor, prioritization_tables$prioritization_tbl_ligand_target, output$exprs_tbl_ligand, output$exprs_tbl_target, lfc_cutoff, ligand_target_matrix, plot_legend = FALSE, heights = NULL, widths = NULL)
exprs_activity_target_plot$combined_plot
Circos plot of prioritized ligand-receptor pairs
filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% top_n(15, prioritization_score) %>% pull(ligand) %>% unique()
prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand, receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup()
colors_sender = brewer.pal(n = prioritized_tbl_oi$sender %>% unique() %>% sort() %>% length(), name = 'Spectral') %>% magrittr::set_names(prioritized_tbl_oi$sender %>% unique() %>% sort())
colors_receiver = c("lavender") %>% magrittr::set_names(prioritized_tbl_oi$receiver %>% unique() %>% sort())
circos_output = make_circos_lr(prioritized_tbl_oi, colors_sender, colors_receiver)
好了券腔，已經(jīng)分享給大家了伏穆，生活很好，有你更好纷纫，百度文庫出現(xiàn)了大量抄襲我的文章枕扫，對(duì)此我深表無奈，我寫的文章辱魁，別人掛上去賺錢烟瞧，抄襲可恥诗鸭，掛到百度文庫的人更可恥
