计算两个单细胞dataset的相关性

日常掰瞎

  对于bulk-RNA来说，计算两个样本的相关性可以选取共有基因的表达值直接来计算相关系数，从而得到两个样本的关系性。那么单细胞如何计算呢？首先，我们需要明白一点，在单细胞层面我们可以将每个细胞看作是一个样本，这样原始样本就可以看成是细胞样本的集合。这样一来，计算两个原始样本的相关性就变成计算对应细胞样本的相关性了，而细胞样本的相关性计算与bulk-RNA一致。方法明白了，下面咱就来看看具体实现代码吧！

代码展示

  这里以10x的数据格式来演示，比如有两个样本A、B，每个样本结果文件有barcodes.tsv、genes.tsv、matrix.mtx三个文件。

library(Seurat)
library(Matrix)
library(ggplot2)

read_count_output <- function(dir) {
  dir <- normalizePath(dir, mustWork = TRUE)
  m <- readMM(paste0(dir, "/", "matrix.mtx"))
  genes <- read.table(paste0(dir, "/", 'genes.tsv'), stringsAsFactors = F,sep='\t',header = F)$V2
  barcodes <- readLines(file(paste0(dir, "/", "barcodes.tsv")))
  colnames(m) <- barcodes
  rownames(m) <- genes
  return(m)
}

countA <- read_count_output('10x_result/sampleA')
countB <- read_count_output('10x_result/sampleB')

over_barcode <- intersect(colnames(countA), colnames(countB))
over_gene <- intersect(rownames(countA), rownames(countB))
sub_countA  <- countA[over_gene, over_barcode]
sub_countB  <- countB[over_gene, over_barcode]

cortest <- function(x,y){
    result <- cor.test(x,y)
    coefficient <- result$estimate
    pvalue <- result$p.value
    return(data.frame(coefficient=coefficient, pvalue=pvalue))
}

countAlist <- as.list(data.frame(countA))
countBlist <- as.list(data.frame(countB))
outcor <- data.frame(t(mapply(cortest, countAlist , countAlist , SIMPLIFY=T)))

ggplot(data.frame(corrlation='corrlation',value=unlist(outcor[,1])),aes(x=corrlation,y=value))+geom_boxplot()+xlab('')+
theme(axis.title = element_text(size=16),axis.text = element_text(size=16))

最后输出相关性箱线图如下：

从箱线图上可以看出，对应细胞样本之间的相关性基本在0.95以上，从而得到样本A、B的相关性很高。

结束语

  这里我们才用的策略是计算原始样本中对应细胞的相关性的方法来估计原始样本的相关性。当然你也可以用其他的方法来估计样本的相关性，比如把scRNA数据变成bulk-RNA数据再算相关性，我觉得这种方法实现起来更简单一些，感兴趣的可以试一试！