基于Seurat结果推断单细胞群肿瘤纯度之ESTIMATE / 四六文摘

Inferring tumour purity and stromal and immune cell admixture from expression data.，NC,3013

单细胞转录组是揭示细胞异质性的的有力武器，鉴于肿瘤的异质性，这一点在肿瘤样本中表现尤为突出。所以肿瘤样本的单细胞转录组就不只是无监督地分个群那么简单，基于我们对肿瘤样本已经积累起来的生物学背景（如TCGA），我们可以从更多侧面来反映和说明肿瘤样本的异质性。

今天给大家介绍一款根据stromal和immune细胞比例估算肿瘤纯度的工具：ESTIMATE。之前是基于bulk表达谱来做的，在简书已经有详细的介绍了：

文章解读：
利用表达数据计算基质打分与免疫打分进而预测肿瘤纯度 --- ESTIMATE

代码实践：
使用ESTIMATE来根据stromal和immune细胞比例估算肿瘤纯度

ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) is a tool for predicting tumor purity, and the presence of infiltrating stromal/immune cells in tumor tissues using gene expression data. ESTIMATE algorithm is based on single sample Gene Set Enrichment Analysis and generates three scores:

stromal score (that captures the presence of stroma in tumor tissue),
immune score (that represents the infiltration of immune cells in tumor tissue), and
estimate score (that infers tumor purity).

可以看到和目前单细胞转录组中有些基于富集的细胞类型定义还是很像的，根据一个基因list通过某种规则（这里是ssGSEA）来对细胞打分，进而推断出细胞的类型。

所以正如生信技能树在简书上所言：

其实对大部分使用该包的的文章来说，需要的反而是该包定义的2个基因集，stromal 和 immune , 列表是：

1StromalSignature estimate DCN PAPPA SFRP4 THBS2 LY86 CXCL14 FOXF1 COL10A1 ACTG2 APBB1IP SH2D1A SULF1 MSR1 C3AR1 FAP PTGIS ITGBL1 BGN CXCL12 ECM2 FCGR2A MS4A4A WISP1 COL1A2 MS4A6A EDNRA VCAM1 GPR124 SCUBE2 AIF1 HEPH LUM PTGER3 RUNX1T1 CDH5 PIK3R5 RAMP3 LDB2 COX7A1 EDIL3 DDR2 FCGR2B LPPR4 COL15A1 AOC3 ITIH3 FMO1 PRKG1 PLXDC1 VSIG4 COL6A3 SGCD COL3A1 F13A1 OLFML1IGSF6 COMP HGF GIMAP5 ABCA6 ITGAM MAF ITM2A CLEC7A ASPN LRRC15 ERG CD86 TRAT1 COL8A2 TCF21 CD93 CD163 GREM1 LMOD1TLR2 ZEB2 C1QB KCNJ8 KDR CD33 RASGRP3 TNFSF4 CCR1 CSF1R BTK MFAP5 MXRA5 ISLR ARHGAP28 ZFPM2 TLR7 ADAM12 OLFML2B ENPP2 CILP SIGLEC1 SPON2 PLXNC1 ADAMTS5 SAMSN1 CH25H COL14A1 EMCN RGS4 PCDH12 RARRES2 CD248 PDGFRB C1QA COL5A3 IGF1 SP140TFEC TNN ATP8B4 ZNF423 FRZB SERPING1 ENPEP CD14 DIO2 FPR1 IL18R1 HDC TXNDC3 PDE2A RSAD2 ITIH5 FASLG MMP3 NOX4 WNT2 LRRC32 CXCL9 ODZ4 FBLN2 EGFL6 IL1B SPON1 CD200 2 3ImmuneSignature estimate LCP2 LSP1 FYB PLEK HCK IL10RA LILRB1 NCKAP1L LAIR1 NCF2 CYBB PTPRC IL7R LAPTM5 CD53 EVI2BSLA ITGB2 GIMAP4 MYO1F HCLS1 MNDA IL2RG CD48 AOAH CCL5 LTB GMFG GIMAP6 GZMK LST1 GPR65 LILRB2 WIPF1 CD37 BIN2 FCER1G IKZF1 TYROBP FGL2 FLI1 IRF8 ARHGAP15 SH2B3 TNFRSF1B DOCK2 CD2 ARHGEF6 CORO1A LY96 LYZ ITGAL TNFAIP3 RNASE6TGFB1 PSTPIP1 CST7 RGS1 FGR SELL MICAL1 TRAF3IP3 ITGA4 MAFB ARHGDIB IL4R RHOH HLA-DPA1 NKG7 NCF4 LPXN ITK SELPLG HLA-DPB1 CD3D CD300A IL2RB ADCY7 PTGER4 SRGN CD247 CCR7 MSN ALOX5AP PTGER2 RAC2 GBP2 VAV1 CLEC2B P2RY14 NFKBIAS100A9 IFI30 MFSD1 RASSF2 TPP1 RHOG CLEC4A GZMB PVRIG S100A8 CASP1 BCL2A1 HLA-E KLRB1 GNLY RAB27A IL18RAP TPST2 EMP3 GMIP LCK IL32 PTPRCAP LGALS9 CCDC69 SAMHD1 TAP1 GBP1 CTSS GZMH ADAM8 GLRX PRF1 CD69 HLA-B HLA-DMA CD74 KLRK1 PTPRE HLA-DRA VNN2 TCIRG1 RABGAP1L CSTA ZAP70 HLA-F HLA-G CD52 CD302 CD27

其实这R包的函数本不多：

1# 安装一下 2library(utils) 3rforge <- "http://r-forge.r-project.org" 4install.packages("estimate", repos=rforge, dependencies=TRUE)

在这里我们就不再?estimateScore来一步一步执行示例数据了，虽然对于新手来说这一步往往是不能省略的。实在想看就看技能树的吧：使用ESTIMATE来根据stromal和immune细胞比例估算肿瘤纯度

我就想，这么好的工具单细胞能不能使用呢？我也是表达谱啊，应该是可以的吧，如果可以，Seurat的数据格式可不可以直接做呢？

带着一系列疑问我们来试试。

无疑，作为表达谱我是合格的。关键就是数据格式的问题了。我们发现这个R包操作的都是路径，就连表达谱要求的都是：?estimate::filterCommonGenes。如果想传入一个Seurat的对象我们是要改造一个函数了。

1myfilterCommonGenes <- edit(estimate::filterCommonGenes)

看到这个函数的文件读入是用read.table的：

1function (input.f, output.f, id = c("GeneSymbol", "EntrezID")) 2{ 3 stopifnot((is.character(input.f) && length(input.f) == 1 && 4 nzchar(input.f)) || (inherits(input.f, "connection") && 5 isOpen(input.f, "r"))) 6 stopifnot((is.character(output.f) && length(output.f) == 7 1 && nzchar(output.f))) 8 id <- match.arg(id) 9 input.df <- read.table(input.f, header = TRUE, row.names = 1, 10 sep = "\t", quote = "", stringsAsFactors = FALSE) 11 merged.df <- merge(common_genes, input.df, by.x = id, by.y = "row.names") 12 rownames(merged.df) <- merged.df$GeneSymbol 13 merged.df <- merged.df[, -1:-ncol(common_genes)] 14 print(sprintf("Merged dataset includes %d genes (%d mismatched).", 15 nrow(merged.df), nrow(common_genes) - nrow(merged.df))) 16 outputGCT(merged.df, output.f) 17}

我想直接把seurat的某个gene-cell矩阵对象给它，于是就把输入改成：

1myfilterCommonGenes <- function(input.f, output.f, id = c("GeneSymbol", "EntrezID")) 2{ 3 4 id <- match.arg(id) 5 input.df <- input.f 6 merged.df <- merge(common_genes, input.df, by.x = id, by.y = "row.names") 7 rownames(merged.df) <- merged.df$GeneSymbol 8 merged.df <- merged.df[, -1:-ncol(common_genes)] 9 print(sprintf("Merged dataset includes %d genes (%d mismatched).", 10 nrow(merged.df), nrow(common_genes) - nrow(merged.df))) 11 outputGCT(merged.df, output.f) 12}

为了让改造后的函数依然在estimated的环境之中：

1environment(myfilterCommonGenes) <- environment(estimate::filterCommonGenes)

我们来试试，读入我们熟悉的pbmc数据，注意这里的数据仅作为Seurat的演示示例：

1pbmc <- readRDS('G:\\Desktop\\Desktop\\RStudio\\single_cell\\filtered_gene_bc_matrices\\hg19pbmc_tutorial.rds') 2 3DefaultAssay(pbmc) <- "RNA" 4myfilterCommonGenes(input.f= as.matrix(pbmc@assays$RNA@data), output.f=paste0("matrixgenes.gct"), id="GeneSymbol")

一行提示：

1[1] "Merged dataset includes 7295 genes (3117 mismatched)."

同时，在当前路径下生成了matrixgenes.gct,然后我们启动打分程序：

1?estimateScore # 还是建议大家多问 2# c("affymetrix", "agilent", "illumina") 3estimateScore(paste0("matrixgenes.gct"),paste0("estimate_score.gct"), platform= "affymetrix") 4 5[1] "1 gene set: StromalSignature overlap= 58" 6[1] "2 gene set: ImmuneSignature overlap= 141" # 可以看到overlap的基因比较少，我们毕竟是10X的数据啊。

我们比较感兴趣的就是这个打分文件estimate_score.gct,我们来看那看是怎样的，以及如何把它写入Seurat对象中。

1estimate_score <- read.table(paste0("estimate_score.gct"),header = F,skip=2,stringsAsFactors = FALSE) 2 3estimate_score <- t(estimate_score) 4rownames(estimate_score) <- estimate_score[,1] 5colnames(estimate_score) <- estimate_score[2,] 6estimate_score <- estimate_score[-c(1:2),] 7estimate_score[1:4,1:4]

读进来的被字符串化了：

1 Description StromalScore ImmuneScore ESTIMATEScore 2AAACATACAACCAC "AAACATACAACCAC" "269.574325901855" "1068.73943272047" "1338.31375862232" 3AAACATTGAGCTAC "AAACATTGAGCTAC" "157.640476166265" "1245.7269262377" "1403.36740240396" 4AAACATTGATCAGC "AAACATTGATCAGC" "399.842094329677" "1208.89031399533" "1608.732408325" 5AAACCGTGCTTCCG "AAACCGTGCTTCCG" "595.490469453362" "1418.44175293577" "2013.93222238914"

所以AddMetaData要处理一下：

1pbmc<-AddMetaData(pbmc,metadata = estimate_score,col.name = colnames(df)) 2pbmc@meta.data$StromalScore <- as.numeric(as.vector(pbmc@meta.data$StromalScore)) 3pbmc@meta.data$ImmuneScore <- as.numeric(as.vector(pbmc@meta.data$ImmuneScore)) 4pbmc@meta.data$TumorPurity <- as.numeric(as.vector(pbmc@meta.data$TumorPurity)) 5pbmc@meta.data$ESTIMATEScore <- as.numeric(as.vector(pbmc@meta.data$ESTIMATEScore))

可以看到给每个细胞都打了分，bulk的一个样本是一个组织，可以用肿瘤纯度，单个细胞也还是有的吗？所以，是不是可以用某一群的平均表达谱来做呢？其实看到这只是根据基因列表的打分机制，这也未尝不可。

1head(pbmc@meta.data) 2 orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters Description 3AAACATACAACCAC pbmc3k 2419 779 3.0177759 1 1 AAACATACAACCAC 4AAACATTGAGCTAC pbmc3k 4903 1352 3.7935958 3 3 AAACATTGAGCTAC 5AAACATTGATCAGC pbmc3k 3147 1129 0.8897363 1 1 AAACATTGATCAGC 6AAACCGTGCTTCCG pbmc3k 2639 960 1.7430845 2 2 AAACCGTGCTTCCG 7AAACCGTGTATGCG pbmc3k 980 521 1.2244898 6 6 AAACCGTGTATGCG 8AAACGCACTGGTAC pbmc3k 2163 781 1.6643551 1 1 AAACGCACTGGTAC 9 StromalScore ImmuneScore ESTIMATEScore TumorPurity 10AAACATACAACCAC 269.5743 1068.739 1338.314 0.6956758 11AAACATTGAGCTAC 157.6405 1245.727 1403.367 0.6887845 12AAACATTGATCAGC 399.8421 1208.890 1608.732 0.6666206 13AAACCGTGCTTCCG 595.4905 1418.442 2013.932 0.6211327 14AAACCGTGTATGCG 540.3476 1019.445 1559.792 0.6719582 15AAACGCACTGGTAC 366.6339 1103.304 1469.938 0.6816675

可视化一下,看看不同分群下，"StromalScore" "ImmuneScore" "ESTIMATEScore" "TumorPurity" 四个分数的比例：

1p1<- RidgePlot(pbmc,features = "StromalScore") 2p2<- RidgePlot(pbmc,features = "ImmuneScore") 3p3<- RidgePlot(pbmc,features = "TumorPurity") 4p4<- RidgePlot(pbmc,features = "ESTIMATEScore") 5 6 # p4<- Seurat::CombinePlots(c(p1 ,p2,p3,p4)) 7 8library(gridExtra) 9grid.arrange(p1 ,p2,p3,p4,ncol = 2 )

就TumorPurity 来绘制熟悉的箱型图：

1ggplot(pbmc@meta.data, aes(x=RNA_snn_res.0.5, y=TumorPurity,fill=RNA_snn_res.0.5)) + 2 geom_boxplot(notch=TRUE)

严格的话，可以做一下显著性检验啊，这里我们就不过了。我比较好奇的是这四个是什么关系？为什么还有TumorPurity ，ESTIMATEScore？

1library(GGally) 2ggpairs(pbmc@meta.data[,c("StromalScore","ImmuneScore","TumorPurity","ESTIMATEScore","seurat_clusters")], 3 mapping = aes(color = seurat_clusters))

可以看到TumorPurity ，ESTIMATEScore是负相关的关系，其实知道一个就可以了。结合这些可视化的结果可以为我们了解哪些群的肿瘤纯度如何，从这个侧面来解释细胞的异质性。

那么有没有其他软件呢？显然是有的，一个可以用的就是Xcell了：https://github.com/dviraran/xCell

xCell is a gene signatures-based methodlearned from thousands of pure cell types from various sources. xCell applies a novel technique for reducing associations between closely related cell types. xCell signatures were validated using extensive in-silico simulations and also cytometry immunophenotyping, and were shown to outperform previous methods. xCell allows researchers to reliably portray the cellular heterogeneity landscape of tissue expression profiles.

还有一个：https://cibersortx.stanford.edu/，需要注册。

官网：
https://bioinformatics.mdanderson.org/public-software/estimate/

References

[1] Inferring tumour purity and stromal and immune cell admixture from expression data: https://links.jianshu.com/go?to=http%3A%2F%2Fwww.nature.com%2Fncomms%2F2013%2F131011%2Fncomms3612%2Ffull%2Fncomms3612.html
[2] 利用表达数据计算基质打分与免疫打分进而预测肿瘤纯度 --- ESTIMATE: https://www.jianshu.com/p/ec5307256ca5
[3] https://github.com/dviraran/xCell
[4] https://cibersortx.stanford.edu/
[5] https://bioinformatics.mdanderson.org/public-software/estimate/

基于Seurat结果推断单细胞群肿瘤纯度之ESTIMATE

References

相关推荐