空間組數(shù)據(jù)和單細胞數(shù)據(jù)的相關(guān)性分析（Seurat）2022-05-20

相似關(guān)鍵詞

  • 單細胞數(shù)據(jù)集相關(guān)性分析
  • 空間組與單細胞數(shù)據(jù)集相關(guān)性分析
  • 空間組數(shù)據(jù)集相關(guān)性分析

適用背景

近年來宾毒，單細胞轉(zhuǎn)錄組測序異惩叮火爆侣诺，而與其類似的空間轉(zhuǎn)錄組測序很有可能成為下一個風(fēng)口拳锚。但是兩者還是存在比較明顯的區(qū)別的益楼，單細胞數(shù)據(jù)的研究分辨率在全細胞或細胞核猾漫，而空間組數(shù)據(jù)則可以從多個細胞到亞細胞的水平，只不過如今空間組技術(shù)還存在較多缺陷感凤，例如mRNA捕獲率低等問題悯周，兩者各有優(yōu)劣，具有一定的可比性陪竿，畢竟它們本質(zhì)上的研究對象都是轉(zhuǎn)錄組禽翼，甚至還可以與bulk RNA進行比較。

類似于單細胞，空間組的注釋也成為具有爭議性的難題闰挡，比較空間組中的單元很大可能不是一個細胞锐墙，有可能是多個細胞或者是幾分之一個細胞，但肯定是可以看成是與某種細胞類型相似的單元长酗。因此構(gòu)建單細胞細胞類型等分組信息與空間組單元的相關(guān)性溪北，可以對空間組的單元進行類似細胞類型的注釋，這對于空間組之后的研究具有重要意義夺脾。

當(dāng)然刻盐，本文寫的函數(shù)可以對任意兩個或多個空間組或單細胞數(shù)據(jù)集進行相關(guān)性分析，而不局限于只對比空間組與單細胞的數(shù)據(jù)劳翰。

之所以會寫這個函數(shù)敦锌，是因為最近接到任務(wù)要畫下面這種圖，于是去看了文章的源代碼佳簸，提取里面的相關(guān)腳本改寫了一下乙墙，如果感興趣也可以訪問其GitHub主頁研究研究。

相關(guān)性圖片示例

主函數(shù)

加載需要的R包

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

本文主要包含3個函數(shù)生均，其中一個函數(shù)cortable用于生成相關(guān)性矩陣和顯著性檢驗听想，而函數(shù)corplot和cor_heatmap是可視化函數(shù).

cortable

這個函數(shù)看起來有點長，主要分為2個部分：相關(guān)性矩陣計算和顯著性檢驗马胧。因為只比較2組汉买，所以數(shù)據(jù)傳入也分為2組，以列表的形式傳入佩脊，一組數(shù)據(jù)里又可以包含多個小組蛙粘，一個小組需要傳入4個參數(shù)，這4個參數(shù)構(gòu)成一個小組代表一個數(shù)據(jù)集：

  • Seurat對象威彰，例如自帶的pbmc_small出牧；
  • Seurat對象的assays名稱，一般單細胞用“RNA”歇盼，而空間組用“Spatial”;
  • 分組列名舔痕，例如細胞類型所在列；
  • 標簽豹缀，對此數(shù)據(jù)集進行標記伯复。
    此外，另外3個參數(shù)是可選項：
  • features可以指定進行相關(guān)性分析的基因列表邢笙，不傳人則代表使用所有基因
  • overlap是指是否要對數(shù)據(jù)集小組內(nèi)部也要計算相關(guān)性系數(shù)啸如，默認不計算
  • corr.method指定相關(guān)性系數(shù)的類型，默認使用spearman相關(guān)性系數(shù)
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 相關(guān)矩陣計算
###分組計算平均表達量
###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個大組的數(shù)據(jù)
###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)
###計算相關(guān)性系數(shù)
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))
}

上面的函數(shù)返回的是一個列表鸣剪，包含2個矩陣组底，相關(guān)性矩陣和對應(yīng)的顯著性矩陣丈积，可以自行根據(jù)這兩個矩陣自行畫圖，當(dāng)然我也寫了下面2個可視化函數(shù)债鸡。

corplot

此函數(shù)是基于corrplot包寫的江滨，是為了復(fù)原文章中的圖，其實感覺這個包很不友好厌均，所以后面寫了另一個函數(shù)唬滑。

  • cor.table 相關(guān)性矩陣
  • pva.table 顯著性矩陣
  • cutoff 對相關(guān)性系數(shù)進行過濾，默認不過濾
  • pf 文件名棺弊，可選項
  • col 調(diào)色板設(shè)置
  • order 排序方法晶密，其它方法可以查看corrplot官網(wǎng)，注意如果用hclust需要設(shè)置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

因為領(lǐng)導(dǎo)說corrplot畫出來的圖太丑了，所以用pheatmap包重新畫了圖侈净。

  • cor.table 相關(guān)性矩陣
  • pva.table 顯著性矩陣
  • cutoff 對相關(guān)性系數(shù)進行過濾尊勿，默認不過濾
  • pf 文件名，可選項
  • col 調(diào)色板設(shè)置
  • res 圖片dpi值畜侦，調(diào)整清晰度
  • 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這個數(shù)據(jù)集展示

> 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)
結(jié)果展示

overlap="true"的結(jié)果


結(jié)果展示

總結(jié)與補充

cortable函數(shù)可以用features參數(shù)指定基因集元扔，因為我們測試發(fā)現(xiàn)用所有基因進行分析效果并不好，用TF或者DEG可能會更好旋膳。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末澎语，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子验懊，更是在濱河造成了極大的恐慌擅羞，老刑警劉巖，帶你破解...
    沈念sama閱讀 216,372評論 6贊 498
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件鲁森，死亡現(xiàn)場離奇詭異祟滴，居然都是意外死亡，警方通過查閱死者的電腦和手機歌溉，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,368評論 3贊 392
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來骑晶，“玉大人痛垛，你說我怎么就攤上這事⊥盎祝” “怎么了匙头？”我有些...
    開封第一講書人閱讀 162,415評論 0贊 353
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵，是天一觀的道長仔雷。 經(jīng)常有香客問我蹂析，道長舔示，這世上最難降的妖魔是什么？ 我笑而不...
    開封第一講書人閱讀 58,157評論 1贊 292
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任电抚，我火速辦了婚禮惕稻，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘蝙叛。我一直安慰自己俺祠，他們只是感情好，可當(dāng)我...
    茶點故事閱讀 67,171評論 6贊 388
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布借帘。 她就那樣靜靜地躺著蜘渣，像睡著了一般。 火紅的嫁衣襯著肌膚如雪肺然。 梳的紋絲不亂的頭發(fā)上蔫缸，一...
    開封第一講書人閱讀 51,125評論 1贊 297
  • 城市分裂傳說
    那天，我揣著相機與錄音际起，去河邊找鬼捂龄。 笑死，一個胖子當(dāng)著我的面吹牛加叁，可吹牛的內(nèi)容都是我干的倦沧。 我是一名探鬼主播，決...
    沈念sama閱讀 40,028評論 3贊 417
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼它匕，長吁一口氣：“原來是場噩夢啊……” “哼展融！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起豫柬，我...
    開封第一講書人閱讀 38,887評論 0贊 274
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤告希，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后烧给，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體燕偶，經(jīng)...
    沈念sama閱讀 45,310評論 1贊 310
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 37,533評論 2贊 332
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年础嫡，在試婚紗的時候發(fā)現(xiàn)自己被綠了指么。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點故事閱讀 39,690評論 1贊 348
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡榴鼎，死狀恐怖伯诬，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情巫财，我是刑警寧澤盗似，帶...
    沈念sama閱讀 35,411評論 5贊 343
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站平项，受9級特大地震影響赫舒，放射性物質(zhì)發(fā)生泄漏悍及。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點故事閱讀 41,004評論 3贊 325
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一接癌、第九天 我趴在偏房一處隱蔽的房頂上張望心赶。 院中可真熱鬧，春花似錦扔涧、人聲如沸园担。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,659評論 0贊 22
  • 一樁弒父案枯夜，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽弯汰。三九已至，卻和暖如春湖雹，著一層夾襖步出監(jiān)牢的瞬間咏闪，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,812評論 1贊 268
  • 情欲美人皮
    我被黑心中介騙來泰國打工摔吏， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留鸽嫂，地道東北人。 一個月前我還...
    沈念sama閱讀 47,693評論 2贊 368
  • 代替公主和親
    正文 我出身青樓征讲，卻偏偏與公主長得像据某，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子诗箍，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 44,577評論 2贊 353

推薦閱讀更多精彩內(nèi)容