常用的細(xì)胞通訊軟件：
- CellphoneDB：是公開的人工校正的，儲(chǔ)存受體、配體以及兩種相互作用的數(shù)據(jù)庫(kù)。此外，還考慮了結(jié)構(gòu)組成，能夠描述異構(gòu)復(fù)合物。（配體-受體+多聚體）
- iTALK：通過平均表達(dá)量方式，篩選高表達(dá)的胚體和受體，根據(jù)結(jié)果作圈圖。（配體-受體）
- CellChat：CellChat將細(xì)胞的基因表達(dá)數(shù)據(jù)作為輸入，并結(jié)合配體受體及其輔助因子的相互作用來(lái)模擬細(xì)胞間通訊。（配體-受體+多聚體+輔因子）
- NicheNet // NicheNet多樣本分析 // 目標(biāo)基因的配體和靶基因活性預(yù)測(cè)：通過將相互作用細(xì)胞的表達(dá)數(shù)據(jù)與信號(hào)和基因調(diào)控網(wǎng)絡(luò)的先驗(yàn)知識(shí)相結(jié)合來(lái)預(yù)測(cè)相互作用細(xì)胞之間的配體-靶標(biāo)聯(lián)系的方法。( 配體-受體+信號(hào)通路）
附：NicheNet使用的常見問題匯總其它細(xì)胞互作軟件還包括
Celltalker，SingleCellSignalR，scTensor和SoptSC（這幾個(gè)也是基于配體-受體相互作用）
官網(wǎng)：https://github.com/saeyslab/nichenetr
Nature Methods原文
1. NicheNet介紹
1.1 NicheNet工作流程：
一般的預(yù)測(cè)細(xì)胞交互的軟件往往只考慮sender細(xì)胞的配體和receiver細(xì)胞的受體表達(dá)，但細(xì)胞交互過程除了配體-受體相互作用以外，還包含了receiver細(xì)胞接受信號(hào)后相關(guān)通路的激活。
NicheNet輸入基因表達(dá)數(shù)據(jù)，并將其與通過整合信號(hào)通路而構(gòu)建的模型相結(jié)合。不止是預(yù)測(cè)配體與受體的相互作用，還整合了細(xì)胞內(nèi)信號(hào)傳導(dǎo)。因此，NicheNet可以預(yù)測(cè)1)來(lái)自一或多種細(xì)胞中的配體（sender）影響了與之相互作用的細(xì)胞中哪些基因的表達(dá)和2)哪些靶基因受每種配體影響以及可能涉及哪些信號(hào)傳導(dǎo)介質(zhì)。
NicheNet工作流程圖
首先從公共數(shù)據(jù)庫(kù)中收集配體-受體配對(duì)信息、信號(hào)通路、基因調(diào)控網(wǎng)絡(luò)等信息，整合成配體主導(dǎo)的權(quán)重配體-靶基因調(diào)控模型。然后將可能受到細(xì)胞通訊影響的差異基因集輸入先驗(yàn)?zāi)Ｐ?，可以?jì)算與這些基因相關(guān)的配體的相關(guān)性系數(shù)。最后挑選根據(jù)相關(guān)性系數(shù)排行靠前的配體，依據(jù)先驗(yàn)數(shù)據(jù)推測(cè)與之匹配的受體、靶基因及下游信號(hào)網(wǎng)絡(luò)等信息。
NicheNet不同于大多數(shù)研究細(xì)胞間通訊的方法，它著眼于配體對(duì)下游基因調(diào)控作用。NicheNet可以預(yù)測(cè)哪些配體影響另一個(gè)細(xì)胞中的表達(dá)，哪些靶基因受到配體的影響以及哪些信號(hào)傳導(dǎo)可能參與其中。
1.2 NicheNet工作pipeline
- 在單細(xì)胞數(shù)據(jù)中定義一個(gè)“sender/niche”細(xì)胞群（配體細(xì)胞）和一個(gè)“receiver/target”細(xì)胞群（受體細(xì)胞）并確定這兩個(gè)細(xì)胞群都表達(dá)哪些基因。
- 定義一個(gè)感興趣的基因集：這些基因來(lái)自受體細(xì)胞群，是可能受到與其相互作用的細(xì)胞配體調(diào)控的基因集。（例如：case-control中的差異表達(dá)基因，也可以是細(xì)胞的signature或其他基因集）
- 定義一個(gè)潛在的配體集：這些配體由配體細(xì)胞群中高表達(dá)（如10%以上的細(xì)胞表達(dá)）并可以與受體細(xì)胞群表達(dá)的受體相結(jié)合（通過先驗(yàn)數(shù)據(jù)推斷）。
- 進(jìn)行NicheNet配體活性分析：其活性主要通過配體與受體細(xì)胞中的差異基因集的相關(guān)性進(jìn)行判斷
- 推斷在配體活性分析中的top-ranked配體所調(diào)控的top-predicted靶基因，以及與配體配對(duì)的受體。
NicheNet提供了一個(gè)三個(gè)功能相似的打包函數(shù)： nichenet_seuratobj_aggregate, nichenet_seuratobj_cluster_de and nichenet_seuratobj_aggregate_cluster_de.它們可以一步完成上述五步seurat對(duì)象的配體調(diào)控網(wǎng)絡(luò)分析。
1.3 NicheNet主要功能
Specific functionalities of this package include:
- assessing how well ligands expressed by a sender cell can predict changes in gene expression in the receiver cell
- prioritizing ligands based on their effect on gene expression
- inferring putative ligand-target links active in the system under study
- inferring potential signaling paths between ligands and target genes of interest: to generate causal hypotheses and check which data sources support the predictions
- validation of the prior ligand-target model
- construction of user-defined prior ligand-target models
- Moreover, we provide instructions on how to make intuitive visualizations of the main predictions (e.g., via circos plots).
2. Perform NicheNet analysis starting from a Seurat object
本文的演示數(shù)據(jù)集和代碼來(lái)自NicheNet官方分析單細(xì)胞數(shù)據(jù)的教程：https://github.com/saeyslab/nichenetr/blob/master/vignettes/seurat_wrapper.md
我們將使用Medaglia等人的小鼠NICHE-seq數(shù)據(jù)，探索淋巴細(xì)胞性脈絡(luò)膜腦膜炎病毒（LCMV）感染之前和之后72小時(shí)的腹股溝淋巴結(jié)T細(xì)胞區(qū)域的細(xì)胞間通訊。在該數(shù)據(jù)集中，觀察到穩(wěn)態(tài)下的CD8 T細(xì)胞與LCMV感染后的CD8 T細(xì)胞之間存在差異表達(dá)。NicheNet可用于觀察淋巴結(jié)中的幾種免疫細(xì)胞群（即單核細(xì)胞，樹突狀細(xì)胞，NK細(xì)胞，B細(xì)胞，CD4 T細(xì)胞）如何調(diào)節(jié)和誘導(dǎo)這些觀察到的基因表達(dá)變化。
#準(zhǔn)備
# devtools::install_github("saeyslab/nichenetr")
library(circlize)
library(nichenetr)
library(Seurat) # please update to Seurat V4
library(tidyverse)
2.1 讀入NicheNet的配體-受體先驗(yàn)?zāi)Ｐ?，配體-受體網(wǎng)絡(luò)和加權(quán)整合網(wǎng)絡(luò)。
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"))
head(lr_network)
## # A tibble: 6 x 4
## from to source database
## <chr> <chr> <chr> <chr>
## 1 CXCL1 CXCR2 kegg_cytokines kegg
## 2 CXCL2 CXCR2 kegg_cytokines kegg
## 3 CXCL3 CXCR2 kegg_cytokines kegg
## 4 CXCL5 CXCR2 kegg_cytokines kegg
## 5 PPBP CXCR2 kegg_cytokines kegg
## 6 CXCL6 CXCR2 kegg_cytokines kegg
weighted_networks = readRDS(url("https://zenodo.org/record/3260758/files/weighted_networks.rds"))
head(weighted_networks$lr_sig) # interactions and their weights in the ligand-receptor + signaling network
## # A tibble: 6 x 3
## from to weight
## <chr> <chr> <dbl>
## 1 A1BG ABCC6 0.422
## 2 A1BG ACE2 0.101
## 3 A1BG ADAM10 0.0970
## 4 A1BG AGO1 0.0525
## 5 A1BG AKT1 0.0855
## 6 A1BG ANXA7 0.457
head(weighted_networks$gr) # interactions and their weights in the gene regulatory network
## # A tibble: 6 x 3
## from to weight
## <chr> <chr> <dbl>
## 1 A1BG A2M 0.0294
## 2 AAAS GFAP 0.0290
## 3 AADAC CYP3A4 0.0422
## 4 AADAC IRF8 0.0275
## 5 AATF ATM 0.0330
## 6 AATF ATR 0.0355
2.2 讀入注釋好的Seurat對(duì)象
seuratObj = readRDS(url("https://zenodo.org/record/3531889/files/seuratObj.rds"))
seuratObj@meta.data %>% head()
## nGene nUMI orig.ident aggregate res.0.6 celltype nCount_RNA nFeature_RNA
## W380370 880 1611 LN_SS SS 1 CD8 T 1607 876
## W380372 541 891 LN_SS SS 0 CD4 T 885 536
## W380374 742 1229 LN_SS SS 0 CD4 T 1223 737
## W380378 847 1546 LN_SS SS 1 CD8 T 1537 838
## W380379 839 1606 LN_SS SS 0 CD4 T 1603 836
## W380381 517 844 LN_SS SS 0 CD4 T 840 513
seuratObj
## An object of class Seurat
## 13541 features across 5027 samples within 1 assay
## Active assay: RNA (13541 features, 1575 variable features)
## 4 dimensional reductions calculated: cca, cca.aligned, tsne, pca
查看存在哪些細(xì)胞群
seuratObj@meta.data$celltype %>% table() # note that the number of cells of some cell types is very low and should preferably be higher for a real application
## .
## B CD4 T CD8 T DC Mono NK Treg
## 382 2562 1645 18 90 131 199
DimPlot(seuratObj, reduction = "tsne")
查看分組
seuratObj@meta.data$aggregate %>% table()
## .
## LCMV SS
## 3886 1141
DimPlot(seuratObj, reduction = "tsne", group.by = "aggregate")
2.3 進(jìn)行NicheNet分析
在這個(gè)演示中，我們希望預(yù)測(cè)哪些配體可能影響CD8 T細(xì)胞在LCMV感染前后的差異表達(dá)基因。因此receiver細(xì)胞群是' CD8 T '細(xì)胞群，而sender細(xì)胞群是' CD4 T '， ' Treg '， ' Mono '， ' NK '， ' B '和' DC '。
我們感興趣的基因是LCMV感染后CD8 T細(xì)胞中差異表達(dá)的基因。因此，將感興趣的條件condition_oi設(shè)置為“LCMV”，而參考/穩(wěn)態(tài)條件condition_reference設(shè)置為“SS”。（計(jì)算差異基因的方法是標(biāo)準(zhǔn)Seurat Wilcoxon檢驗(yàn)）
用于預(yù)測(cè)活性靶基因和構(gòu)建活性配體-受體網(wǎng)絡(luò)的top-ranked配體的數(shù)量默認(rèn)是20個(gè)。（top_n_ligands參數(shù)指定用于后續(xù)分析的高活性配體的數(shù)量 ）
# indicated cell types should be cell class identities
# check via:
# seuratObj %>% Idents() %>% table()
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = seuratObj,
receiver = "CD8 T",
condition_colname = "aggregate", condition_oi = "LCMV", condition_reference = "SS",
sender = c("CD4 T","Treg", "Mono", "NK", "B", "DC"),
ligand_target_matrix = ligand_target_matrix, lr_network = lr_network, weighted_networks = weighted_networks, organism = "mouse")
## [1] "Read in and process NicheNet's networks"
## [1] "Define expressed ligands and receptors in receiver and sender cells"
## [1] "Perform DE analysis in receiver cell"
## [1] "Perform NicheNet ligand activity analysis"
## [1] "Infer active target genes of the prioritized ligands"
## [1] "Infer receptors of the prioritized ligands"
# 輸出的是一個(gè)列表：
nichenet_output %>% names()
## [1] "ligand_activities" "top_ligands" "top_targets"
## [4] "top_receptors" "ligand_target_matrix" "ligand_target_heatmap"
## [7] "ligand_target_df" "ligand_activity_target_heatmap" "ligand_receptor_matrix"
##[10] "ligand_receptor_heatmap" "ligand_receptor_df" "ligand_receptor_matrix_bonafide"
##[13] "ligand_receptor_heatmap_bonafide" "ligand_receptor_df_bonafide" "geneset_oi"
##[16] "background_expressed_genes"
View(nichenet_output)
- 查看配體活性分析結(jié)果
NicheNet做的第一件事是根據(jù)預(yù)測(cè)的配體活性來(lái)確定配體的優(yōu)先級(jí)。使用如下命令查看這些配體的排名:
nichenet_output$ligand_activities
## # A tibble: 44 x 6
## test_ligand auroc aupr pearson rank bona_fide_ligand
## <chr> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 Ebi3 0.662 0.238 0.219 1 FALSE
## 2 Il15 0.596 0.160 0.109 2 TRUE
## 3 Crlf2 0.560 0.160 0.0890 3 FALSE
## 4 App 0.499 0.134 0.0750 4 TRUE
## 5 Tgfb1 0.498 0.134 0.0631 5 TRUE
## 6 Ptprc 0.539 0.142 0.0602 6 TRUE
## 7 H2-M3 0.526 0.149 0.0533 7 TRUE
## 8 Icam1 0.544 0.134 0.0496 8 TRUE
## 9 Cxcl10 0.536 0.134 0.0457 9 TRUE
## 10 Adam17 0.517 0.129 0.0378 10 TRUE
## # ... with 34 more rows
不同的配體活性檢測(cè)值（auroc, aupr,pearson相關(guān)系數(shù)）是用來(lái)評(píng)估配體對(duì)觀測(cè)到的差異表達(dá)基因的預(yù)測(cè)能力。NicheNet主要參考pearson相關(guān)系數(shù)對(duì)這些配體進(jìn)行排序（效果最好）。
‘bona_fide_ligand’這一列的意思是這個(gè)配體是否存在于受體-配體數(shù)據(jù)庫(kù)中（有的話為TRUE）。
查看top20配體
nichenet_output$top_ligands
[1] "Ebi3" "Il15" "Crlf2" "App" "Tgfb1" "Ptprc" "H2-M3"
[8] "Icam1" "Cxcl10" "Adam17" "Cxcl11" "Cxcl9" "H2-T23" "Sema4d"
[15] "Ccl5" "C3" "Cxcl16" "Itgb1" "Anxa1" "Sell"
查看哪個(gè)細(xì)胞群表達(dá)了這些配體
p = DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()
ggsave("top20_ligands.png", p, width = 12, height = 6)
#??%>% rev()這一步是將橫坐標(biāo)的基因反過來(lái)排序
大多數(shù)top-ranked配體似乎主要由樹突狀細(xì)胞和單核細(xì)胞表達(dá)。
觀察這些配體在LCMV感染后是否有有差異表達(dá)
p=DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), split.by = "aggregate") + RotatedAxis()
ggsave("top20_ligands_compare.png", p, width = 12, height = 8)
用小提琴圖對(duì)比配體的表達(dá)情況
p=VlnPlot(seuratObj, features = nichenet_output$top_ligands, split.by = "aggregate", pt.size = 0, combine = T)
ggsave("VlnPlot_ligands_compare.png", p, width = 24, height = 16)
feature后面沒有加%>% rev()，出圖的基因順序就和前兩個(gè)圖相反
- 查看配體調(diào)控靶基因
推斷活躍的配體-靶標(biāo)連接
p = nichenet_output$ligand_target_heatmap
ggsave("Heatmap_ligand-target.png", p, width = 12, height = 6)
這是一個(gè)普通的ggplot2對(duì)象，可以改變顏色和軸標(biāo)簽等等（如下圖）
p = nichenet_output$ligand_target_heatmap + scale_fill_gradient2(low = "whitesmoke", high = "royalblue", breaks = c(0,0.0045,0.009)) + xlab("anti-LCMV response genes in CD8 T cells") + ylab("Prioritized immmune cell ligands")
ggsave("Heatmap_ligand-target2.png", p, width = 12, height = 6)
查看20個(gè) top-ranked配體的top-predicted靶基因。
x = nichenet_output$top_targets
#x2 <- nichenet_output$ligand_target_df
write.csv(x, "ligand_target.csv", row.names = F)
## [1] "Cd274" "Cd53" "Ddit4" "Id3" "Ifit3" "Irf1" "Irf7" "Irf9" "Parp14" "Pdcd4"
## [11] "Pml" "Psmb9" "Rnf213" "Stat1" "Stat2" "Tap1" "Ubc" "Zbp1" "Cd69" "Gbp4"
## [21] "Basp1" "Casp8" "Cxcl10" "Nlrc5" "Vim" "Actb" "Ifih1" "Myh9" "B2m" "H2-T23"
## [31] "Rpl13a" "Cxcr4"
查看這些受體基因在病毒感染前后CD8中的表達(dá)
p = DotPlot(seuratObj %>% subset(idents = "CD8 T"), features = nichenet_output$top_targets %>% rev(), split.by = "aggregate") + RotatedAxis()
ggsave("Targets_Expression_dotplot.png", p, width = 12, height = 6)
查看部分配體調(diào)控靶基因的表達(dá)情況
p=VlnPlot(seuratObj %>% subset(idents = "CD8 T"), features = c("Zbp1","Ifit3","Irf7"), split.by = "aggregate", pt.size = 0, combine = T)
ggsave("Targets_Expression_vlnplot.png", p, width = 12, height = 4)
若設(shè)置combine=F，就會(huì)畫出三張圖
- 查看受體情況
可視化配體和靶基因活性
p = nichenet_output$ligand_activity_target_heatmap
ggsave("Heatmap_ligand_activity_target.png", p, width = 12, height = 6)
查看配體-受體互作
p = nichenet_output$ligand_receptor_heatmap
ggsave("Heatmap_ligand-receptor.png", p, width = 12, height = 6)
x <- nichenet_output$ligand_receptor_matrix
#x <- nichenet_output$ligand_receptor_df
write.csv(x, "ligand_receptor.csv", row.names = F)
查看受體表達(dá)情況
p = DotPlot(seuratObj %>% subset(idents = "CD8 T"),
features = nichenet_output$top_receptors,
split.by = "aggregate") + RotatedAxis()
ggsave("Receptors_Expression_dotplot.png", p, width = 12, height = 6)
p = VlnPlot(seuratObj %>% subset(idents = "CD8 T"), features = nichenet_output$top_receptors,
split.by = "aggregate", pt.size = 0, combine = T, ncol = 8)
ggsave("Receptors_Expression_vlnplot.png", p, width = 12, height = 8)
有文獻(xiàn)報(bào)道的配體-受體
# Show ‘bona fide’ ligand-receptor links that are described in the literature and not predicted based on PPI
p = nichenet_output$ligand_receptor_heatmap_bonafide
ggsave("Heatmap_ligand-receptor_bonafide.png", p, width = 8, height = 4)
x <- nichenet_output$ligand_receptor_matrix_bonafide
#x <- nichenet_output$ligand_receptor_df_bonafide
write.csv(x, "ligand_receptor_bonafide.csv", row.names = F)
3. Circos繪圖來(lái)可視化配體-靶標(biāo)和配體-受體的相互作用。
參考：https://github.com/saeyslab/nichenetr/blob/master/vignettes/circos.md
這一可視化分組根據(jù)最強(qiáng)表達(dá)的細(xì)胞類型預(yù)測(cè)活性配體。因此，我們需要確定每種細(xì)胞類型，它們表達(dá)的配體比其他細(xì)胞類型更強(qiáng)。計(jì)算發(fā)送細(xì)胞中平均配體表達(dá)量。
# avg_expression_ligands = AverageExpression(seuratObj %>% subset(subset = aggregate == "LCMV"),features = nichenet_output$top_ligands) # if want to look specifically in LCMV-only cells
avg_expression_ligands = AverageExpression(seuratObj, features = nichenet_output$top_ligands)
分配配體給發(fā)送細(xì)胞
為了給發(fā)送端細(xì)胞類型分配配體，我們可以查找哪個(gè)發(fā)送端細(xì)胞類型的表達(dá)式高于平均值+ SD。
sender_ligand_assignment = avg_expression_ligands$RNA %>% apply(1, function(ligand_expression){
ligand_expression > (ligand_expression %>% mean() + ligand_expression %>% sd())
}) %>% t()
sender_ligand_assignment[1:4,1:4]
# CD8 T CD4 T Treg B
# Ebi3 FALSE FALSE FALSE FALSE
# Il15 FALSE FALSE FALSE FALSE
# Crlf2 FALSE FALSE FALSE FALSE
# App FALSE FALSE FALSE FALSE
sender_ligand_assignment = sender_ligand_assignment %>% apply(2, function(x){x[x == TRUE]}) %>% purrr::keep(function(x){length(x) > 0})
names(sender_ligand_assignment)
## [1] "B" "NK" "Mono" "DC"
(sender_ligand_assignment)
# $B
# H2-M3
# TRUE
# $NK
# Ptprc Itgb1
# TRUE TRUE
# $Mono
# Ebi3 Crlf2 App Tgfb1 Cxcl10 Adam17 Cxcl11 Cxcl9 Sema4d C3 Anxa1
# TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# $DC
# Il15 Icam1 H2-T23 Ccl5 Cxcl16 Itgb1
# TRUE TRUE TRUE TRUE TRUE TRUE
頂部的配體似乎在B細(xì)胞、NK細(xì)胞、單核細(xì)胞和DCs中表達(dá)最強(qiáng)烈。我們也會(huì)知道在多種細(xì)胞類型中哪些配體是常見的(=特定于> 1細(xì)胞類型的配體，或在之前的代碼塊中未指定給某個(gè)細(xì)胞類型的配體)?，F(xiàn)在確定哪些優(yōu)先配體是由CAFs或內(nèi)皮細(xì)胞表達(dá)的。
all_assigned_ligands = sender_ligand_assignment %>% lapply(function(x){names(x)}) %>% unlist()
unique_ligands = all_assigned_ligands %>% table() %>% .[. == 1] %>% names()
general_ligands = nichenet_output$top_ligands %>% setdiff(unique_ligands)
B_specific_ligands = sender_ligand_assignment$B %>% names() %>% setdiff(general_ligands)
NK_specific_ligands = sender_ligand_assignment$NK %>% names() %>% setdiff(general_ligands)
Mono_specific_ligands = sender_ligand_assignment$Mono %>% names() %>% setdiff(general_ligands)
DC_specific_ligands = sender_ligand_assignment$DC %>% names() %>% setdiff(general_ligands)
ligand_type_indication_df = tibble(
ligand_type = c(rep("B-specific", times = B_specific_ligands %>% length()),
rep("NK-specific", times = NK_specific_ligands %>% length()),
rep("Mono-specific", times = Mono_specific_ligands %>% length()),
rep("DC-specific", times = DC_specific_ligands %>% length()),
rep("General", times = general_ligands %>% length())),
ligand = c(B_specific_ligands, NK_specific_ligands, Mono_specific_ligands, DC_specific_ligands, general_ligands))
ligand_type_indication_df %>% head
## A tibble: 6 x 2
# ligand_type ligand
# <chr> <chr>
#1 B-specific H2-M3
#2 NK-specific Ptprc
#3 Mono-specific Ebi3
#4 Mono-specific Crlf2
#5 Mono-specific App
#6 Mono-specific Tgfb1
定義感興趣的配體-目標(biāo)鏈接
為了避免circos圖中有太多配體目標(biāo)鏈接，我們將只顯示權(quán)重高于預(yù)定義截止值的鏈接:屬于最低分?jǐn)?shù)的40%的鏈接被刪除。這并不是說用于這種可視化的邊界和其他邊界可以根據(jù)用戶的需要進(jìn)行更改。
active_ligand_target_links_df = nichenet_output$ligand_target_df %>% mutate(target_type = "LCMV-DE") %>% inner_join(ligand_type_indication_df) # if you want ot make circos plots for multiple gene sets, combine the different data frames and differentiate which target belongs to which gene set via the target type
cutoff_include_all_ligands = active_ligand_target_links_df$weight %>% quantile(0.40)
active_ligand_target_links_df_circos = active_ligand_target_links_df %>% filter(weight > cutoff_include_all_ligands)
ligands_to_remove = setdiff(active_ligand_target_links_df$ligand %>% unique(), active_ligand_target_links_df_circos$ligand %>% unique())
targets_to_remove = setdiff(active_ligand_target_links_df$target %>% unique(), active_ligand_target_links_df_circos$target %>% unique())
circos_links = active_ligand_target_links_df %>% filter(!target %in% targets_to_remove &!ligand %in% ligands_to_remove)
circos_links
## A tibble: 125 x 5
# ligand target weight target_type ligand_type
# <chr> <chr> <dbl> <chr> <chr>
# 1 Ebi3 Cd274 0.00325 LCMV-DE Mono-specific
# 2 Ebi3 Cd53 0.00321 LCMV-DE Mono-specific
# 3 Ebi3 Ddit4 0.00335 LCMV-DE Mono-specific
# 4 Ebi3 Id3 0.00373 LCMV-DE Mono-specific
# 5 Ebi3 Ifit3 0.00320 LCMV-DE Mono-specific
# 6 Ebi3 Irf1 0.00692 LCMV-DE Mono-specific
# 7 Ebi3 Irf7 0.00312 LCMV-DE Mono-specific
# 8 Ebi3 Irf9 0.00543 LCMV-DE Mono-specific
# 9 Ebi3 Parp14 0.00336 LCMV-DE Mono-specific
#10 Ebi3 Pdcd4 0.00335 LCMV-DE Mono-specific
## … with 115 more rows
準(zhǔn)備circos可視化:給每個(gè)片段配體和目標(biāo)特定的顏色和順序
grid_col_ligand =c("General" = "lawngreen",
"NK-specific" = "royalblue",
"B-specific" = "darkgreen",
"Mono-specific" = "violet",
"DC-specific" = "steelblue2")
grid_col_target =c(
"LCMV-DE" = "tomato")
grid_col_tbl_ligand = tibble(ligand_type = grid_col_ligand %>% names(), color_ligand_type = grid_col_ligand)
grid_col_tbl_target = tibble(target_type = grid_col_target %>% names(), color_target_type = grid_col_target)
circos_links = circos_links %>% mutate(ligand = paste(ligand," ")) # extra space: make a difference between a gene as ligand and a gene as target!
circos_links = circos_links %>% inner_join(grid_col_tbl_ligand) %>% inner_join(grid_col_tbl_target)
links_circle = circos_links %>% select(ligand,target, weight)
ligand_color = circos_links %>% distinct(ligand,color_ligand_type)
grid_ligand_color = ligand_color$color_ligand_type %>% set_names(ligand_color$ligand)
target_color = circos_links %>% distinct(target,color_target_type)
grid_target_color = target_color$color_target_type %>% set_names(target_color$target)
grid_col =c(grid_ligand_color,grid_target_color)
# give the option that links in the circos plot will be transparant ~ ligand-target potential score
transparency = circos_links %>% mutate(weight =(weight-min(weight))/(max(weight)-min(weight))) %>% mutate(transparency = 1-weight) %>% .$transparency
準(zhǔn)備可視化的circos:排序配體和目標(biāo)
target_order = circos_links$target %>% unique()
ligand_order = c(Mono_specific_ligands, DC_specific_ligands, NK_specific_ligands,B_specific_ligands, general_ligands) %>% c(paste(.," ")) %>% intersect(circos_links$ligand)
order = c(ligand_order,target_order)
準(zhǔn)備circos可視化:定義不同片段之間的間隙
width_same_cell_same_ligand_type = 0.5
width_different_cell = 6
width_ligand_target = 15
width_same_cell_same_target_type = 0.5
gaps = c(
# width_ligand_target,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "Mono-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "DC-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "NK-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "B-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "General") %>% distinct(ligand) %>% nrow() -1)),
width_ligand_target,
rep(width_same_cell_same_target_type, times = (circos_links %>% filter(target_type == "LCMV-DE") %>% distinct(target) %>% nrow() -1)),
width_ligand_target
)
渲染circos的情節(jié)(所有鏈接相同的透明度)。只有表明每個(gè)靶基因的阻滯的寬度與配體-靶的調(diào)控電位成正比(~支持調(diào)控相互作用的先驗(yàn)知識(shí))。
circos.par(gap.degree = gaps)
chordDiagram(links_circle, directional = 1,order=order,link.sort = TRUE, link.decreasing = FALSE, grid.col = grid_col,transparency = 0, diffHeight = 0.005, direction.type = c("diffHeight", "arrows"),link.arr.type = "big.arrow", link.visible = links_circle$weight >= cutoff_include_all_ligands,annotationTrack = "grid",
preAllocateTracks = list(track.height = 0.075))
# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 1)
}, bg.border = NA) #
circos.clear()
繪制circos圖(透明度由配體-靶標(biāo)相互作用的調(diào)控潛力決定)
circos.par(gap.degree = gaps)
chordDiagram(links_circle, directional = 1,order=order,link.sort = TRUE, link.decreasing = FALSE, grid.col = grid_col,transparency = transparency, diffHeight = 0.005, direction.type = c("diffHeight", "arrows"),link.arr.type = "big.arrow", link.visible = links_circle$weight >= cutoff_include_all_ligands,annotationTrack = "grid",
preAllocateTracks = list(track.height = 0.075))
# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 1)
}, bg.border = NA) #
circos.clear()
svg("ligand_target_circos.svg", width = 10, height = 10)
circos.par(gap.degree = gaps)
chordDiagram(links_circle, directional = 1,order=order,link.sort = TRUE, link.decreasing = FALSE, grid.col = grid_col,transparency = transparency, diffHeight = 0.005, direction.type = c("diffHeight", "arrows"),link.arr.type = "big.arrow", link.visible = links_circle$weight >= cutoff_include_all_ligands,annotationTrack = "grid",
preAllocateTracks = list(track.height = 0.075))
# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 1)
}, bg.border = NA) #
circos.clear()
dev.off()
在circos圖中可視化優(yōu)先配體與受體的相互作用
lr_network_top_df = nichenet_output$ligand_receptor_df %>% mutate(receptor_type = "LCMV_CD8T_receptor") %>% inner_join(ligand_type_indication_df)
grid_col_ligand =c("General" = "lawngreen",
"NK-specific" = "royalblue",
"B-specific" = "darkgreen",
"Mono-specific" = "violet",
"DC-specific" = "steelblue2")
grid_col_receptor =c(
"LCMV_CD8T_receptor" = "darkred")
grid_col_tbl_ligand = tibble(ligand_type = grid_col_ligand %>% names(), color_ligand_type = grid_col_ligand)
grid_col_tbl_receptor = tibble(receptor_type = grid_col_receptor %>% names(), color_receptor_type = grid_col_receptor)
circos_links = lr_network_top_df %>% mutate(ligand = paste(ligand," ")) # extra space: make a difference between a gene as ligand and a gene as receptor!
circos_links = circos_links %>% inner_join(grid_col_tbl_ligand) %>% inner_join(grid_col_tbl_receptor)
links_circle = circos_links %>% select(ligand,receptor, weight)
ligand_color = circos_links %>% distinct(ligand,color_ligand_type)
grid_ligand_color = ligand_color$color_ligand_type %>% set_names(ligand_color$ligand)
receptor_color = circos_links %>% distinct(receptor,color_receptor_type)
grid_receptor_color = receptor_color$color_receptor_type %>% set_names(receptor_color$receptor)
grid_col =c(grid_ligand_color,grid_receptor_color)
# give the option that links in the circos plot will be transparant ~ ligand-receptor potential score
transparency = circos_links %>% mutate(weight =(weight-min(weight))/(max(weight)-min(weight))) %>% mutate(transparency = 1-weight) %>% .$transparency
制備可視化的circos:有序配體和受體
receptor_order = circos_links$receptor %>% unique()
ligand_order = c(Mono_specific_ligands, DC_specific_ligands, NK_specific_ligands,B_specific_ligands, general_ligands) %>% c(paste(.," ")) %>% intersect(circos_links$ligand)
order = c(ligand_order,receptor_order)
準(zhǔn)備馬戲團(tuán)可視化:定義不同片段之間的間隙
width_same_cell_same_ligand_type = 0.5
width_different_cell = 6
width_ligand_receptor = 15
width_same_cell_same_receptor_type = 0.5
gaps = c(
# width_ligand_target,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "Mono-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "DC-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "NK-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "B-specific") %>% distinct(ligand) %>% nrow() -1)),
width_different_cell,
rep(width_same_cell_same_ligand_type, times = (circos_links %>% filter(ligand_type == "General") %>% distinct(ligand) %>% nrow() -1)),
width_ligand_receptor,
rep(width_same_cell_same_receptor_type, times = (circos_links %>% filter(receptor_type == "LCMV_CD8T_receptor") %>% distinct(receptor) %>% nrow() -1)),
width_ligand_receptor
)
渲染馬戲團(tuán)的情節(jié)(所有鏈接相同的透明度)。只有表明每個(gè)受體的阻滯的寬度與配體-受體相互作用的重量成比例(~支持相互作用的先驗(yàn)知識(shí))。
circos.par(gap.degree = gaps)
chordDiagram(links_circle, directional = 1,order=order,link.sort = TRUE, link.decreasing = FALSE, grid.col = grid_col,transparency = 0, diffHeight = 0.005, direction.type = c("diffHeight", "arrows"),link.arr.type = "big.arrow", link.visible = links_circle$weight >= cutoff_include_all_ligands,annotationTrack = "grid",
preAllocateTracks = list(track.height = 0.075))
# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 0.8)
}, bg.border = NA) #
circos.clear()
渲染circos圖(透明程度由配體-受體相互作用的先驗(yàn)相互作用權(quán)重決定——正如指示每個(gè)受體的塊的寬度)
circos.par(gap.degree = gaps)
chordDiagram(links_circle, directional = 1,order=order,link.sort = TRUE, link.decreasing = FALSE, grid.col = grid_col,transparency = transparency, diffHeight = 0.005, direction.type = c("diffHeight", "arrows"),link.arr.type = "big.arrow", link.visible = links_circle$weight >= cutoff_include_all_ligands,annotationTrack = "grid",
preAllocateTracks = list(track.height = 0.075))
# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.55), cex = 0.8)
}, bg.border = NA) #
circos.clear()
