ixxmu
commented
1 year ago 单细胞高级分析(23): 单细胞联合bulk转录组（上） by TOP生物信息
目前单细胞联合bulk转录组的分析思路有很多，今天介绍其中用得最多的套路，之所以说它是套路，是因为很多数据挖掘为主的生信文章都喜欢用这种分析思路。
大致思路就是：
找到合适的公开的单细胞数据集，按照scRNA-seq常规流程分析一遍，降维聚类注释。需要单细胞数据的目的：单细胞分析能精确得到不同celltype的表达特征，而这层信息正是第二步所参考的； 下载 TCGA某种/些癌种的bulk转录组表达数据，借助细胞类型占比估计软件估计bulk肿瘤样本中各种celltype的占比；将celltype的占比和患者 临床表型（最常见的就是生存时间）联系起来，得出celltype与表型相关等等结论。
而说到细胞类型占比估计软件就不得不提CIBERSORT（最新版叫CIBERSORTx，网址https://cibersortx.stanford.edu/），作者Aaron Newman是计算生物领域的大牛，其团队还开发了单细胞、空间组学相关工具，如CytoTRACE、CytoSPACE
接下来会以CIBERSORTx为例，走完上述的step1 2 3。
1.准备单细胞数据集
2. 基于CIBERSORTx的细胞占比估计
2.1 如何构建signature matrix
2.2 推断TCGA样本的细胞类型占比
3. 联系生存预后
题外话
参考
因为这个专题的篇幅很长，我分了上下两篇。第一篇讲到2.1 如何构建signature matrix，第二篇讲余下部分。
如何获取这个专题所有的代码和示例数据？公众号后台回复8位日期可知：20230602
1.准备单细胞数据集
这次的dataset仍然来自那篇我在B站讲过很多次的鼻咽癌单细胞文章（Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma）。
跟着Cell Research学单细胞分析系列视频（UP主：TOP菌）
早在2021年，我讲单细胞基础绘图的时候，就把这篇文章的数据从原始fastq文件重新跑了一遍。
考虑到这套数据的特点：
成纤维细胞、内皮细胞等基质细胞非常少； 肿瘤细胞由于病人间/内异质性，很难存在共有的表达特征（不同于基质细胞，基质细胞会有 足够的共有的表达谱）。
因此，在演示时删掉非免疫细胞，下文的讨论集中在免疫细胞的占比估计（实际分析中可以灵活处理）。
删掉细胞数很少的6个样本，剩余10个样本。最终的图如下：
2. 基于CIBERSORTx的细胞占比估计
CIBERSORTx主要涉及的是step2：细胞类型占比估计。
我不推荐用R版本的CIBERSORT（是一个叫CIBERSORT.R的代码文件，GitHub能搜到很多包都调用了这个文件），原因是代码非常旧，我在一些包中看到的这个CIBERSORT.R文件甚至是2016年的。
# CIBERSORT R script v1.04 (last updated 10-24-2016)
另一方面，CIBERSORTx的网站足够好用。
如果你无法登录CIBERSORTx官网（https://cibersortx.stanford.edu/），可以试试借助魔法。以下链接就是我用的魔法：
https://tagss01.pro#/register?invite=0LZOXWz72.1 如何构建signature matrix
下图是step1+step2的流程图，最左边的step1准备单细胞数据集我们已经做完了，剩下的都是step2。我们可以看到构建细胞类型参考矩阵这一环节有两种做法：自己找marker构建；利用cibersortx网页构建。在下文我会比较两种构建signature matrix的做法，并给出我所推荐的做法（过程有点长，不关心这个问题的话，你可以直接跳到2.1节的结尾、2.2节的开头看结论）。
为了比较两种做法的好坏，我们开启
上帝视角，做一个已知正确答案的模拟测试：我们将step1中涉及到的10个样本的单细胞数据做成pesudobulk样本。具体做法是把每个样本中所有细胞在各个基因上的umi加和，得到20000基因10列的新矩阵。借助两种构建参考矩阵方法得到的两个signature matrix（一次用一个），把新的伪bulk样本当作常规的bulk去做细胞占比估计。由于真实的细胞占比是知道的，因此就能比较两种方法的优劣。
2.1.1 CIBERSORTx构建
我们需要从单细胞数据集对应的Seurat对象导出表达谱矩阵（single cell reference sample file，后续会上传CIBERSORTx网页），该文件的格式需要满足一些条件：
数值不能取对数。 以制表符分隔的txt文件，没有双引号和缺失值NA。 第1列基因；第1行是细胞类型label。相同类型的细胞应具有相同的label；删除未定义的细胞。 CIBERSORTx会自动对输入数据进行标准化，使得每个样本内部的数值之和相同。如果提供基因长度矫正的表达矩阵（例如RPKM），则signature matrix将以TPM（每百万转录本）表示。如果提供count矩阵，则signature matrix将以CPM（每百万计数）表示。 无论输入什么，请确保signature matrix和后续bulk转录组矩阵都尽可能使用相同的标准化方式。
第4点可能不好理解，稍微解释一下：
第4点提到的RPKM,TPM,CPM是几种常见的bulk转录组标准化方式。而且
RPKM/sum(RPKM)*10^6可以转换为TPM；
count/sum(count)*10^6可以转换为CPM
因为在这次演示中，后续的实际TCGA样本得到的是TPM矩阵，所以此处我们仍然希望得到TPM形式的signature matrix。
按照第4点，我们需要RPKM形式的单细胞表达矩阵。那10X单细胞数据能不能得到类似的矩阵？答案是可以，参考我之前写的答读者问(6)：单细胞TPM矩阵如何分析？
说白了就是，10X单细胞数据不需要对长度矫正，所以CP10K，CPM, TPM, RPKM这几个标准化方法
意义是相同的，只是倍数不同。
library(Seurat)
library(tidyverse)
allseu = readRDS("10个样本的免疫细胞.rds")
### single cell reference sample file
mat.RPKM.like=as.data.frame(allseu[["RNA"]]@counts[,sample(1:26931,5000)]) %>%
apply(2,function(x){x/sum(x) * (10^9)})
mat.RPKM.like=mat.RPKM.like[rowSums(mat.RPKM.like) > 0,]
allanno = allseu@meta.data[,c("CB","newtype")]
smallanno = allanno[colnames(mat.RPKM.like),]
colnames(mat.RPKM.like) = smallanno$newtype
write.table(mat.RPKM.like,file = "mat.RPKM.like.txt",quote = F,sep = "\t",row.names = T,col.names = T)
这个文本文件的左上角没有"Gene"，也就是第一列的列名，用电脑文本编辑器加上即可。最终的单细胞表达谱矩阵是这样的：
随后导入CIBERSORTx生成signature matrix
上传数据
之后会进入新的界面（https://cibersortx.stanford.edu/upload.php）
网页生成signature matrix
进入CIBERSORTx的分析界面（https://cibersortx.stanford.edu/runcibersortx.php）
等待几分钟，结果出来后，下载即可
下载后的signature matrix长这样
2.1.2 自己构建
如果是自己根据marker来构建，就省事很多。
library(Seurat)
library(tidyverse)
allseu = readRDS("10个样本的免疫细胞.rds")
#更改默认的细胞identity
Idents(allseu) = "newtype"
ref_marker=FindAllMarkers(allseu,logfc.threshold = 0.8,min.pct = 0.1,only.pos = T,test.use = "wilcox")
ref_marker=ref_marker%>%filter(p_val_adj < 0.01)
ref_marker$d=ref_marker$pct.1 - ref_marker$pct.2
#ref_marker %>% ggplot(aes(x=d,fill = cluster))+geom_density()+facet_grid(cluster~.)+geom_vline(xintercept = 0.1)
#借助这张图选择一个合适的阈值
ref_marker=ref_marker%>%filter(d > 0.1)
ref_marker=ref_marker%>%arrange(cluster,desc(avg_log2FC))
ref_marker=as.data.frame(ref_marker)
used_gene=sort(unique(ref_marker$gene))
sig_matrix=AverageExpression(allseu, slot = 'data' ,verbose = FALSE,features = used_gene)
sig_matrix=sig_matrix$RNA * 100
#这里要乘以100是因为seurat的scale.factor是10000，而TPM/CPM的定义是1百万，为了保持一致这里再乘以100
sig_matrix=as.data.frame(sig_matrix)
write.table(sig_matrix,file = "ref_sig_by_DEG.txt",quote = F,sep = "\t",row.names = T,col.names = T)
结果文件补上第一列的列名之后就是这样的：
2.1.3 哪种方法好
我们先构建pesudobulk样本，同时统计每个样本真实的细胞类型占比
library(Seurat)
library(tidyverse)
allseu = readRDS("10个样本的免疫细胞.rds")
### 真实占比
ground_truth = prop.table(table(allseu@meta.data$newtype,allseu@meta.data$sample),margin = 2) %>% as.data.frame()
colnames(ground_truth) = c("newtype","sample","ratio")
### 构建pesudobulk样本
bulkmat = data.frame()
for (si in unique(allseu@meta.data$sample)) {
onesample = allseu@meta.data[allseu@meta.data$sample == si,"CB"]
onemat = allseu[["RNA"]]@counts[,onesample]
onemat = rowSums(onemat) %>% as.data.frame()
colnames(onemat) = si
if(which(unique(allseu@meta.data$sample) == si) == 1){
bulkmat = onemat
} else {
bulkmat = cbind(bulkmat,onemat)
}
}
bulkmat = bulkmat[rowSums(bulkmat) > 0,]
随后将count矩阵标准化为TPM
由于该pesudobulk数据的特殊性（每一个count数值都不受gene length的影响），在计算TPM的时候，公式里不加gene length。
###转换为TPM
bulkmat.tpm = bulkmat %>% apply(2,function(x){x/sum(x) * (10^6)})
write.table(bulkmat.tpm,file = "bulkmat.tpm.txt",quote = F,sep = "\t",row.names = T,col.names = T)
补上基因列名，就是这样的
接着我们用CIBERSORTx来推断细胞类型占比（真实占比我们已经知道）
上传数据
上传页面（https://cibersortx.stanford.edu/upload.php）红色是文件类型，要分别选择选择：
Impute Cell Fractions
分析页面（https://cibersortx.stanford.edu/runcibersortx.php）
做两次：一次用方法1生成的signature matrix，一次用方法2的signature matrix。这里仅演示第一次运行，第二次运行同理。
结果出来后，点击下载即可，用于DIY画图
准确性评估
将真实结果和推断结果画成散点图就能看出来了，主要看截距、斜率、相关性等指标。从结果来看，自己找marker构建signature matrix得到的结果略好一些（当然前提是marker找得好）。
https://mp.weixin.qq.com/s/4W5kav6r4UcX-haQ2rv3NA