基于非负矩阵分解的单细胞降维聚类分群

生信技能树

发布于 2022-03-03 13:09:16

2.3K0

发布于 2022-03-03 13:09:16

文章被收录于专栏：生信技能树

最近学徒在每天的单细胞分享腾讯会议上面讲解了一个2021发表在CELL杂志的脑胶质瘤的单细胞文献，标题是：《Inhibitory CD161 receptor identified in glioma-infiltrating T cells by single-cell analysis.》，其数据的链接在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163108 ，公开可以获取全部的表达量矩阵文件：

GSE163108_RAW.tar 135.4 Mb TAR (of CSV, H5)
GSE163108_metadata_10x.csv.gz 526.2 Kb   CSV
GSE163108_metadata_Smart-Seq2.csv.gz 98.2 Kb   CSV

里面关于CD4和CD8的T细胞的细分亚群蛮有意思的，如下所示：

CD4和CD8的T细胞的细分亚群

可以看到，在CD4和CD8的T细胞的各自矩阵内部降维聚类分群，这6个细分亚群都并不是泾渭分明的界限。听完分享才知道，原来作者这个时候的细分亚群其实并不关心它们内部是不是有不同的独立的单细胞亚群，仅仅是有这6个不同状态或者说发挥不同功能单细胞亚群。而区分它们的手段是非负矩阵分解，并不需要有很清晰的界限，只需要各个亚群的核心功能基因集有差异即可。

我们仍然是以 pbmc3k 数据集为例子给大家展现一下基于非负矩阵分解的单细胞降维聚类分群；

library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")  
data("pbmc3k")  
sce <- pbmc3k.final   
library(Seurat)
table(Idents(sce))
p1=DimPlot(sce,label = T)

因为pbmc3k 数据集的细胞种类比较多，我们这里仅仅是挑选单核单细胞亚群进行后续非负矩阵分解的演示，代码如下所示：

sub_sce = sce[,grepl('Mono',Idents(sce))]
dim(sub_sce)
table(Idents(sub_sce))
sub_sce$celltype = Idents(sub_sce)
p2=DimPlot(sub_sce,label = T)
library(patchwork)
p1+p2

可以看到，挑选单细胞子集的情况：

> table(Idents(sce))

 Naive CD4 T Memory CD4 T   CD14+ Mono            B 
         697          483          480          344 
       CD8 T FCGR3A+ Mono           NK           DC 
         271          162          155           32 
    Platelet 
          14 
> table(Idents(sub_sce))

  CD14+ Mono FCGR3A+ Mono 
         480          162

也可以可视化对比：

从pbmc3k 数据集挑选单核细胞

接下来的分析就是针对单核单细胞亚群矩阵。

先找各个亚群的特征基因

需要首先对单核单细胞亚群矩阵进行归一化，代码如下所示：

sub_sce=CreateSeuratObject(
  counts = sub_sce@assays$RNA@counts,
  meta.data = sub_sce@meta.data
)
library(dplyr)
sub_sce = NormalizeData(sub_sce) %>% FindVariableFeatures() %>% ScaleData(do.center = F)

这样得到了归一化的矩阵，进行后续非负矩阵分解分析，直接使用NMF包即可：

suppressPackageStartupMessages(library(NMF))
vm <- sub_sce@assays$RNA@scale.data
res <- nmf(vm,rank=2,method = "snmf/r") 
# 默认交替最小二乘法(Alternating Least Squares(ALS))——snmf/r  
# 参数rank=2，是期望的细胞亚群数量

可以看到，使用方法超级简单，就是NMF包的nmf函数，因为前面提到了这个单核细胞就是 CD14+ Mono 和FCGR3A+ Mono ，所以这里我们设置参数rank=2，是期望的细胞亚群数量，这个NMF包的nmf函数针对我们的矩阵进行了非负矩阵分解分析，得到了一个NMFfit的对象，里面的元素超级多。我们简单的看看：

plot(t(coef(res) ), 
     col=sub_sce$celltype)
head(basis(res))

可以看到，这个NMFfit的对象里面有我们的单细胞的降维坐标啦，而且很明显的CD14+ Mono 和FCGR3A+ Mono区分的很好。然后也有每个细胞亚群的各自基因权重占比：

> head(basis(res))
                    [,1]       [,2]
ISG15        0.073035389 0.09687721
B3GALT6      0.007342527 0.02587450
PUSL1        0.020592560 0.01623260
RP11-345P4.9 0.006339823 0.02125048
CDK11B       0.017043189 0.01725070
C1orf86      0.039707480 0.05986485

接下来我们就需要提取每个细胞亚群的权重排名靠前的基因，这里我们使用extractFeatures函数，也可以自己对上面的 basis(res) 筛选。：

# 期望单细胞分成2群，拿出来每个单细胞亚群各自的特征基因
fs <- extractFeatures(res,10L)
fs <- lapply(fs,function(x)rownames(res)[x])
 
> fs
[[1]]
 [1] "FCGR3A"   "RHOC"     "MS4A7"    "IFITM2"  
 [5] "CD79B"    "ABI3"     "SIGLEC10" "CTSC"    
 [9] "CKB"      "IFITM3"  

[[2]]
 [1] "S100A8" "S100A9" "LGALS2" "LYZ"    "GSTP1" 
 [6] "FCN1"   "RPL19"  "GNB2L1" "FOS"    "RPL12"

有这两个特征基因列表，就可以去可视化：

library(ggplot2) 
th=theme(axis.text.x = element_text(angle = 45, 
                                    vjust = 0.5, hjust=0.5)) 
DotPlot(sce[,grepl('Mono',Idents(sce))],
        features =   fs)  + coord_flip() +th

如下所示的可视化也可以看到：

DotPlot

然后降维聚类分群可视化

前面的非负矩阵分解相当于是替代了PCA操作，但是它的结果需要导入到seurat对象里面。所以进行如下所示的代码：

## 选择后续分析的因子，使用 NMF 运行的结果进行降维和聚类 
sub_sce
sub_sce <- RunPCA(sub_sce,verbose = F)
sub_sce@reductions$nmf <- sub_sce@reductions$pca
sub_sce@reductions$nmf@cell.embeddings <- t(coef(res) )
sub_sce@reductions$nmf@feature.loadings <- basis(res) 
sub_sce <- RunUMAP(sub_sce,reduction = "nmf",dims = 1:2)

我们的RunUMAP函数是基于非负矩阵分解后的结果哦，接下来进行分群：

sub_sce <- FindNeighbors(sub_sce,reduction = "nmf",dims = 1:2) %>% FindClusters(resolution = 0.1)

sub_sce$cluster <- apply(NMF::coefficients(res)[1:2,],2,which.max)

table( sub_sce$celltype ,sub_sce$cluster)
table(Idents(sub_sce) ,sub_sce$cluster)

可以看到我们这里有3种分群结果，而且一致性都挺好的：

> table( sub_sce$celltype ,sub_sce$cluster)
              
                 1   2
  CD14+ Mono    20 460
  FCGR3A+ Mono 158   4
> table(Idents(sub_sce) ,sub_sce$cluster)
   
      1   2
  0   0 383
  1 178   1
  2   0  80

最开始的CD14+ Mono 和FCGR3A+ Mono毫无疑问是金标准，然后我们的非负矩阵分解指定区分了两个亚群，最后基于非负矩阵分解的结果重新进行FindNeighbors和FindClusters根据resolution = 0.1又区分成为了3个亚群。

可以看到 FCGR3A+ Mono是比较纯粹的单细胞亚群了，但是 CD14+ Mono 是可以有细分的空间的。这些都是对自己的单细胞数据的进一步的理解。

非负矩阵分解的其它应用

从上面的演示来看，我们的基于非负矩阵分解的单细胞降维聚类分群特殊性在于，预先就指定了待分解的单细胞亚群数量，而且可以找到每个单细胞亚群的各自的特征基因，而无需走常规的降维聚类分群流程。

基于这个特性，我们的非负矩阵分解还有另外一个应用，也是在很多肿瘤单细胞文献里面可以看到，绝大部分的肿瘤研究单细胞研究我介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这个第一次分群规则是：

immune (CD45+,PTPRC),
epithelial/cancer (EpCAM+,EPCAM),
stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)

这样拿到了上皮细胞我们通常是走inferCNV去看其恶性与否，这个时候就有一个难题，我们判断出来了很多恶性的肿瘤细胞，但是它们其实是肿瘤的不同恶性程度，不同状态，虽然我们可以从算法是进行降维聚类分群，并且给出各个亚群的高表达量基因，但是它们会大量受肿瘤病人个体异质性的影响，因为如果不抹除病人特异性出来的结果就是各个病人的恶性肿瘤细胞独自成为一个亚群。

两年前的鼻咽癌单细胞文章：《Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma》，对其恶性的肿瘤上皮细胞就进行了非负矩阵分解，如下所示：

image-20220221144840145

尽管每个病人的肿瘤细胞都是大不一样，但是它们仍然是可以有生物学功能的共性，作者使用NMF包对11个肿瘤的7581个恶性上皮细胞进行处理：

We selected k = 4 as the number of factors because it yielded a high cophenetic correlation coefficient and effectively decomposed the dataset of each tumor.
A total of 44 metagenes were identified across the 11 tumors, we listed the top-ranked genes according to their loadings of the NMF factor Table S4.
The 44 metagenes were then compared by hierarchical clustering, using one minus the Pearson correlation coefficient over all gene scores as a distance metric. Five clusters of signatures were identified manually

也就是说这11个肿瘤病人，每个肿瘤病人独立的NMF处理，各自内部都是假定4个特征基因（metagenes），得到了 44个 metagenes，但是简单的相关性计算后层次聚类就可以看到其实是 5个基因集。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-02-21，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法