黄思源

V1

2022/07/10阅读:60主题:雁栖湖

答读者问(6):单细胞TPM矩阵如何分析?

之前在平台发布了几个问题,接下来我会依次给出我的想法,仅供参考。本文有些观点读者朋友可能是第一次看到,有不同意见的话,欢迎给我发邮件讨论:huangsiyuan1001@163.com

1问题

  • 一、有的文章只提供TPM的单细胞表达矩阵,可以用seurat分析吗?
  • 二、分析流程和用count矩阵有什么不同?
  • 三、10X的单细胞转录组数据的标准化需要考虑基因长度吗?

2先来看看第3个小问题

10X的单细胞转录组数据的标准化需要考虑基因长度吗?

答案是不需要。

我们看一下seurat里面NormalizeData()函数是如何做标准化,然后求Log。

test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)

LogNormalize: Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor(默认是10000). This is then natural-log transformed using log1p.(每一列都是先除以一个细胞UMI的总和,再乘以10000,再加1,再求对数)

的确没有考虑基因长度 这要从测序原理去解释。以10X的3'端测序为例,磁珠上有很多短序列(它们的cellular barcode是一样的,UMI各不相同),在PolyT结合mRNA的3'端PolyA之后,尽管转录本有长有短,但因为测序序列长度有所限制,这些转录本分子会在靠近3'端的地方随机打断(只留下最靠边的序列),

然后两边连上接头进行后续的测序步骤...

所以表达定量的多少和基因长度关系不大。另外,多说一句,UMI的作用是消除PCR的影响,只要是来源于一个转录本,不管扩增多少次,最后定量值只会加1。

那传统二代转录组测序为什么考虑基因长度呢? 还是从测序原理去解释。一个完整的转录本分子会被随机打断,转录本越长,片段会越多,而这些片段最终都会被测序。因为这个偏好,标准化的时候才需要考虑转录本长度的影响。


3回到第1个小问题

TPM的单细胞表达矩阵,可以用seurat分析吗?(可以)

先来看一下TPM的公式,顺带把RPKM的也看了: 上述公式把转录本长度这一项去掉,就会发现和前面讲的seurat的LogNormalize方法中的标准化这一步是如此类似

"每一列都是先除以一个细胞UMI的总和,再乘以10000"

只是乘上的倍数不同,形式是一样的。把umi/sum(umi) * 10000当作TPM、RPKM的特例。既然如此,TPM、RPKM这一类标准化形式的数据应该也是能走seurat标准流程的,只要把文库大小、转录本长度这些因素考虑了就行。代码略有不同,待会儿再说。


类似的问题,比如从10X单细胞数据中导出类似TPM的矩阵,是不是也很简单了呢?从seurat的count矩阵和data矩阵都能推算出来。 为啥会问这个问题?之前有读者问我运行cibersort的时候,需要导入单细胞得到的参考矩阵,需要类似TPM的数据(要和bulk数据保持一致),行是基因,列是细胞类型。

在得到类TPM矩阵后,手动计算每种细胞类型的均值,当然也可以用seurat的AverageExpression这个函数,结果一样。

### 10X导出类tpm矩阵 #####
### 第一种方法
test.seu=readRDS("D:/hsy/bioinformatics/PTJ/011_paper/CR/malignant.seu.1129.rds")
mat.tpm.like=as.data.frame(test.seu[["RNA"]]@counts) %>% apply(2,function(x){x/sum(x) * 10000})
#减少计算量,只选择一个亚群
selectCB=test.seu@meta.data[test.seu@meta.data$sample == "P14","CB"]
mat.tpm.like=mat.tpm.like[,selectCB]
#手动计算均值
selectCB.tpm.mean=as.data.frame(rowMeans(mat.tpm.like))
colnames(selectCB.tpm.mean)="method1_P14"

### 第二种方法
Idents(test.seu)="sample"
averageexp=AverageExpression(test.seu)
averageexp=averageexp$RNA
averageexp=as.data.frame(averageexp[,"P14"])
colnames(averageexp)="method2_P14"

### 比较两种方法
twomethod=cbind(averageexp,selectCB.tpm.mean)
twomethod$d=twomethod$method2_P14-twomethod$method1_P14
twomethod$d=round(twomethod$d,5)

> table(twomethod$d)
    0 
33538
# 差值都是0,说明两种方法求出来是一样的

相信你现在对AverageExpression函数有了一个更深的理解,现在知道它返回的是什么了吧。

我想这也是为什么早期有几篇文章描述10X的数据会用到TPM这个词:

Our data were transformed into log 2 (TPM+1), where transcripts per million (TPM) was calculated as the proportion of UMIs of a gene within each cell multiplied by 1,000,000. (来自文献: Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma)

OK, 上述都是我个人的理解。seurat官方GitHub是这么说的:

  • If you have TPM data, you can simply manually log transform the gene expression matrix in the object@data slot before scaling the data.
  • You have to replace your object@data slot with the desired gene expression matrix as follows: pbmc@data = log(x = norm + 1))

(分别来自https://github.com/satijalab/seurat/issues/668和https://github.com/satijalab/seurat/issues/747)


4第2个小问题

在分析流程上,这种矩阵和用count矩阵有什么不同?

只有前几步有点区别,

### 常规流程
test.seu=CreateSeuratObject(counts = testdf)
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
......

### TPM数据
test.seu1=CreateSeuratObject(counts = log(mat_tpm+1))
test.seu1 <- FindVariableFeatures(test.seu1, selection.method = "vst", nfeatures = 2000)
# 或者
test.seu=CreateSeuratObject(counts = mat_tpm)
test.seu[["RNA"]]@data=as(as.matrix(log(mat_tpm + 1)), "dgCMatrix")
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
......

后面的流程一样。有人会担心,提供TPM,不能找差异基因。事实上,seurat默认找差异基因的方法并没有用count,而是data矩阵,所以还是可以用FindMarkers()的。


本文图片来自谷歌和10X官网;本文有些观点各位读者朋友大概率是第一次看到,仅代表我个人观点。

欢迎点赞、转发、转载

分类:

其他

标签:

其他

作者介绍

黄思源
V1