當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

哇！单细胞测序-配体受体互作分析原来可以这么简单又高大上！

發布時間：2025/3/15 编程问答 31 豆豆

生活随笔收集整理的這篇文章主要介紹了哇！单细胞测序-配体受体互作分析原来可以这么简单又高大上！小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

NGS系列文章包括NGS基礎、轉錄組分析（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析（ChIP-seq基本分析流程）、單細胞測序分析 (重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程（原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。

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

于是，我一不小心發現了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 metadataser.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）)

# 為細胞類型增加metadatacell.types <- vector("logical",length=ncol(ser.obj)) names(cell.types) <- colnames(ser.obj)cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="0"] <- "cd4.tconv" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="1"] <- "cd4.tconv" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="2"] <- "b.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="3"] <- "b.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="4"] <- "cd14.monocytes" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="5"] <- "cd8.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="6"] <- "cd4.tconv" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="7"] <- "cd4.tconv" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="8"] <- "b.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="9"] <- "b.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="10"] <- "nk.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="11"] <- "cd8.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="12"] <- "plasma.cells" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="13"] <- "cd14.monocytes" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="14"] <- "cd16.monocytes" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="15"] <- "pdc" cell.types[ser.obj@meta.data$RNA_snn_res.0.5=="16"] <- "RBCs"ser.obj[["cell.types"]] <- cell.types#去除紅細胞rbc.cell.names <- names(cell.types)[ser.obj@meta.data$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 celltalkerexpr.mat <- GetAssayData(ser.obj,slot="counts") defined.clusters <- ser.obj@meta.data$cell.types defined.groups <- ser.obj@meta.data$sample.type defined.replicates <- ser.obj@meta.data$sample.idreshaped.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組特有的受體-配體互作

CIRCOS圈圖繪制 - circos安裝
CIRCOS圈圖繪制 - 最簡單繪圖和解釋
CIRCOS圈圖繪制 - 染色體信息展示和調整
CIRCOS增加熱圖、點圖、線圖和區塊屬性

作者：May（協和醫院）

編輯：生信寶典

總結

以上是生活随笔為你收集整理的哇！单细胞测序-配体受体互作分析原来可以这么简单又高大上！的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇： UCI机器学习数据集
下一篇：最后一周|高级转录组分析和R语言数据可视