上周一位讀者聯(lián)系我，讓我?guī)兔Πl(fā)一個繪圖的單子。在朋友圈發(fā)單后，感興趣的朋友很多，有十幾位還私聊我讓我分享一下代碼，可見大家還是很感興趣的。不過等了兩天，依舊沒有勇士接單，可能是因為這種圖比較少見，大家畫得少。

image

先來理解一下這張圖，在b圖中：

• 左邊是EC細胞表達的ligand，右邊是mNEUR細胞表達的receptor。
• ligand這一列對應(yīng)的基因會排序，依據(jù)是兩個group(比如young和old兩組)在EC細胞中找完差異基因后，能知道這些基因的log(FC_young_to_old)，從大到小依次往下排，并涂色。同時，找差異基因也能知道p值，原文表示顯著性用的是圓環(huán)的顏色，越顯著，圓環(huán)顏色越深。
• receptor這一列類似ligand列，不過是在mNEUR細胞中，兩組之間找差異基因。
• 圖中的線段，表示ligand-receptor對，而這個信息（誰和誰構(gòu)成受配體）是已知的，數(shù)據(jù)庫可以下載到（也就是說，這個圖即便不做細胞通訊分析，只要能下載到ligand-receptor信息，也能畫）。cellphonedb的輸出結(jié)果也間接含有這個信息，我的代碼基于這個。
• 至于連線的顏色，我看了一下，沒怎么變化，所以我推測線段的顏色取決于線段兩側(cè)ligand-receptor是高表達還是低表達，高表達時線段用淺紅色，低表達時線段用淺藍色。

以上解讀均為我的理解，我沒有看原文。我的圖會略作更改，下文會說。

另外，因為一些互作關(guān)系沒法用一個ligand一個receptor去表示，這樣的pair不適合用這個圖展示。

此篇推文的代碼一共有3個：CCC_compare2.R、CCC_line.R、CCC_line2.R。每個代碼都能出一張有意義的圖。

準備數(shù)據(jù)

source("CCC_compare2.R")
CCC_compare2(group1.name = "Old",group2.name = "Young",
             group1.pfile = "cellphonedb/Old/pvalues.txt",group1.mfile="cellphonedb/Old/means.txt",
             group2.pfile="cellphonedb/Young/pvalues.txt",group2.mfile="cellphonedb/Young/means.txt",
             p.threshold = 0.05,thre=0.5,
             cell.pair="EC|APC", #指定ligand產(chǎn)生的細胞|receptor產(chǎn)生的細胞
             plot.width=15,plot.height=30,filename = "test3_"
)

前幾行代碼跟上一節(jié)是類似的，也能輸出一個表格（test3_Old2Young.xlsx）和對應(yīng)的氣泡圖。
后續(xù)的連線圖基于這個表格的數(shù)據(jù)。

連線圖有個特點，就是"誰是ligand，誰是receptor，ligand和receptor分別由什么細胞產(chǎn)生"非常明確，CCC_compare2.R這個代碼比CCC_compare.R (在上一節(jié))增加一百多行去理清這個事情，輸出的圖和表格pair左邊的就是ligand(產(chǎn)生的細胞)，右邊的就是receptor(產(chǎn)生的細胞)

image

連線圖的繪制 第1種方法

source("CCC_line.R")
CCC_line(table.path="test3_Old2Young.xlsx",ligand.cell="EC",receptor.cell="APC",
         group1.name = "Old",group2.name = "Young",#這五個參數(shù)和上一步對應(yīng)
         ligand.color="#4dbbd6",receptor.color="#90d1c1",
         pt.size=6,
         line.thre1=0.5,line.thre2=6,#line.thre1和上一步的"thre"參數(shù)一致，line.thre2可以用來調(diào)整線的粗細，值越大，線越細
         file.name="test3_",plot.width=25,plot.height=20)

image

第一種方法不涉及差異基因，因此左右兩列是統(tǒng)一的圓點。線段粗顏色紅，表示（相較于group2）group1的互作強；線段粗顏色藍，表示（相較于group2）group1的互作弱。

連線圖的繪制 第2種方法

library(Seurat)
testseu=readRDS("testseu.rds")
# 此次演示為了加快運行速度，人為減少了數(shù)據(jù)量，實際分析中找差異基因不建議這么做
selectedCB=sample(testseu@meta.data$CB,1000)
testseu=testseu%>%subset(CB %in% selectedCB)

# 基于分組找差異基因
marker_group=data.frame()
Idents(testseu)="celltype_age"
for ( ci in c("EC","APC") ) {
  tmp.marker <- FindMarkers(
    testseu, logfc.threshold = 0, min.pct = 0.01,
    only.pos = F, test.use = "wilcox",
    ident.1=paste0(ci,"_Old"),ident.2=paste0(ci,"_Young")
  )
  
  tmp.marker$gene=rownames(tmp.marker)
  tmp.marker$cluster_group=ifelse(tmp.marker$avg_log2FC > 0,paste0(ci,"_Old"),paste0(ci,"_Young"))
  tmp.marker$cluster=ci
  tmp.marker=tmp.marker%>%arrange(desc(avg_log2FC))
  
  marker_group=marker_group%>%rbind(tmp.marker)
}
#本次演示的數(shù)據(jù)集為小鼠數(shù)據(jù)集，在運行cellphonedb時，進行了基因symbol的轉(zhuǎn)換。
#此處找差異基因得到的symbol為真實基因名，為了讓兩個分析匹配，DEG表格也應(yīng)該做基因名轉(zhuǎn)換。
#但是為了簡化，此處只是簡單地將小鼠基因名轉(zhuǎn)為大寫，不是很精確。大家在分析的時候建議嚴格一點。
marker_group$gene=marker_group$gene %>% toupper()

第二種方法需要不設(shè)置篩選條件去找差異基因，然后用CCC_line2函數(shù)就可以了

source("CCC_line2.R")
CCC_line2(cpdb.table.path = "test3_Old2Young.xlsx",marker_group = marker_group,
          ligand.cell = "EC",receptor.cell = "APC",
          group1.name = "Old",group2.name = "Young",
          line.size = 2,file.name = "test3_",plot.width = 25,plot.height = 20
          )

image

這種方法得到的圖，左右兩列添加了差異基因相關(guān)的信息，p值、log2FC。線段的粗細固定了，線段的顏色表示（相較于group2）group1的互作強弱。

本文代碼編寫(3個函數(shù))花費大量時間，故不無償提供，有需要的朋友可以后-臺回復(fù)2022B