相似關(guān)鍵詞

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

適用背景

近年來，單細胞轉(zhuǎn)錄組測序異?；鸨c其類似的空間轉(zhuǎn)錄組測序很有可能成為下一個風口。但是兩者還是存在比較明顯的區(qū)別的，單細胞數(shù)據(jù)的研究分辨率在全細胞或細胞核，而空間組數(shù)據(jù)則可以從多個細胞到亞細胞的水平，只不過如今空間組技術(shù)還存在較多缺陷，例如mRNA捕獲率低等問題，兩者各有優(yōu)劣，具有一定的可比性，畢竟它們本質(zhì)上的研究對象都是轉(zhuǎn)錄組，甚至還可以與bulk RNA進行比較。

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

當然，本文寫的函數(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ù)這兩個矩陣自行畫圖，當然我也寫了下面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可能會更好。