chuiqin / irGSEA

The integration of single cell rank-based gene set enrichment analysis
Other
106 stars 17 forks source link

irGSEA.density.scatterplot #1

Open TPOW-001 opened 2 years ago

TPOW-001 commented 2 years ago

您好, 近期在尋找single cell scoring相關方式, 發現了這麼方便的整合性package.

但目前使用上有一個問題是, 在irGSEA.density.scatterplot畫出的scatter plot如圖呈現: https://ibb.co/GfZYsnY

scatterplot <- irGSEA.density.scatterplot(object = da1, method = "ssgsea", show.geneset = c("MSC-1","MSC-2","EC-1"), reduction = "umap")

但如果直接調出ssGSEA data的話, 會如下圖呈現: https://ibb.co/pbkrL8k

DefaultAssay(da1) <- 'ssgsea' FeaturePlot(da1, features = c("MSC-1","MSC-2","EC-1"))

想詢問您在這兩張圖是如何進行換算的, 謝謝~

再次感謝您整合了這麼多種scoring方式, 方便讓我們對有興趣的signatures進行比較.

chuiqin commented 2 years ago

很感谢你使用我的工具。irGSEA.density.scatterplot函数反映的是基因集打分分数的密度分数。它主要是通过Nebulosa包实现的,Nebulosa包结合每个细胞在二维空间的投影(例如TSNE,UMAP)以及细胞自身的基因集打分分数,估算了一个加权的核密度估计分数。这个密度分数就是irGSEA.density.scatterplot图中展示的分数。然后,FeaturePlot函数主要反映的是基因集打分的实际分数。ssgsea的打分分数看着有点大,是因为我取消了ssgsea中最后一不的标准化,为的是让这个基因集富集分析算法不受样本组成和批次的影响。其余几种算法(例如AUCell,UCell和singscore)也是基于单个细胞基因表达秩次排名计算。更多关于irGSEA包算法原理以及每一步实现的方法的中文教程,可以参考我在简书上的教程:https://www.jianshu.com/p/463dd6e2986f。 最后,感谢你的提问,我也在最近的1.1.2版本中修复了irGSEA.density.scatterplot函数(主要是取消不区分基因正负向情况下,singscore分数的center化。这种center化不利于irGSEA.density.scatterplot函数的展示。取消singscore分数的center化,对irGSEA.integrate的计算结果有轻微影响,需要重新计算)和irGSEA.density.scatterplot函数(主要是调整了可视化中的参数选择,获得更好的可视化体验)中存在的一些小问题。假如你想获得更好的irGSEA.density.scatterplot函数可视化体验,可以更新到1.1.2版本。 下面3个代码块可以达到一样的可视化效果:

library(irGSEA)
irGSEA.density.scatterplot(object = pbmc3k.final, 
                                          method = "ssgsea",  
                                          show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE", 
                                          reduction = "umap")
library(viridis)
library(Nebulosa)
library(ggplot2)
DefaultAssay(pbmc3k.final) <- "ssgsea"
test <- FeaturePlot(pbmc3k.final, features = "HALLMARK-INFLAMMATORY-RESPONSE",slot = "scale.data", reduction = "umap")
test <- test$data
test$HALLMARK.INFLAMMATORY.RESPONSE <- as.numeric(test$HALLMARK.INFLAMMATORY.RESPONSE)
dens <- Nebulosa:::calculate_density(test[, 4], test[,1:2], method = "wkde")
test$dens <- dens
ggplot(test, mapping = aes(x = UMAP_1, y = UMAP_2, colour = dens)) +
  geom_point(size = 0.5)+
  scale_color_viridis()+
  theme_classic()+
  ggtitle("HALLMARK-INFLAMMATORY-RESPONSE")+
  theme(plot.title = element_text(hjust = 0.5))
library(Nebulosa)
library(ggplot2)
DefaultAssay(pbmc3k.final) <- "ssgsea"
Nebulosa::plot_density(pbmc3k.final,  features = "HALLMARK-INFLAMMATORY-RESPONSE",
                       slot = "scale.data", reduction = "umap",
                       method = "wkde", joint = T) +
  ggplot2::theme(plot.title = ggplot2::element_text(size = 10, hjust = 0.5),
                 axis.title = ggplot2::element_text(size = 10))
TPOW-001 commented 2 years ago

真的非常謝謝您這麼詳盡的回覆, 為了理解Neublosa的作用, 所以花了幾天理解Weighted 2D kernel density estimation, 也條 列一下目前使用您這麼方便使用的irGSEA心得:

  1. 中間發現了使用PCA帶入座標的話, 會造成跟UMAP以及tSNE完全不同的結果, 觀測的gene sets/signatures會表現在不 同細胞族群中. Montage

去看了一下Neublosa作者發在Bioinformatics的paper (PMID:33459785), 結果在sup Fig 1A-1F也會看到相同現象. image

目前僅能猜測PC_1, PC_2的座標並不能代表整體細胞狀況, 也因此tSNE, 以及UMAP這種比較可以代表更多dimension(包含更多基因數)的呈現方式, 會有較好的呈現效果. 也因此FindNeighbors包含更多dimensions應該會更準確.

2. 另外, 也借用您的code以及idea, 測試microarray以及bulk RNA-seq是否能呈現特別gene sets/signatures在特定細胞上. 目前從http://biocc.hrbmu.edu.cn/CellMarker/ 網站參考取出Hematopoietic stem cell (HSC), Monocyte (MO), Lung epithelial cell (EpiC), Endothelial cell (EC)的相關gene set. 的確也能在microarray及bulk RNA-seq看到各自細胞的signatures出現在各自細胞上. (目前測試: Microarray要RMA格式, bulk RNA-seq要TPM格式) Montage

3. 目前呈現最好狀況的是ssgsea的方式, 其他計分方式使用範例細胞去做的話, 似乎會出現一定的偏移 (例如HSC signatures不一定在HSC細胞族群中表現)

4. 覺得目前您開發的irGSEA非常好用, 假設能整合microarray以及bulk RNA-seq資料也能同時應用的話, 覺得會更加潛力無限. 以下是自己使用microarray data應用的相關資料, 也許能讓您參考: https://1drv.ms/u/s!Ali-uEPmqcGojcJPXPRmZzGJUHdU_g?e=ZdedD0