哇!單細胞測序-配體受體互作分析原來可以這麼簡單又高大上!

單細胞轉錄組測序發展至今，我們發現許多文章的最後一部分都會落到配受體結合，可是如何挑配受體，哪些基因可能是配受體，我一臉懵逼。。。

於是，我一不小心發現了celltalker(https://arc85.github.io/celltalker/articles/celltalker.html#vignette-overview)，大家可以嘗試一下哦，嘻嘻，廢話不多說。

Introduction

對單細胞RNAseq數據可能進行的多種分析之一是評估細胞間的交流 (cell-cell communication)。celltalker通過尋找細胞群內和細胞群之間已知的配體和受體對的表達來評估細胞之間的交流。 我們採用的配受體資料庫來自Ramilowski等人於2015年在Nature communication上發表的A draft network of ligand-receptor-mediated multicellular signalling in human描述的一組配體和受體。我們建議使用此數據集作為起點，並整理自己的已知配體和受體列表。另外Tormo 2018年發表的Nature文章Single-cell reconstruction of the early maternal-fetal interface in humans擴展了受體和配體對，也會應用於cellTalker的更新版中。

為了獲得可靠的結果，我們要求每個組中都有多個重複樣品，並且只對不同組間一致性表達的配體和受體感興趣（而僅在單個重複中發現的互作可信度低）。我們通過評估每組中各個重複的表達矩陣並僅對滿足一定閾值（這個閾值隨意性也比較強）的相互作用進行提取。

配體和受體相互作用的差異至少在三種方面具有生物學意義：

在一組細胞中獨特地存在；

各個cluster間配體或受體的互作差異；

參與組間配體和受體相互作用的細胞網絡的差異。

我們提供了評估每種潛在生物學差異的方法，並提供了具體示例。

在這個vignette中，我們展示了cellTalker在評估健康捐獻者外周血（N = 2）和扁桃體（N = 3）中鑑定配體/受體相互作用的基本應用。該數據可從我們最近發布的數據集GSE139324中獲得 (Cillo et al, Immunity 2020)。

Vignette overview

展示Celltalker應用於10X Genomics數據的的標準用法。具體分為下面幾步：

使用標準的Seurat工作流程（v.3.1.1）對數據進行聚類;

使用Celltalker建立樣品組中穩定表達的配體和受體的列表;

確定配體/受體相互作用;

評估組之間特異表達的配體/受體對;

識別和可視化組特異的配體/受體對;

Installation

library(devtools)
install_github("arc85/celltalker")
library(celltalker)

Clustering data with Seurat

使用Seurat進行標準的聚類分析和免疫譜系識別（假設已從GEO下載了raw matrix）。（重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述））

suppressMessages({
library(Seurat)
library(celltalker)
})

# 設置可重複性的種子數字

set.seed(02221989)

# 讀取raw data

# path.to.working：自行修改路徑
path.to.working <- 「」
setwd(paste(path.to.working,"/data_matrices/",sep=""))
data.paths <- list.files()
# GRCh38 根據需要調整為其它基因組版本
specific.paths <- paste(path.to.working,"data_matrices",data.paths,"GRCh38",sep="/")
setwd(path.to.working)

raw.data <- Read10X(specific.paths)

# 設置metadata （這一部分是這個數據特異的，實際還是自己整理個metadata文件更為方便）

sample.data <- data.frame(matrix(data=NA,nrow=ncol(raw.data),ncol=2))
rownames(sample.data) <- colnames(raw.data)
colnames(sample.data) <- c("sample.id","sample.type")

sample.data[grep("^[A-z]",rownames(sample.data)),"sample.id"] <- "pbmc_1"
sample.data[grep("^2",rownames(sample.data)),"sample.id"] <- "tonsil_1"
sample.data[grep("^3",rownames(sample.data)),"sample.id"] <- "pbmc_2"
sample.data[grep("^4",rownames(sample.data)),"sample.id"] <- "tonsil_2"
sample.data[grep("^5",rownames(sample.data)),"sample.id"] <- "pbmc_3"
sample.data[grep("^6",rownames(sample.data)),"sample.id"] <- "tonsil_3"

sample.data[,"sample.type"] <- sapply(strsplit(sample.data$sample.id,split="_"),function(x) x[1])

## 使用原始數據創建Seurat對象，並添加sample-specific metadata

ser.obj <- CreateSeuratObject(counts=raw.data,meta.data=sample.data)

Seurat三部曲

#標準Seurat工作流程
ser.obj <- NormalizeData(ser.obj)
ser.obj <- FindVariableFeatures(ser.obj)
ser.obj <- ScaleData(ser.obj)

PCA分析

ser.obj <- RunPCA(ser.obj)

獲得對各個主成分貢獻比較大的基因 (用了這麼多年的PCA可視化竟然是錯的！！！)

## PC_ 1
## Positive: S100A6, IL32, S100A4, ANXA1, VIM, FTL, TRBC1, SRGN, S100A9, S100A8
## TYROBP, LYZ, CTSW, XIST, NEAT1, VCAN, S100A12, FCER1G, S100A11, FCN1
## PLAC8, ID2, CCL5, NKG7, CST3, CSTA, ZFP36, IL1B, MT2A, KLRB1
## Negative: RGS13, KIAA0101, NUSAP1, AURKB, MKI67, BIRC5, TYMS, TOP2A, TK1, CDKN3
## UBE2C, PTTG1, CDK1, STMN1, CCNB2, GTSE1, BIK, RRM2, TCL1A, SHCBP1
## CDCA3, CDC20, TPX2, LRMP, CCNA2, MND1, CCNB1, PBK, ZWINT, RMI2
## PC_ 2
## Positive: CST3, LYZ, FCN1, CSTA, S100A9, S100A8, TYROBP, LST1, FGL2, VCAN
## S100A12, SERPINA1, MNDA, FCER1G, CLEC7A, MS4A6A, CD14, CFD, IL1B, TYMP
## LGALS1, RP11-1143G9.4, AIF1, CTSS, NAMPT, CFP, TNFSF13B, CSF3R, MPEG1, TMEM176B
## Negative: IL32, NPM1, CD69, TRBC1, ISG20, ITM2A, IGKC, IGHA1, HSP90AB1, DDIT4
## HIST1H4C, PSIP1, AQP3, MYC, PIM2, HMGN1, PASK, NUCB2, HSPA1B, HSPB1
## CD79A, SUSD3, KLRB1, SYNE2, CHI3L2, IGHG3, IGLC2, FKBP11, IGHG1, SH2D1A
## PC_ 3
## Positive: IL32, NKG7, CTSW, TRBC1, GZMA, CST7, GNLY, MKI67, ANXA1, TOP2A
## CCL5, PRF1, BIRC5, S100A4, KLRB1, CCNA2, AURKB, CENPF, GTSE1, CDKN3
## KLRD1, UBE2C, CDK1, TYMS, TPX2, RRM2, ID2, S100A6, FGFBP2, CDC20
## Negative: HLA-DRA, HLA-DQA1, HLA-DQB1, CD79A, HLA-DRB1, MS4A1, CD74, HLA-DPA1, HLA-DPB1, CD79B
## HLA-DMA, HLA-DMB, BANK1, VPREB3, IGKC, HLA-DRB5, MEF2C, CD22, IRF8, CD19
## SMIM14, FCRLA, HLA-DOB, CD24, CD40, FCER2, BLK, HLA-DQA2, IGHD, CTSH
## PC_ 4
## Positive: TOP2A, UBE2C, MKI67, GTSE1, CENPF, AURKB, PLK1, CCNA2, CDK1, CDCA8
## HMMR, CDCA3, CDC20, TPX2, CDKN3, DLGAP5, CENPE, BIRC5, CCNB2, CENPA
## KIF2C, CKAP2L, PBK, NUSAP1, KIFC1, AURKA, SPC25, NUF2, KIF23, ASPM
## Negative: NKG7, GNLY, CST7, GZMB, GZMA, PRF1, KLRD1, FGFBP2, CCL5, KLRF1
## HOPX, CTSW, GZMH, TRDC, FCGR3A, SPON2, CLIC3, MATK, ADGRG1, S1PR5
## CCL4, CMC1, XCL2, PFN1, CD160, FCRL6, IL2RB, TRGC1, KLRC1, C12orf75
## PC_ 5
## Positive: ICA1, PDCD1, TBC1D4, ITM2A, ICOS, MAF, TOX2, IL32, TNFRSF4, PASK
## PKM, SMCO4, ACTG1, CORO1B, CTLA4, NPM1, TRBC1, PCAT29, TIGIT, AC133644.2
## TOX, ANP32B, ENO1, GBP2, COTL1, GAPDH, SUSD3, PIM2, AQP3, SERPINA9
## Negative: NKG7, GNLY, KLRD1, FGFBP2, GZMB, GZMA, KLRF1, CCL5, PRF1, TRDC
## GZMH, CST7, CTSW, BANK1, MATK, PLK1, HMMR, HLA-DPB1, CENPA, CLIC3
## GTSE1, CENPE, CCL4, SPON2, PDLIM1, HLA-DPA1, CDCA8, DLGAP5, TPX2, IGHD

拐點法尋找top可用的主成分用於後續分析 (具體選擇方式見：（重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述））)

ElbowPlot(ser.obj)

降維可視化

＃我們選擇top 15 PCs用於後續分析
ser.obj <- RunUMAP(ser.obj,reduction="pca", dims=1:15)

聚類

## Computing nearest neighbor graph
ser.obj <- FindNeighbors(ser.obj,reduction="pca",dims=1:15)
## Computing SNN
ser.obj <- FindClusters(ser.obj,resolution=0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15524
## Number of edges: 543084
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9185
## Number of communities: 17
## Elapsed time: 2 seconds

畫圖看一看！ggplot2高效實用指南 (可視化腳本、工具、套路、配色)

p1 <- DimPlot(ser.obj,reduction="umap",group.by="sample.id")
p2 <- DimPlot(ser.obj,reduction="umap",group.by="sample.type")
p3 <- DimPlot(ser.obj,reduction="umap",group.by="RNA_snn_res.0.5",label=T) + NoLegend()
cowplot::plot_grid(p1,p2,p3)

讓我們看看部分基因的表達情況！

FeaturePlot(ser.obj,reduction="umap",features=c("CD3D","CD8A","CD4","CD14","MS4A1","FCGR3A","IL3RA"))

命名細胞簇並移除紅細胞 (cellassign：用於腫瘤微環境分析的單細胞注釋工具（9月Nature）)

# 為細胞類型增加metadata

cell.types <- vector("logical",length=ncol(ser.obj))
names(cell.types) <- colnames(ser.obj)

cell.types[[e-mail]$RNA_snn_res.0.5=="0"] <- "cd4.tconv"
cell.types[[e-mail]$RNA_snn_res.0.5=="1"] <- "cd4.tconv"
cell.types[[e-mail]$RNA_snn_res.0.5=="2"] <- "b.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="3"] <- "b.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="4"] <- "cd14.monocytes"
cell.types[[e-mail]$RNA_snn_res.0.5=="5"] <- "cd8.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="6"] <- "cd4.tconv"
cell.types[[e-mail]$RNA_snn_res.0.5=="7"] <- "cd4.tconv"
cell.types[[e-mail]$RNA_snn_res.0.5=="8"] <- "b.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="9"] <- "b.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="10"] <- "nk.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="11"] <- "cd8.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="12"] <- "plasma.cells"
cell.types[[e-mail]$RNA_snn_res.0.5=="13"] <- "cd14.monocytes"
cell.types[[e-mail]$RNA_snn_res.0.5=="14"] <- "cd16.monocytes"
cell.types[[e-mail]$RNA_snn_res.0.5=="15"] <- "pdc"
cell.types[[e-mail]$RNA_snn_res.0.5=="16"] <- "RBCs"

ser.obj[["cell.types"]] <- cell.types

#去除紅細胞

rbc.cell.names <- names(cell.types)[[e-mail]$RNA_snn_res.0.5=="16"]

ser.obj <- ser.obj[,!colnames(ser.obj) %in% rbc.cell.names]

Consistently expressed ligands and receptors

現在，我們已經在數據中識別並命名了cluster，我們將繼續進行celltalker分析。隨該軟體包一起提供的有一個ramilowski_pairs，它是一個由配體、受體和推測的配體受體對組成的data.frame。

首先，根據通用型配體和受體從我們的數據集中選出對應的基因，然後進行差異基因分析，只保留在樣品組之間差異表達的配體受體。

然後，我們將為每個重複樣本創建單獨的數據矩陣，存儲為tibble格式，以便於使用tidyverse進行後續處理。

(生信寶典註：如果我們自己有受體配體對，也可以整理成這樣一個三列的格式，導入進來，替換掉數據包原有的配體受體對信息，實現更加定製的分析。)

#查閱 ramilowski_pairs 數據集

head(ramilowski_pairs)

## ligand receptor pair
## 1 A2M LRP1 A2M_LRP1
## 2 AANAT MTNR1A AANAT_MTNR1A
## 3 AANAT MTNR1B AANAT_MTNR1B
## 4 ACE AGTR2 ACE_AGTR2
## 5 ACE BDKRB2 ACE_BDKRB2
## 6 ADAM10 AXL ADAM10_AXL

dim(ramilowski_pairs)

## [1] 2557 3

#該數據集中有2,557個特異的配體/受體對

#鑑定差異表達的配體和受體

#在我們的數據集中識別配體和受體

ligs <- as.character(unique(ramilowski_pairs$ligand))
recs <- as.character(unique(ramilowski_pairs$receptor))

ligs.present <- rownames(ser.obj)[rownames(ser.obj) %in% ligs]
recs.present <- rownames(ser.obj)[rownames(ser.obj) %in% recs]

genes.to.use <- union(ligs.present,recs.present) # union用於合併子集

# 使用FindAllMarkers區分組之間差異表達的配體和受體

Idents(ser.obj) <- "sample.type"
markers <- FindAllMarkers(ser.obj, assay="RNA",features=genes.to.use,only.pos=TRUE)

ligs.recs.use <- unique(markers$gene)
length(ligs.recs.use)

## [1] 61

#產生61個配體和受體以進行評估

#過濾ramilowski配受對

interactions.forward1 <- ramilowski_pairs[as.character(ramilowski_pairs$ligand) %in% ligs.recs.use,]
interactions.forward2 <- ramilowski_pairs[as.character(ramilowski_pairs$receptor) %in% ligs.recs.use,]
interact.for <- rbind(interactions.forward1, interactions.forward2)
dim(interact.for)

## [1] 241 3

生成celltalker的輸入數據

#產生241個配體和受體相互作用

#Create data for celltalker

expr.mat <- GetAssayData(ser.obj,slot="counts")
defined.clusters <- [e-mail]$cell.types
defined.groups <- [e-mail]$sample.type
defined.replicates <- [e-mail]$sample.id

reshaped.matrices <- reshape_matrices(count.matrix=expr.mat,clusters=defined.clusters,groups=defined.groups,replicates=defined.replicates,ligands.and.receptors=interact.for)

#查看tibble的層次結構
reshaped.matrices

## # A tibble: 2 x 2
## group samples
## <chr> <list>
## 1 pbmc <tibble [3 × 2]>
## 2 tonsil <tibble [3 × 2]>

數據展開為單個樣品展示

unnest(reshaped.matrices,cols="samples")

## # A tibble: 6 x 3
## group sample expr.matrices
## <chr> <chr> <list>
## 1 pbmc pbmc_1 <named list [8]>
## 2 pbmc pbmc_2 <named list [8]>
## 3 pbmc pbmc_3 <named list [8]>
## 4 tonsil tonsil_1 <named list [8]>
## 5 tonsil tonsil_2 <named list [8]>
## 6 tonsil tonsil_3 <named list [8]>

names(pull(unnest(reshaped.matrices,cols="samples"))[[1]])

## [1] "b.cells" "cd14.monocytes" "cd16.monocytes" "cd4.tconv"
## [5] "cd8.cells" "nk.cells" "pdc" "plasma.cells"

在此初始步驟中，我們要做的是將我們的整體表達矩陣中給每個樣本中分離出單獨的表達矩陣。

接下來，使用create_lig_rec_tib函數為每個組創建一組一致表達的配體和受體。

# cells.reqd=10：每個cluster中至少有10個細胞表達了配體/受體
# freq.pos.reqd=0.5：至少有50％重複個體中表達的配體/受體
consistent.lig.recs <- create_lig_rec_tib(exp.tib=reshaped.matrices,clusters=defined.clusters,groups=defined.groups,replicates=defined.replicates,cells.reqd=10,freq.pos.reqd=0.5,ligands.and.receptors=interact.for)

consistent.lig.recs

## # A tibble: 2 x 2
## group lig.rec.exp
## <chr> <list>
## 1 pbmc <tibble [8 × 2]>
## 2 tonsil <tibble [8 × 2]>

unnest(consistent.lig.recs[1,2],cols="lig.rec.exp")

## # A tibble: 8 x 2
## cluster.id ligands.and.receptors
## <chr> <named list>
## 1 b.cells <named list [2]>
## 2 cd14.monocytes <named list [2]>
## 3 cd16.monocytes <named list [2]>
## 4 cd4.tconv <named list [2]>
## 5 cd8.cells <named list [2]>
## 6 nk.cells <named list [2]>
## 7 pdc <named list [2]>
## 8 plasma.cells <named list [2]>

pull(unnest(consistent.lig.recs[1,2],cols="lig.rec.exp")[1,2])[[1]]

根據上面指定的標準，我們現在已經從每個實驗組的每個cluster中獲得了一致表達的配體和受體的列表。

## $ligands
## [1] "S100A9" "S100A8" "IL1B" "FN1" "BTLA" "SPON2"
## [7] "LRPAP1" "VCAN" "CD14" "LY86" "HLA-G" "HLA-A"
## [13] "HLA-E" "HLA-B" "LTB" "HSPA1A" "CD24" "NAMPT"
## [19] "TIMP1" "CD40LG" "ADAM28" "PNOC" "IL7" "ANXA1"
## [25] "SEMA4D" "VIM" "PSAP" "LYZ" "SELPLG" "HMGB1"
## [31] "TNFSF13B" "GZMB" "CALM1" "SERPINA1" "HSP90AA1" "B2M"
## [37] "PKM" "IL16" "CCL5" "CCL3" "ICAM2" "CD70"
## [43] "ICAM1" "ICAM3" "TGFB1" "FLT3LG" "APP"
##
## $receptors
## [1] "CSF3R" "TGFBR3" "KCNA3" "CD53" "CD1A" "CD247"
## [7] "SELL" "CXCR4" "ITGA4" "GRM7" "TGFBR2" "CCR5"
## [13] "TFRC" "TLR1" "IL7R" "CD180" "ADRB2" "CD74"
## [19] "HMMR" "HLA-F" "KCNQ5" "IGF2R" "CCR6" "CD36"
## [25] "CXCR3" "PGRMC1" "CD72" "TGFBR1" "ABCA1" "IFITM1"
## [31] "CD81" "KCNQ1" "CD44" "CD82" "IL10RA" "CD3D"
## [37] "CD3G" "CXCR5" "SORL1" "APLP2" "ITGB1" "FAS"
## [43] "CD27" "CD4" "KLRG1" "CLEC2B" "KLRD1" "KLRC1"
## [49] "PDE1B" "CD63" "TNFRSF17" "CD19" "ITGAL" "ITGAM"
## [55] "TNFRSF13B" "ERBB2" "PTPRA" "CD40" "OPRL1" "INSR"
## [61] "TYROBP" "CD79A" "KCNN4" "FPR2" "LILRB2" "LILRB1"
## [67] "KIR2DL3" "KIR2DL1" "KIR3DL2" "TNFRSF13C" "CELSR1" "ITGB2"

Determine putative ligand/receptor pairs

獲得穩定表達的受體-配體後，鑑定給定樣品組細胞簇內和細胞簇間可能存在的互作 (基於ramilowski_pairs$pair數據)。獲得的Tibble數據 (put.int)中包含樣本組和每個組的配體/受體對列表，以及參與配體/受體相互作用的cluster。

# freq.group.in.cluster: 只對包含細胞數大於總細胞數5%的簇進行互作分析
put.int <- putative_interactions(ligand.receptor.tibble=consistent.lig.recs,clusters=defined.clusters,groups=defined.groups,freq.group.in.cluster=0.05,ligands.and.receptors=interact.for)

## Warning: `cols` is now required.
## Please use `cols = c(lig.rec.exp)`

## Warning: `cols` is now required.
## Please use `cols = c(lig.rec.exp)`

Identifying and visualizing unique ligand/receptor pairs in a group

現在我們有了配體/受體相互作用的列表，我們可以使用unique_interactions函數鑑定組特異的互作，並使用circos_plot函數可視化組之間的差異。

#Identify unique ligand/receptor interactions present in each sample
unique.ints <- unique_interactions(put.int,group1="pbmc",group2="tonsil",interact.for)

#Get data to plot circos for PBMC
pbmc.to.plot <- pull(unique.ints[1,2])[[1]]
for.circos.pbmc <- pull(put.int[1,2])[[1]][pbmc.to.plot]

circos_plot(interactions=for.circos.pbmc,clusters=defined.clusters)

PBMC組特有的受體-配體互作

從以上圖中我們可以看出研究人員用黃色表示配體基因，綠色表示受體基因，然後將可以相互配對的基因連在一起構成簇之間的互作關係。

Tonsil組特有的受體-配體互作

#Get data to plot circos for tonsil
tonsil.to.plot <- pull(unique.ints[2,2])[[1]]
for.circos.tonsil <- pull(put.int[2,2])[[1]][tonsil.to.plot]

circos_plot(interactions=for.circos.tonsil,clusters=defined.clusters)

Tonsil組特有的受體-配體互作

作者：May（協和醫院）

編輯：生信寶典

Reference

https://arc85.github.io/celltalker/articles/celltalker.html#vignette-overview

你可能還想看

微信公眾號近日更改推送機制，推文不再按照時間線顯示，將生信寶典設置星標，多點在看或留言，這樣就不會錯過精彩教程和送書信息啦～

往期精品(點擊關鍵詞直達對應教程)

後臺回復「生信寶典福利第一波」或點擊閱讀原文獲取教程合集