celltalker配体受体互作
2020-04-25 本文已影响0人
生信编程日常
对单细胞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)
image.png
欢迎关注~
公众号二维码.jpg
参考：
https://arc85.github.io/celltalker/articles/celltalker.html#clustering-data-with-seurat
