celltalker配体受体互作

生信编程日常

发布于 2020-05-04 21:42:09

1.7K0

发布于 2020-05-04 21:42:09

文章被收录于专栏：生物信息学、python、R、linux生物信息学、python、R、linux

对单细胞RNAseq数据常规的分析经常会有对细胞间的相互作用的分析(cell-cell communication)。这里介绍一个简单的分析包celltalker，通过寻找细胞群内和细胞群之间已知的配体和受体对的表达来评估细胞之间的交流。

主要流程

使用标准的Seurat工作流程（v.3.1.1）对数据进行聚类使用Celltalker建立一个组一致表达的配体和受体的列表确定与一组的推定配体/受体相互作用评估不同组间的配体/受体作用识别和可视化一组中独特的配体/受体对

安装

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

使用seurat进行常规分析

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

#Set seed for reproducibility

set.seed(02221989)

#Read in raw data

setwd(paste(path.to.working,"/data_matrices/",sep=""))
data.paths <- list.files()
specific.paths <- paste(path.to.working,"data_matrices",data.paths,"GRCh38",sep="/")
setwd(path.to.working)

raw.data <- Read10X(specific.paths)

#Create 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])

#Create a Seurat object with associated metadata

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

#Standard Seurat workflow

ser.obj <- NormalizeData(ser.obj)
ser.obj <- FindVariableFeatures(ser.obj)
ser.obj <- ScaleData(ser.obj)

ser.obj <- RunPCA(ser.obj)
ElbowPlot(ser.obj)

#We will select the first 15 PCs to use
ser.obj <- RunUMAP(ser.obj,reduction="pca",dims=1:15)

ser.obj <- FindNeighbors(ser.obj,reduction="pca",dims=1:15)
ser.obj <- FindClusters(ser.obj,resolution=0.5)

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)

image.png

#定义细胞类型
#Add metadata for cell types

cell.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

#Let's remove RBCs from our analysis

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]

celltalker分析

完成分群并定义细胞类型后，将继续进行celltalker分析。随该软件包一起提供的是data.frame“ ramilowski_pairs”，它是一个由配体，受体和推定的proxy_receptor对组成的data.frame。如果对如何构造此data.frame感兴趣，请参考“ data-raw”文件夹和“ create_ramolowski_pairs_data.R”文件。

首先，我们将进行差异基因表达分析，以将潜在配体和受体的范围缩小到两个样品组之间差异表达的范围。我们确定总体数据集中存在的配体和受体（来自ramilowski_pairs），然后使用这些配体和受体进行差异表达。

我们将为每个重复样本创建单独的数据矩阵，存储为tibble格式，以便于使用tidyverse进行后续处理。

#Check out ramilowski_pairs data.frame
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

#There are 2,557 unique ligand/receptor interactions in this dataset
#Identification of differentially expressed ligands and receptors
#Identify ligands and receptors in our dataset

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)

#Use FindAllMarkers for differentially expressed ligands and receptors between groups

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)

#Yields 61 ligands and receptors to evaluate
#Filter ramilowski pairs
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)

#Yields 241 ligand and receptor interactions to evaluate
#Create data for celltalker

expr.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.id

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

#Check out the hierarchy of the tibble
reshaped.matrices

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

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

每个组创建一组一致表达的配体和受体

#接下来，我们可以使用create_lig_rec_tib函数为每个组创建一组一致表达的配体和受体。
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

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

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

确定推定的配体/受体对

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)

识别和可视化一组中独特的配体/受体对

#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)

image.png

#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)

参考：

https://arc85.github.io/celltalker/articles/celltalker.html#clustering-data-with-seurat

http://blog.sciencenet.cn/blog-118204-1220235.html

本文参与腾讯云自媒体分享计划，分享自作者个人站点/博客。

如有侵权请联系 cloudcommunity@tencent.com 删除

html

本文分享自作者个人站点/博客前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

html

登录后参与评论

0 条评论

热度