空间组数据和单细胞数据的相关性分析（Seurat）2022-05-20

相似关键词

  • 单细胞数据集相关性分析
  • 空间组与单细胞数据集相关性分析
  • 空间组数据集相关性分析

适用背景

近年来，单细胞转录组测序异常火爆，而与其类似的空间转录组测序很有可能成为下一个风口。但是两者还是存在比较明显的区别的，单细胞数据的研究分辨率在全细胞或细胞核，而空间组数据则可以从多个细胞到亚细胞的水平，只不过如今空间组技术还存在较多缺陷，例如mRNA捕获率低等问题，两者各有优劣，具有一定的可比性，毕竟它们本质上的研究对象都是转录组，甚至还可以与bulk RNA进行比较。

类似于单细胞，空间组的注释也成为具有争议性的难题，比较空间组中的单元很大可能不是一个细胞，有可能是多个细胞或者是几分之一个细胞，但肯定是可以看成是与某种细胞类型相似的单元。因此构建单细胞细胞类型等分组信息与空间组单元的相关性，可以对空间组的单元进行类似细胞类型的注释，这对于空间组之后的研究具有重要意义。

当然，本文写的函数可以对任意两个或多个空间组或单细胞数据集进行相关性分析，而不局限于只对比空间组与单细胞的数据。

之所以会写这个函数，是因为最近接到任务要画下面这种图，于是去看了文章的源代码，提取里面的相关脚本改写了一下，如果感兴趣也可以访问其GitHub主页研究研究。

相关性图片示例

主函数

加载需要的R包

library(ggplot2)
library(Seurat)
library(corrplot)
library(pheatmap)

本文主要包含3个函数，其中一个函数cortable用于生成相关性矩阵和显著性检验，而函数corplot和cor_heatmap是可视化函数.

cortable

这个函数看起来有点长，主要分为2个部分：相关性矩阵计算和显著性检验。因为只比较2组，所以数据传入也分为2组，以列表的形式传入，一组数据里又可以包含多个小组，一个小组需要传入4个参数，这4个参数构成一个小组代表一个数据集：

  • Seurat对象，例如自带的pbmc_small；
  • Seurat对象的assays名称，一般单细胞用“RNA”，而空间组用“Spatial”;
  • 分组列名，例如细胞类型所在列；
  • 标签，对此数据集进行标记。
    此外，另外3个参数是可选项：
  • features可以指定进行相关性分析的基因列表，不传人则代表使用所有基因
  • overlap是指是否要对数据集小组内部也要计算相关性系数，默认不计算
  • corr.method指定相关性系数的类型，默认使用spearman相关性系数
cortable <- function(ref=list(c(pbmc_small,"RNA","groups","label1")),
                    que=list(c(pbmc_small,"RNA","groups","label2")),
                    features=NULL,overlap="false",corr.method = 'spearman') {
##Part 1 相关矩阵计算
###分组计算平均表达量
###FUN1 AVERAGEEXPRESSION 
preave <- function(obj,
                   features=c(rownames(pbmc_small)),
                   group="groups",
                   assay="RNA",
                   label="label1") {
ob1 <- obj
DefaultAssay(ob1) <- assay
grp1 <- group
Idents(ob1)<-ob1[[grp1]]
ave1 <- AverageExpression(ob1,features = features,assays=assay)
ave1 <- ave1[[1]]
colnames(ave1) <- paste0(label,"_",colnames(ave1))

Sp1 = ave1
  avg = rowMeans(Sp1)
  Sp1 = sweep(Sp1,1,avg,"/")
  rm(avg)

Sp1[is.nan(Sp1)] <- 0
ave1 <- Sp1
return(ave1)
}
###获取所有基因
###FUN2 GET FEATURES
getfeatures <- function(ref=list(c(pbmc_small,"RNA","groups","label1")),
                        que=list(c(pbmc_small,"RNA","groups","label2"))) {
lenref <- length(ref)
lenque <- length(que)

reflist <- list()
for (i in 1:lenref) {
tmp <- ref[[i]][[1]]
DefaultAssay(tmp) <- ref[[i]][[2]]
reflist[[i]] <- tmp
}
geneset1 <- lapply(reflist[1:lenref],rownames)
gene1 <- Reduce(intersect, geneset1)

quelist <- list()
for (i in 1:lenque) {
tmp <- que[[i]][[1]]
DefaultAssay(tmp) <- que[[i]][[2]]

quelist[[i]] <- tmp
}
geneset2 <- lapply(quelist[1:lenque],rownames)
gene2 <- Reduce(intersect, geneset2)

final_featues <- intersect(gene1,gene2)
return(final_featues)
}
if (is.null(features)) {
features <- getfeatures(ref,que)
}else{
features1 <- getfeatures(ref,que)
features <- intersect(features,features1)
}
###合并2个大组的数据
###FUN3 MERGE AVERAGEEXPRESSION
preave_list <- function(inlist, features=NULL) {
len <- length(inlist)
matlist <- list()
for (i in 1:length(inlist)) {
matlist1 <- preave(obj=inlist[[i]][[1]],
                  features=features,
                  assay=inlist[[i]][[2]],
                  group=inlist[[i]][[3]],
                  label=inlist[[i]][[4]])
matlist[[i]] <- matlist1
}
mat <- Reduce(cbind,matlist)
return(mat)
#lapply(inlist,preave(),features=features,group=)
}

###Compute cor values
Sp1 = preave_list(ref,features=features)
colnames(Sp1) <- paste0(colnames(Sp1),"_ref")
Sp2 = preave_list(que,features=features)
colnames(Sp2) <- paste0(colnames(Sp2),"_que")

geTable = merge(Sp1,Sp2, by='row.names', all=F)
###计算相关性系数
rownames(geTable) = geTable$Row.names
geTable = geTable[,2:ncol(geTable)]
# corr.method = c('spearman', 'pearson') etc.
corr.method <- corr.method
Corr.Coeff.Table = cor(geTable,method=corr.method)

##Part 2 显著性检验
###Estimate p-value
nPermutations = 1000
shuffled.cor.list = list()
  pb   <- txtProgressBar(1, 100, style=3)

  for (i in 1:nPermutations){
    shuffled = apply(geTable[,1:ncol(Sp1)],1,sample)
    shuffled2 = apply(geTable[,(ncol(Sp1)+1):ncol(geTable)],1,sample)
    shuffled = cbind(t(shuffled),t(shuffled2))
    shuffled.cor = cor(shuffled,method=corr.method)
    shuffled.cor.list[[i]] = shuffled.cor
    rm(list=c('shuffled','shuffled2','shuffled.cor'))
    if ((i %% 100) ==0){
      setTxtProgressBar(pb, (i*100)/nPermutations)
    }
  }

  p.value.table = matrix(ncol=ncol(geTable), nrow = ncol(geTable))

  rownames(p.value.table) = colnames(geTable)
  colnames(p.value.table) = colnames(geTable)

  shuffled.mean.table = matrix(ncol=ncol(geTable), nrow = ncol(geTable))
  rownames(shuffled.mean.table) = colnames(geTable)
  colnames(shuffled.mean.table) = colnames(geTable)

  a = combn(1:ncol(geTable),2)
  for (i in 1:ncol(a)){
    cor.scores = sapply(shuffled.cor.list,"[",a[1,i],a[2,i])
    shuffled.mean.table[a[1,i],a[2,i]] = mean(cor.scores)
    shuffled.mean.table[a[2,i],a[1,i]] = mean(cor.scores)
    p.value = mean(abs(cor.scores)>=abs(Corr.Coeff.Table[a[1,i],a[2,i]]))
    p.value.table[a[1,i],a[2,i]] = p.value
    p.value.table[a[2,i],a[1,i]] = p.value
    rm(list=c('cor.scores','p.value'))
    setTxtProgressBar(pb, (i*100)/ncol(a))
  }
if (overlap=="false") {
M <- p.value.table
mat <- M[,grep('ref',colnames(M))]
mat <- mat[grep('que',rownames(M)),]
p.value.table  <- mat

M <- Corr.Coeff.Table
mat <- M[,grep('ref',colnames(M))]
mat <- mat[grep('que',rownames(M)),]
Corr.Coeff.Table <- mat
}

return(list(Corr.Coeff.Table,p.value.table))
}

上面的函数返回的是一个列表，包含2个矩阵，相关性矩阵和对应的显著性矩阵，可以自行根据这两个矩阵自行画图，当然我也写了下面2个可视化函数。

corplot

此函数是基于corrplot包写的，是为了复原文章中的图，其实感觉这个包很不友好，所以后面写了另一个函数。

  • cor.table 相关性矩阵
  • pva.table 显著性矩阵
  • cutoff 对相关性系数进行过滤，默认不过滤
  • pf 文件名，可选项
  • col 调色板设置
  • order 排序方法，其它方法可以查看corrplot官网，注意如果用hclust需要设置overlap=“true”，不然会报错，这也是其一个bug吧
  • wid 宽值
  • hei 高值
  • label.size 文字大小
corplot <- function(cor.table=NULL,
                    pva.table=NULL,
                    cutoff=0,
                    pf=NULL,
                    col=colorRampPalette(c("darkblue", "white","darkred")),
                    order="original",
                    wid=1000,
                    hei=1000,
                    label.size=1) {

cor.table[cor.table<cutoff] <- 0
if (is.null(pf)) {
pf <- paste0("corrplot_filter",cutoff,"_",order,".jpg")
}
jpeg(pf,wid,hei)
print(
      corrplot(cor.table, order=order,tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(cor.table),max(cor.table)), is.corr=F,tl.cex=label.size, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=pva.table, col=col(200), main= paste(pf),mar=c(3,1,5,1),cl.align.text="l")
      )
dev.off()

}

cor_heatmap

因为领导说corrplot画出来的图太丑了，所以用pheatmap包重新画了图。

  • cor.table 相关性矩阵
  • pva.table 显著性矩阵
  • cutoff 对相关性系数进行过滤，默认不过滤
  • pf 文件名，可选项
  • col 调色板设置
  • res 图片dpi值，调整清晰度
  • wid 宽值
  • hei 高值
  • scale 是否对坐标轴进行缩放
cor_heatmap <- function(cor.table,pva.table=NULL,
                        col= colorRampPalette(c("darkblue", "white","darkred"))(256),
                        pf=NULL,wid=1000,hei=1000,res=1200,
                        cutoff=0,
                        scale="none"){

cor.table[cor.table<cutoff] <- 0
if (is.null(pva.table)) {
p1<-pheatmap(cor.table,scale=scale,show_colnames = T,show_rownames = T,fontsize = 10,
             cluster_rows = T,cluster_cols = T,border_color = "NA",color = col,res=res,filename=NA)
}else{
pmt <- pva.table
ssmt <- pmt< 0.01
pmt[ssmt] <-'**'
smt <- pmt >0.01& pmt <0.05
pmt[smt] <- '*'
pmt[!ssmt&!smt]<- ''

p1<-pheatmap(cor.table,scale=scale,show_colnames = T,show_rownames = T,fontsize = 10,
             cluster_rows = T,cluster_cols = T,border_color = "NA",color = col, res=res,filename=NA,
             display_numbers = pmt,fontsize_number = 12, number_color = "black")

}
if (is.null(pf)) {
pf <- paste0("corheatmap_filter",cutoff,".jpg")
}
jpeg(pf,wid,hei)
print(p1)
dev.off()

}

运行示例

使用pbmc_small这个数据集展示

> head(pbmc_small)
                  orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.8
ATGCCAGAACGACT SeuratProject         70           47               0
CATGGCCTGTGCAT SeuratProject         85           52               0
GAACCTGATGAACC SeuratProject         87           50               1
TGACTGGATTCTCA SeuratProject        127           56               0
AGTCAGACTGCACA SeuratProject        173           53               0
TCTGATACACGTGT SeuratProject         70           48               0
TGGTATCTAAACAG SeuratProject         64           36               0
GCAGCTCTGTTTCT SeuratProject         72           45               0
GATATAACACGCAT SeuratProject         52           36               0
AATGTTGACAGTCA SeuratProject        100           41               0
               letter.idents groups RNA_snn_res.1
ATGCCAGAACGACT             A     g2             0
CATGGCCTGTGCAT             A     g1             0
GAACCTGATGAACC             B     g2             0
TGACTGGATTCTCA             A     g2             0
AGTCAGACTGCACA             A     g2             0
TCTGATACACGTGT             A     g1             0
TGGTATCTAAACAG             A     g1             0
GCAGCTCTGTTTCT             A     g1             0
GATATAACACGCAT             A     g1             0
AATGTTGACAGTCA             A     g1             0

代码示例

ob1 <- pbmc_small
ob2 <- pbmc_small
ob3 <- pbmc_small

ref <- list(c(ob1,"RNA","letter.idents","ob1"),
            c(ob3,"RNA","RNA_snn_res.1","ob3"))
que <- list(c(ob2,"RNA","groups","ob2"),
            c(ob3,"RNA","RNA_snn_res.1","ob3"))

tmp <- cortable(ref,que)
Corr.Coeff.Table <- tmp[[1]]
p.value.table <- tmp[[2]]

comp_table.species1 <- Corr.Coeff.Table
p_table.species1 <- p.value.table

corplot(cor.table=Corr.Coeff.Table,pva.table=p.value.table,cutoff=0,hei=1800)
cor_heatmap(Corr.Coeff.Table,p.value.table)
结果展示

overlap="true"的结果


结果展示

总结与补充

cortable函数可以用features参数指定基因集，因为我们测试发现用所有基因进行分析效果并不好，用TF或者DEG可能会更好。

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 211,265评论 6赞 490
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,078评论 2赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 156,852评论 0赞 347
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,408评论 1赞 283
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,445评论 5赞 384
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,772评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,921评论 3赞 406
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,688评论 0赞 266
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,130评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,467评论 2赞 325
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,617评论 1赞 340
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,276评论 4赞 329
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,882评论 3赞 312
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,740评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,967评论 1赞 265
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,315评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,486评论 2赞 348

推荐阅读更多精彩内容