常用的細胞通訊軟件：

• CellphoneDB：是公開的人工校正的铅鲤，儲存受體身坐、配體以及兩種相互作用的數據庫巴帮。此外晶密，還考慮了結構組成，能夠描述異構復合物秀鞭。（配體-受體+多聚體）
• iTALK：通過平均表達量方式趋观，篩選高表達的胚體和受體，根據結果作圈圖锋边。（配體-受體）
• CellChat：CellChat將細胞的基因表達數據作為輸入皱坛，并結合配體受體及其輔助因子的相互作用來模擬細胞間通訊咱枉。（配體-受體+多聚體+輔因子）
• NicheNet // NicheNet多樣本分析 // 目標基因的配體和靶基因活性預測：通過將相互作用細胞的表達數據與信號和基因調控網絡的先驗知識相結合來預測相互作用細胞之間的配體-靶標聯系的方法彤路。( 配體-受體+信號通路）
  附：NicheNet使用的常見問題匯總

其它細胞互作軟件還包括Celltalker，SingleCellSignalR躯保，scTensor和SoptSC（這幾個也是基于配體-受體相互作用）

官網：https://github.com/saeyslab/nichenetr

Nature Methods原文

1. NicheNet介紹

1.1 NicheNet工作流程：

一般的預測細胞交互的軟件往往只考慮sender細胞的配體和receiver細胞的受體表達往扔，但細胞交互過程除了配體-受體相互作用以外贩猎，還包含了receiver細胞接受信號后相關通路的激活。
NicheNet輸入基因表達數據瓤球，并將其與通過整合信號通路而構建的模型相結合融欧。不止是預測配體與受體的相互作用敏弃，還整合了細胞內信號傳導卦羡。因此，NicheNet可以預測1)來自一或多種細胞中的配體（sender）影響了與之相互作用的細胞中哪些基因的表達和2)哪些靶基因受每種配體影響以及可能涉及哪些信號傳導介質。

NicheNet工作流程圖

首先從公共數據庫中收集配體-受體配對信息绿饵、信號通路欠肾、基因調控網絡等信息，整合成配體主導的權重配體-靶基因調控模型拟赊。然后將可能受到細胞通訊影響的差異基因集輸入先驗模型刺桃，可以計算與這些基因相關的配體的相關性系數。最后挑選根據相關性系數排行靠前的配體吸祟，依據先驗數據推測與之匹配的受體瑟慈、靶基因及下游信號網絡等信息。

NicheNet不同于大多數研究細胞間通訊的方法屋匕，它著眼于配體對下游基因調控作用葛碧。NicheNet可以預測哪些配體影響另一個細胞中的表達，哪些靶基因受到配體的影響以及哪些信號傳導可能參與其中过吻。

1.2 NicheNet工作pipeline

• 在單細胞數據中定義一個“sender/niche”細胞群（配體細胞）和一個“receiver/target”細胞群（受體細胞）并確定這兩個細胞群都表達哪些基因进泼。
• 定義一個感興趣的基因集：這些基因來自受體細胞群，是可能受到與其相互作用的細胞配體調控的基因集纤虽。（例如：case-control中的差異表達基因乳绕，也可以是細胞的signature或其他基因集）
• 定義一個潛在的配體集：這些配體由配體細胞群中高表達（如10%以上的細胞表達）并可以與受體細胞群表達的受體相結合（通過先驗數據推斷）。
• 進行NicheNet配體活性分析：其活性主要通過配體與受體細胞中的差異基因集的相關性進行判斷
• 推斷在配體活性分析中的top-ranked配體所調控的top-predicted靶基因逼纸，以及與配體配對的受體洋措。

NicheNet提供了一個三個功能相似的打包函數： nichenet_seuratobj_aggregate, nichenet_seuratobj_cluster_de and nichenet_seuratobj_aggregate_cluster_de.它們可以一步完成上述五步seurat對象的配體調控網絡分析。

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

本文的演示數據集和代碼來自NicheNet官方分析單細胞數據的教程：https://github.com/saeyslab/nichenetr/blob/master/vignettes/seurat_wrapper.md

我們將使用Medaglia等人的小鼠NICHE-seq數據杰刽，探索淋巴細胞性脈絡膜腦膜炎病毒（LCMV）感染之前和之后72小時的腹股溝淋巴結T細胞區(qū)域的細胞間通訊呻纹。在該數據集中，觀察到穩(wěn)態(tài)下的CD8 T細胞與LCMV感染后的CD8 T細胞之間存在差異表達专缠。NicheNet可用于觀察淋巴結中的幾種免疫細胞群（即單核細胞雷酪，樹突狀細胞，NK細胞涝婉，B細胞哥力，CD4 T細胞）如何調節(jié)和誘導這些觀察到的基因表達變化。

#準備
# devtools::install_github("saeyslab/nichenetr")
library(circlize)
library(nichenetr)
library(Seurat) # please update to Seurat V4
library(tidyverse)

2.1 讀入NicheNet的配體-受體先驗模型墩弯，配體-受體網絡和加權整合網絡吩跋。

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對象

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

查看存在哪些細胞群

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 進行NicheNet分析

在這個演示中，我們希望預測哪些配體可能影響CD8 T細胞在LCMV感染前后的差異表達基因渔工。因此receiver細胞群是' CD8 T '細胞群锌钮，而sender細胞群是' CD4 T '， ' Treg '引矩， ' Mono '梁丘， ' NK '侵浸， ' B '和' DC '。
我們感興趣的基因是LCMV感染后CD8 T細胞中差異表達的基因氛谜。因此掏觉，將感興趣的條件condition_oi設置為“LCMV”，而參考/穩(wěn)態(tài)條件condition_reference設置為“SS”值漫。（計算差異基因的方法是標準Seurat Wilcoxon檢驗）
用于預測活性靶基因和構建活性配體-受體網絡的top-ranked配體的數量默認是20個澳腹。（top_n_ligands參數指定用于后續(xù)分析的高活性配體的數量 ）

# 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"

# 輸出的是一個列表：
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)

• 查看配體活性分析結果
  NicheNet做的第一件事是根據預測的配體活性來確定配體的優(yōu)先級。使用如下命令查看這些配體的排名:

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

不同的配體活性檢測值（auroc, aupr,pearson相關系數）是用來評估配體對觀測到的差異表達基因的預測能力杨何。NicheNet主要參考pearson相關系數對這些配體進行排序（效果最好）酱塔。
‘bona_fide_ligand’這一列的意思是這個配體是否存在于受體-配體數據庫中（有的話為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"  

查看哪個細胞群表達了這些配體

p = DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()
ggsave("top20_ligands.png", p, width = 12, height = 6)
#??%>% rev()這一步是將橫坐標的基因反過來排序

大多數top-ranked配體似乎主要由樹突狀細胞和單核細胞表達危虱。

觀察這些配體在LCMV感染后是否有有差異表達

p=DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), split.by = "aggregate") + RotatedAxis()
ggsave("top20_ligands_compare.png", p, width = 12, height = 8)

用小提琴圖對比配體的表達情況

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()延旧，出圖的基因順序就和前兩個圖相反

• 查看配體調控靶基因

推斷活躍的配體-靶標連接

p = nichenet_output$ligand_target_heatmap
ggsave("Heatmap_ligand-target.png", p, width = 12, height = 6)

這是一個普通的ggplot2對象，可以改變顏色和軸標簽等等（如下圖）

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個 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中的表達

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)

查看部分配體調控靶基因的表達情況

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)

若設置combine=F迁沫，就會畫出三張圖

• 查看受體情況

可視化配體和靶基因活性

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)

查看受體表達情況

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)

有文獻報道的配體-受體

# 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繪圖來可視化配體-靶標和配體-受體的相互作用。

參考：https://github.com/saeyslab/nichenetr/blob/master/vignettes/circos.md
這一可視化分組根據最強表達的細胞類型預測活性配體捌蚊。因此集畅，我們需要確定每種細胞類型，它們表達的配體比其他細胞類型更強缅糟。計算發(fā)送細胞中平均配體表達量挺智。

# 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ā)送細胞
為了給發(fā)送端細胞類型分配配體，我們可以查找哪個發(fā)送端細胞類型的表達式高于平均值+ 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細胞赦颇、NK細胞、單核細胞和DCs中表達最強烈赴涵。我們也會知道在多種細胞類型中哪些配體是常見的(=特定于> 1細胞類型的配體媒怯，或在之前的代碼塊中未指定給某個細胞類型的配體)。現在確定哪些優(yōu)先配體是由CAFs或內皮細胞表達的髓窜。

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 

定義感興趣的配體-目標鏈接
為了避免circos圖中有太多配體目標鏈接扇苞，我們將只顯示權重高于預定義截止值的鏈接:屬于最低分數的40%的鏈接被刪除。這并不是說用于這種可視化的邊界和其他邊界可以根據用戶的需要進行更改寄纵。

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

準備circos可視化:給每個片段配體和目標特定的顏色和順序

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 

準備可視化的circos:排序配體和目標

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)

準備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é)(所有鏈接相同的透明度)鳖敷。只有表明每個靶基因的阻滯的寬度與配體-靶的調控電位成正比(~支持調控相互作用的先驗知識)。

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圖(透明度由配體-靶標相互作用的調控潛力決定)

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)

準備馬戲團可視化:定義不同片段之間的間隙

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
  )

渲染馬戲團的情節(jié)(所有鏈接相同的透明度)程拭。只有表明每個受體的阻滯的寬度與配體-受體相互作用的重量成比例(~支持相互作用的先驗知識)定踱。

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圖(透明程度由配體-受體相互作用的先驗相互作用權重決定——正如指示每個受體的塊的寬度)

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()

參考：http://www.reibang.com/p/46f487f58a9c