本期內(nèi)容已同步分享至單細胞天地

細胞通訊是指細胞與細胞之間的聯(lián)系趟妥。細胞和人類一樣钝侠，多細胞生物的很多細胞會相互作用，形成“細胞社會”伐弹，在這個社會里拉馋，細胞與細胞之間會發(fā)生相互作用和信息的傳遞，細胞建立通訊聯(lián)絡是必需的。如生物體的生長發(fā)育煌茴、分化柠逞、各種組織器官的形成、組織的維持以及它們各種生理活動的協(xié)調(diào)景馁。經(jīng)典的例子莫過于 神經(jīng)細胞之間的神經(jīng)遞質(zhì)的傳遞與接收。

細胞與細胞之間的通訊有三種形式：

• 細胞之間的直接接觸逗鸣；
• 細胞之間通過其間的胞外基質(zhì)相互聯(lián)系合住；
• 細胞之間通過分子間的相互作用產(chǎn)生聯(lián)系。

在這個過程中撒璧，“受體-配體”的概念十分重要透葛，什么是受體，什么是配體卿樱？受體與配體之間結(jié)合的結(jié)果是受體被激活僚害，并產(chǎn)生受體激活后續(xù)信號傳遞的基本步驟。在單細胞分析當中繁调，不同細胞類型之間的通訊可能會對某些生物學過程具有重要意義萨蚕，因此利用單細胞數(shù)據(jù)進行細胞通訊分析是單細胞高級分析的一大重點。

CellChat

CellChat是一款2021年發(fā)表于Nature Communications的單細胞細胞通訊分析工具蹄胰。

CellChat上游分析

1. 安裝 CellChat

devtools::install_github("sqjin/CellChat")

2. 細胞通訊分析

這里我們使用 pbmc3k 的公共數(shù)據(jù)集為例來一起學習如何利用 CellChat 進行細胞通訊分析岳遥。

2.1 數(shù)據(jù)輸入

CellChat 需要的輸入文件包括：

• 細胞的基因表達矩陣（已經(jīng)經(jīng)過normalize的）
  不同的基因作為行名，細胞ID作為列名裕寨。這一點和Seurat對象的結(jié)構(gòu)是保持一致的浩蓉。
• 細胞的metadata信息
  細胞層面的信息，這里可以是我們上游分析的注釋信息宾袜，可以直接把SeuratObject的metadata提取出來使用捻艳。

所以我們的輸入文件可以這樣準備：

library(CellChat)
library(Seurat)
library(SeuratData)

data("pbmc3k.final")

data <- GetAssayData(object = pbmc3k.final, slot = 'data')
meta <- pbmc3k.final@meta.data

這里涉及到一個小知識點：SeuratObject中RNA assay中不同slot所存儲的是什么值：

• counts：原始的基因counts數(shù)，也就是簡單reads計數(shù)的結(jié)果庆猫，對于10X Genomics平臺數(shù)據(jù)來說认轨，這是cellranger運行的結(jié)果；
• data：這是原始表達矩陣經(jīng)過質(zhì)控之后進行NormalizeData()的數(shù)據(jù)月培，NormalizeData去除了不同細胞之間測序深度的差異好渠，同時對結(jié)果進行了對數(shù)化，這個數(shù)據(jù)是CellChat想要的數(shù)據(jù)节视。
• scale.data：這是函數(shù)ScaleData()運行的結(jié)果拳锚，主要是將每個基因的表達量轉(zhuǎn)換成了符合標準正態(tài)分布的數(shù)據(jù)，從而降低部分細胞異常表達值的影響寻行。

因為實際上我們在使用CellChat時是都已經(jīng)默認上游的處理已經(jīng)完成了霍掺，所以我在這里不打算介紹CellChat本身自帶的函數(shù)去normalize原始counts的分析。

2.2 創(chuàng)建CellChat對象

CellChat對象和Seurat對象很像，我們可以這樣來進行創(chuàng)建：

cellchat <- createCellChat(object = data,
                           meta = meta,
                           group.by = 'seurat_annotations')

這里的group.by參數(shù)是不可少的杆烁，來自于我們的metadata的列名牙丽，我們一般指定為細胞類型注釋的結(jié)果，這是有利于我們的分析結(jié)果的兔魂。

除此之外烤芦，不得不介紹的還有這個函數(shù)不僅僅支持將你自己提取的表達矩陣作為輸入，還支持直接用SeuratObject作為輸入析校，不過需要注意的是构罗，如果是多組樣本整合的SeuratObject，我們不能利用integrated assay作為輸入智玻，因為其含有負值遂唧，所以我們可用下面的這段代碼實現(xiàn)和前面代碼一樣的目的：

cellchat <- createCellChat(object = pbmc3k.final,
                           meta = meta,
                           group.by = 'seurat_annotations',
                           assay = 'RNA')

此外還有一個參數(shù)do.sparse不要去改動它（默認為TRUE），用稀疏矩陣處理單細胞數(shù)據(jù)能夠節(jié)省更多的空間和時間吊奢。

簡單介紹一下CellChat數(shù)據(jù)的結(jié)構(gòu)：

str(cellchat)

Formal class 'CellChat' [package "CellChat"] with 14 slots
  ..@ data.raw      : num[0 , 0 ] 
  ..@ data          :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. ..@ i       : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. ..@ p       : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. ..@ Dim     : int [1:2] 13714 2638
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC" "AAACATTGAGCTAC" "AAACATTGATCAGC" "AAACCGTGCTTCCG" ...
  .. .. ..@ x       : num [1:2238732] 1.64 1.64 2.23 1.64 1.64 ...
  .. .. ..@ factors : list()
  ..@ data.signaling: num[0 , 0 ] 
  ..@ data.scale    : num[0 , 0 ] 
  ..@ data.project  : num[0 , 0 ] 
  ..@ net           : list()
  ..@ netP          : list()
  ..@ meta          :'data.frame':  2638 obs. of  7 variables:
  .. ..$ orig.ident        : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ nCount_RNA        : num [1:2638] 2419 4903 3147 2639 980 ...
  .. ..$ nFeature_RNA      : int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ...
  .. ..$ seurat_annotations: Factor w/ 9 levels "Naive CD4 T",..: 2 4 2 3 7 2 5 5 1 6 ...
  .. ..$ percent.mt        : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ...
  .. ..$ RNA_snn_res.0.5   : Factor w/ 9 levels "0","1","2","3",..: 2 4 2 3 7 2 5 5 1 6 ...
  .. ..$ seurat_clusters   : Factor w/ 9 levels "0","1","2","3",..: 2 4 2 3 7 2 5 5 1 6 ...
  ..@ idents        : Factor w/ 9 levels "Naive CD4 T",..: 2 4 2 3 7 2 5 5 1 6 ...
  ..@ DB            : list()
  ..@ LR            : list()
  ..@ var.features  : list()
  ..@ dr            : list()
  ..@ options       :List of 1
  .. ..$ mode: chr "single"

截止到目前為止盖彭，我們的CellChat對象結(jié)構(gòu)如上所示，可以看到：

• 我們輸入的數(shù)據(jù)存在了data這個稀疏矩陣里面页滚，而后面的data.signaling等對象還有待后面的分析進行填充召边；
• 此外，我們應該重點關(guān)注的就是meta部分裹驰，這個部分是一個數(shù)據(jù)框掌实，所以我們還是可以在創(chuàng)建CellChat對象之后來更改細胞的信息，這是比較方便的邦马，可以看到目前CellChat對象對細胞的分析是基于seurat_annotations的贱鼻，這是我們前面設置好的，如果我們想要改變只需要setIdent(cellchat, ident.use = "labels")就可以了滋将，同時我們還可以通過levels(cellchat@idents)來查看當前的細胞分組信息邻悬。

2.3 選擇受體配體數(shù)據(jù)庫

CellChat 有一個專門的數(shù)據(jù)庫，叫做CellChatDB随闽，這個數(shù)據(jù)庫是 CellChat 的作者們通過閱讀大量文獻父丰，手動整理出來的“受體-配體”對，目前有人掘宪、鼠以及斑馬魚的版本蛾扇。其中人的叫做 CellChatDB.human，鼠的叫做 CellChatDB.mouse魏滚，斑馬魚的叫做 CellChatDB.zebrafish镀首。關(guān)于這兩個數(shù)據(jù)庫中具體的“受體-配體”對信息可以通過showDatabaseCategory()函數(shù)獲得，在這里就不贅述了鼠次。

• mouse
  2,021 validated molecular interactions, including 60% of secrete autocrine/paracrine signaling interactions, 21% of extracellular matrix (ECM)-receptor interactions and 19% of cell-cell contact interactions.
• human
  1,939 validated molecular interactions, including 61.8% of paracrine/autocrine signaling interactions, 21.7% of extracellular matrix (ECM)-receptor interactions and 16.5% of cell-cell contact interactions.

加載數(shù)據(jù)庫：

CellChatDB <- CellChatDB.human

來看一下數(shù)據(jù)庫的結(jié)構(gòu)（option）：

head(CellChatDB$interaction)

最后更哄，千萬不要忘了把數(shù)據(jù)庫嵌入到 CellChat 對象中芋齿，否則后面的分析會報錯：

cellchat@DB <- CellChatDB

有的時候我們并不想分析所有的細胞通訊，例如我們只關(guān)心Secreted Signaling成翩，或者我們只相信經(jīng)過KEGG分析驗證過的細胞通訊觅捆，為了節(jié)約運行空間和時間，我們可以使用subsetDB()函數(shù)來取數(shù)據(jù)庫的子集麻敌，簡單查看一下數(shù)據(jù)庫的數(shù)據(jù)結(jié)構(gòu)：

dplyr::glimpse(CellChatDB$interaction)

Rows: 1,939
Columns: 11
$ interaction_name   <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB3_TGFBR1_TGFBR2", "TGFB1_ACVR1B_TGFB…
$ pathway_name       <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb",…
$ ligand             <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB2", "TGFB2", "TGFB3", "TGFB3", "TGFB1…
$ receptor           <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFbR2", "ACVR1C_TGFbR2", "ACVR1B_TGFbR2…
$ agonist            <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb a…
$ antagonist         <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antag…
$ co_A_receptor      <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",…
$ co_I_receptor      <chr> "TGFb inhibition receptor", "TGFb inhibition receptor", "TGFb inhibition receptor", "TG…
$ evidence           <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa04350", "PMID: 27449815", "PMID: 27449815…
$ annotation         <chr> "Secreted Signaling", "Secreted Signaling", "Secreted Signaling", "Secreted Signaling",…
$ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2)", "TGFB3 - (TGFBR1+TGFBR2)", "TGFB1…

我們只選擇Secreted Signaling相關(guān)的細胞間通訊信息：

subsetDB(CellChatDB.human, search = 'Secreted Signaling')

當然我們能做的絕不僅這些栅炒，我們可以通過更改subsetDB()函數(shù)的key參數(shù)來實現(xiàn)任何標準的篩選，例如我們想篩選出evidence為KEGG的細胞間通訊：

subsetDB(CellChatDB.human, search = '^KEGG', key = 'evidence') #regular expression here!!!

2.4 細胞通訊鑒定

在這個部分我們?yōu)榱私档瓦\行的內(nèi)存和時間消耗术羔，我們會抽取部分細胞進行分析赢赊，分析的流程也很好理解：首先，鑒定出細胞中高表達的受體和配體編碼基因聂示，然后再依據(jù)這個表達譜構(gòu)建細胞之間的受體與配體作用對。和單細胞分析一樣簇秒，這里我們也可以使用 future 包來進行并行計算鱼喉，加速我們的分析過程。

library(future)
plan(strategy = 'multiprocess', workers = 4)
#subset
cellchat <- subsetData(cellchat)

為什么這里要做subset趋观，因為為了節(jié)省運算時間和空間扛禽，通過這一步，我們后面的分析將只關(guān)注于與細胞通訊有關(guān)的基因（從數(shù)據(jù)框中提取出來的）皱坛，所以在這里針對于表達矩陣的基因取了子集编曼。

cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

上面兩個函數(shù)分別鑒定除了每個細胞群中高表達的細胞通訊相關(guān)基因和細胞群之間根據(jù)細胞通訊相關(guān)基因的表達情況（identifyOverExpressedGenes()）最終鑒定出的細胞間通訊關(guān)系（identifyOverExpressedInteractions）。

cellchat <- computeCommunProb(cellchat)
cellchat <- filterCommunication(cellchat, min.cells = 10)

注意到為了獲得更可信的細胞間通訊剩辟，方便后續(xù)的驗證掐场，這里首先對每個通訊計算了相應的概率值，同時進行了基于每個細胞類群中支持該通訊的細胞數(shù)量的過濾贩猎。

需要注意的是熊户，預測出的細胞間通訊的數(shù)量是和每個細胞群的基因表達量有關(guān)的，所以對于基因平均表達量的計算就顯得很重要吭服，在默認情況下嚷堡，computeCommunProb()函數(shù)使用的是trimean方法，該方法默認如果一群細胞里面不足25%的細胞表達某個基因的話艇棕，這個基因在該群細胞里面的平均表達量就是0蝌戒。不過你可以通過設置type = "truncatedMean"和trim=來自己指定這個閾值，比如trim=0.1就表示閾值為10%沼琉。顯然默認方法trimean能幫我們篩選出更少的北苟、更可信的細胞間通訊。此外打瘪，考慮到細胞數(shù)更多的細胞群之間的細胞通訊信號往往會更強粹淋，為了去除不同細胞群之間的細胞數(shù)量差異吸祟，我們可以嘗試使用population.size = TRUE來輔助我們發(fā)現(xiàn)稀有細胞類群之間的細胞通訊。你可以通過computeAveExpr(cellchat, features = c("CXCL12","CXCR4"), type = "truncatedMean", trim = 0.1)來提取你所感興趣的細胞通訊相關(guān)基因的平均表達量信息桃移。

2.5 細胞通訊與信號通路關(guān)系的構(gòu)建

細胞之間的通訊往往會是信號通路的重要組成部分屋匕，把細胞通訊放到信號通路中進行理解可能會更有利于我們理解生物學過程。

cellchat <- computeCommunProbPathway(cellchat)

每對受配體的細胞間通訊網(wǎng)絡會被存儲到net的slot中借杰，而每個信號通過的細胞間通訊網(wǎng)絡信息將會被存儲到netP網(wǎng)絡中过吻。

2.6 細胞通訊網(wǎng)絡的構(gòu)建

進一步的，我們將細胞間的通訊進行整合蔗衡，就能構(gòu)建出細胞間的通訊網(wǎng)絡纤虽。

cellchat <- aggregateNet(cellchat)

當然，你可能只關(guān)心部分細胞之間的通訊绞惦，你可以通過指定 信號發(fā)出細胞 和 信號作用細胞 來進行個性化的分析：

?aggregateNet
#output
aggregateNet(
  object,
  sources.use = NULL,
  targets.use = NULL,
  signaling = NULL,
  pairLR.use = NULL,
  remove.isolate = TRUE,
  thresh = 0.05,
  return.object = TRUE
)

也就是這里的 sources.use 和 targets.use逼纸。那么指定成什么呢？這個時候metadata信息的作用就來了济蝉，指定的就是前面 group.by 所包含的細胞類群信息杰刽。

至此，上游分析已經(jīng)算是結(jié)束王滤，下面就是如何對信息進行展示贺嫂，即 可視化。