上周一位讀者聯(lián)系我但两，讓我?guī)兔Πl(fā)一個(gè)繪圖的單子。在朋友圈發(fā)單后供置，感興趣的朋友很多镜遣，有十幾位還私聊我讓我分享一下代碼，可見(jiàn)大家還是很感興趣的士袄。不過(guò)等了兩天悲关，依舊沒(méi)有勇士接單，可能是因?yàn)檫@種圖比較少見(jiàn)娄柳，大家畫(huà)得少寓辱。
先來(lái)理解一下這張圖，在b圖中：
- 左邊是
EC細(xì)胞表達(dá)的ligand赤拒，右邊是mNEUR細(xì)胞表達(dá)的receptor秫筏。 -
ligand這一列對(duì)應(yīng)的基因會(huì)排序，依據(jù)是兩個(gè)group(比如young和old兩組)在EC細(xì)胞中找完差異基因后挎挖，能知道這些基因的log(FC_young_to_old)这敬，從大到小依次往下排，并涂色蕉朵。同時(shí)崔涂，找差異基因也能知道p值，原文表示顯著性用的是圓環(huán)的顏色始衅，越顯著冷蚂，圓環(huán)顏色越深。 -
receptor這一列類(lèi)似ligand列汛闸，不過(guò)是在mNEUR細(xì)胞中蝙茶，兩組之間找差異基因。 - 圖中的線(xiàn)段诸老，表示
ligand-receptor對(duì)隆夯，而這個(gè)信息（誰(shuí)和誰(shuí)構(gòu)成受配體）是已知的，數(shù)據(jù)庫(kù)可以下載到（也就是說(shuō)，這個(gè)圖即便不做細(xì)胞通訊分析蹄衷，只要能下載到ligand-receptor信息苞尝，也能畫(huà)）。cellphonedb的輸出結(jié)果也間接含有這個(gè)信息宦芦，我的代碼基于這個(gè)。 - 至于連線(xiàn)的顏色轴脐，我看了一下调卑，沒(méi)怎么變化，所以我推測(cè)線(xiàn)段的顏色取決于線(xiàn)段兩側(cè)ligand-receptor是高表達(dá)還是低表達(dá)大咱，高表達(dá)時(shí)線(xiàn)段用淺紅色恬涧，低表達(dá)時(shí)線(xiàn)段用淺藍(lán)色。
以上解讀均為我的理解碴巾，我沒(méi)有看原文溯捆。我的圖會(huì)略作更改，下文會(huì)說(shuō)厦瓢。
另外提揍，因?yàn)橐恍┗プ麝P(guān)系沒(méi)法用一個(gè)ligand一個(gè)receptor去表示，這樣的pair不適合用這個(gè)圖展示煮仇。
此篇推文的代碼一共有3個(gè)：CCC_compare2.R劳跃、CCC_line.R、CCC_line2.R浙垫。每個(gè)代碼都能出一張有意義的圖刨仑。
準(zhǔn)備數(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)生的細(xì)胞|receptor產(chǎn)生的細(xì)胞
plot.width=15,plot.height=30,filename = "test3_"
)
前幾行代碼跟上一節(jié)是類(lèi)似的，也能輸出一個(gè)表格（test3_Old2Young.xlsx）和對(duì)應(yīng)的氣泡圖夹姥。
后續(xù)的連線(xiàn)圖基于這個(gè)表格的數(shù)據(jù)杉武。
連線(xiàn)圖有個(gè)特點(diǎn)，就是"誰(shuí)是ligand辙售，誰(shuí)是receptor轻抱，ligand和receptor分別由什么細(xì)胞產(chǎn)生"非常明確，CCC_compare2.R這個(gè)代碼比CCC_compare.R (在上一節(jié))增加一百多行去理清這個(gè)事情旦部，輸出的圖和表格pair左邊的就是ligand(產(chǎn)生的細(xì)胞)十拣，右邊的就是receptor(產(chǎn)生的細(xì)胞)
連線(xiàn)圖的繪制 第1種方法
source("CCC_line.R")
CCC_line(table.path="test3_Old2Young.xlsx",ligand.cell="EC",receptor.cell="APC",
group1.name = "Old",group2.name = "Young",#這五個(gè)參數(shù)和上一步對(duì)應(yīng)
ligand.color="#4dbbd6",receptor.color="#90d1c1",
pt.size=6,
line.thre1=0.5,line.thre2=6,#line.thre1和上一步的"thre"參數(shù)一致，line.thre2可以用來(lái)調(diào)整線(xiàn)的粗細(xì)志鹃，值越大际起，線(xiàn)越細(xì)
file.name="test3_",plot.width=25,plot.height=20)
第一種方法不涉及差異基因，因此左右兩列是統(tǒng)一的圓點(diǎn)缝彬。線(xiàn)段粗顏色紅财忽，表示（相較于group2）group1的互作強(qiáng)；線(xiàn)段粗顏色藍(lán)，表示（相較于group2）group1的互作弱秘血。
連線(xiàn)圖的繪制 第2種方法
library(Seurat)
testseu=readRDS("testseu.rds")
# 此次演示為了加快運(yùn)行速度味抖，人為減少了數(shù)據(jù)量，實(shí)際分析中找差異基因不建議這么做
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ù)集灰粮，在運(yùn)行cellphonedb時(shí)仔涩，進(jìn)行了基因symbol的轉(zhuǎn)換。
#此處找差異基因得到的symbol為真實(shí)基因名粘舟，為了讓兩個(gè)分析匹配熔脂，DEG表格也應(yīng)該做基因名轉(zhuǎn)換。
#但是為了簡(jiǎn)化柑肴，此處只是簡(jiǎn)單地將小鼠基因名轉(zhuǎn)為大寫(xiě)霞揉，不是很精確。大家在分析的時(shí)候建議嚴(yán)格一點(diǎ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
)
這種方法得到的圖，左右兩列添加了差異基因相關(guān)的信息硕舆，p值秽荞、log2FC。線(xiàn)段的粗細(xì)固定了抚官，線(xiàn)段的顏色表示（相較于group2）group1的互作強(qiáng)弱蚂会。
本文代碼編寫(xiě)(3個(gè)函數(shù))花費(fèi)大量時(shí)間，故不無(wú)償提供耗式，有需要的朋友可以后-臺(tái)回復(fù)2022B
