内容复习----NMF在单细胞空间中的运用

原创

追风少年i

发布于 2025-02-13 15:49:34

8700

代码可运行

运行总次数：0

代码可运行

作者，Evil Genius

大家觉得对于单细胞空间的数据分析，掌握的方法越多越好么？

如果你是初学者，这个问题的答案是，对；但如果你是有很多分析经验的科研工作者，这个问题的答案是很片面。

举一个简单的例子，单细胞主流的通讯分析软件有CCCExplorer,cell2cell,CellCall,CellChat,CellComNet,CellDialog,CellEnBoost,CellGiQ,CellPhoneDB,celltalker,CLARIFY,CommPath,COMUNET,CrossChat,CrossTalkeR,CytoTalk,DcjComm,DIISCO,exFINDER,FlowSig,GraphComm,iCELLNET,InterCellar,iTALK,LRLoop,MDIC3,MEBOCOST,mistyR,multinichenetr,NATMI,NeuronChat,NicheNet,NICHES,PyMINEr,RaCInG,RSoptSC,scCrossTalk,scDCA,scDiffCom,scHyper,scMLnet,scriabin,scTenifoldXct,scTensor,SEGCECO,SEnSCA,SingleCellSignalR,SPRUCE,TraSig，每一个都有高分文章引用。

空间主流的通讯分析软件有BulkSignalR,CCPLS,CellPhoneDB,COMMOT,Copulacci,DeepCOLOR,DeepLinc,DeepTalk,GCNG,Giotto,HoloNet,IGAN,ncem,NeST,NeuronChat,NICHES,Renoir,SCVA,SpaCET,spaCI,Spacia,SpaOTsc,spARC,SpaTalk,SpatialDM,STcomm,stMLnet,VGAECCI。

那是不是说明我们做单细胞空间项目每个软件都分析一遍然后解析结果呢？

这显然是不可能的，所以对于有经验的工作者，理解背后的原理，根据自己的课题有选择的进行分析才是正解。

好了，开始我们今天的复习，NMF讲过很多了，其意义不再过多赘述，核心原理如下：

这里补充一点，如下图。

如果我们寻找细胞的异质性，通过标记基因来注释亚群，那Aty和Cla是无法区别的，因为特征基因两者都表达，只是表达高低的不同，如果想区分，NMF是一个很好的办法。

空间也经常使用，cell2location就提供了NMF分析细胞类型关系关系的例子。

今日目标，寻找NMF因子以及识别因子关键的基因

代码很简单

library(Seurat)
library(ggplot2)
library(UCell)
library(patchwork)
library(tidyr)
library(dplyr)
library(RColorBrewer)

####install.packages("GeneNMF")
###remotes::install_github("carmonalab/GeneNMF")  #from GitHub
library(GeneNMF)

seu <- readRDS('test.rds')

Consistent NMF programs across multiple samples

识别稳健的基因程序需要在多个样本和输入参数变化中的检测。对于NMF来说，最重要的参数可能是维度k，其对应于低维矩阵的program数量。为了确定稳健的program，需要运行NMF多次，并确定在这些运行中始终存在的program。

DefaultAssay(seu) <- "RNA"
seu.list <- SplitObject(seu, split.by = "Sample")

geneNMF.programs <- multiNMF(seu.list, assay="RNA", k=4:10, min.exp = 0.05)

可视化构成meta-programs的单个基因programs之间的成对相似性(以余弦相似度或Jaccard指数衡量)

ph <- plotMetaPrograms(geneNMF.metaprograms,
                       similarity.cutoff = c(0.1,1))
ph

geneNMF.metaprograms$metaprograms.metrics
##      sampleCoverage   silhouette meanSimilarity numberGenes numberPrograms
## MP1       1.0000000  0.696882673          0.723          21             50
## MP2       1.0000000  0.297392424          0.433          33             60
## MP3       0.9090909  0.129408777          0.295         142             60
## MP4       0.9090909 -0.073756249          0.161          45            109
## MP5       0.8181818  0.509947198          0.624          23             45
## MP6       0.5454545  0.459199536          0.560          11             20
## MP7       0.5454545  0.448717232          0.543          19             17
## MP8       0.5454545  0.008033304          0.233          52             44
## MP9       0.3636364  0.549435393          0.617          45             12
## MP10      0.2727273  0.530347451          0.599          36             12

Interpretation of gene programs

lapply(geneNMF.metaprograms$metaprograms.genes, head)
## $MP1
## [1] "TOP2A"  "CENPF"  "NUSAP1" "PCLAF"  "MAD2L1" "TK1"   
## 
## $MP2
## [1] "KRT6B"   "CALML5"  "S100A7"  "DMKN"    "LGALS7B" "KRT1"   
## 
## $MP3
## [1] "NIN"    "FAT3"   "PEAK1"  "NAV2"   "SLC7A2" "INTS6" 
## 
## $MP4
## [1] "FOSB"   "FOS"    "EGR1"   "ZFP36"  "DNAJB1" "NR4A1" 
## 
## $MP5
## [1] "ACTA2" "TAGLN" "MYLK"  "MYL9"  "CSRP1" "ANXA3"
## 
## $MP6
## [1] "S100A8" "S100A9" "CHI3L1" "TIMP1"  "CCL2"   "S100A7"
## 
## $MP7
## [1] "ATF3"  "ZFP36" "SOCS3" "IRF1"  "BTG2"  "NR4A1"
## 
## $MP8
## [1] "CRCT1"     "CRNN"      "RBP2"      "LINC00964" "FOS"       "SOX18"    
## 
## $MP9
## [1] "SFRP1"   "ALK"     "COL14A1" "DCN"     "MEGF6"   "CHGA"
geneNMF.metaprograms$metaprograms.genes.weights$MP6
##      S100A8      S100A9      CHI3L1       TIMP1        CCL2      S100A7 
## 0.215830597 0.210983115 0.115345508 0.105439904 0.090343059 0.051177038 
##       TAGLN        MT2A       INHBA       KRT6A       KRT6B      S100A2 
## 0.037815897 0.029383116 0.023189284 0.022741619 0.019063225 0.013392001 
##      IFITM3       ZFP36       KRT75         MGP       SOCS3       ANXA3 
## 0.011160279 0.010551846 0.009011202 0.006894923 0.006609389 0.005611245 
##       KLF10     TMEM45A       SULF2 
## 0.005390284 0.005292033 0.004774439

因子解读

top_p <- lapply(geneNMF.metaprograms$metaprograms.genes, function(program) {
  runGSEA(program, universe=rownames(seu), category = "C5", subcategory = "GO:BP")
})

head(top_p$MP1)
##                            pathway         pval         padj overlap  size
##                             <char>        <num>        <num>   <int> <int>
## 1:         GOBP_CELL_CYCLE_PROCESS 2.071463e-27 1.585705e-23      27  1175
## 2:         GOBP_MITOTIC_CELL_CYCLE 9.610340e-26 3.678358e-22      24   866
## 3:                 GOBP_CELL_CYCLE 1.444483e-24 3.685839e-21      28  1716
## 4:          GOBP_ORGANELLE_FISSION 2.410380e-24 3.842722e-21      20   490
## 5: GOBP_MITOTIC_CELL_CYCLE_PROCESS 2.509942e-24 3.842722e-21      22   714
## 6:              GOBP_CELL_DIVISION 5.394203e-24 6.882104e-21      21   618
##    overlapGenes
##          <list>
## 1: ASPM, BR....
## 2: BRCA1, B....
## 3: ASPM, BR....
## 4: ASPM, BR....
## 5: BRCA1, B....
## 6: ASPM, BR....

Signature scores for gene programs

mp.genes <- geneNMF.metaprograms$metaprograms.genes
seu <- AddModuleScore_UCell(seu, features = mp.genes, ncores=4, name = "")

VlnPlot(seu, features=names(mp.genes), group.by = "patient_bcc",
        pt.size = 0, ncol=5)