10X单细胞（10X空间转录组）数据分析之NMF寻找转录programs

hello，大家好，今天我们分享的内容来自文章Spatially organized multicellular immune hubs in human colorectal cancer,很多的内容都非常的经典，其中最重要的部分就是NMF对单细胞数据分析的影响，当然，也涉及到空间位置信息，很好的文章，我们把其中的重点提取出来，供大家参考。

前言部分

1、Recently, imaging-based studies have highlighted cellular interaction networks based on recurrent co-localization of different cells in neighborhoods (说白了还是要做临近通讯，而不是单细胞通讯的分析方法）。

2、revealed multicellular interaction networks based on co-variation of gene program activities（这个地方是今天的重点） in different cell subsets across individual tumors and imaged key molecules for predicted cell subsets and programs to localize these interaction networks in matched tissues from affected individuals.（空间网络与临近）。

单细胞分析部分

质控部分需要注意的地方，

（1）dropletUtils去除empty cell droplets，by testing against a background generated from barcodes with 1,000 to 10 UMIs, with cutoffs determined dynamically based on channel-specific characteristics.

(2) UMI and gene saturation was estimated in individual cells by sub-sampling reads without replacement in each cell barcode, in incremental fractions of 2%, with 20 repeats. A saturation function of the form y = ax/(x + b) + c was fit based on the number of UMIs observed while sampling reads at different depths.(这个地方对测序深度的计算跟10X官网的计算方法是类似的，目的在于尽量测到尽可能多的基因)。

（3）1、Fewer than 200 genes; 2、Fewer than 1,000 read（这个过滤我们可以忽略）; 3、Fewer than 500 UMIs; 4、 More than 50% of UMIs mapping to the mitochondrial genome（线粒体基因的阈值是50%）;5、 Non-empty droplet with false discovery rate (FDR) less than 0.1; 6、 Over 5% of reads estimated as coming from swapped barcodes/chimeric reads 这几个过滤标准大家要注意了。

（4）对每个通道和每个通道细胞类型组合分别使用异常值排除，然后根据以下标准标记偏离中位数 > 2 个四分位距 (IQR) 的细胞，1、 log10(total transcript UMI), 2、 Fraction of barcode swaps, 3、 Gene saturation estimate,4、 UMI saturation estimate, 5、 Fraction of UMI supported by > 1 read（异常值排除的地方也值得大家借鉴）。

（5）如果细胞显着偏离拟合，则根据以下关系进一步标记细胞: 1、 Total reads versus total UMI, 2、 Total UMI versus log likelihood of being empty; 3、 Total UMI versus total number of genes. A cell was excluded if it was flagged by at least two of these criteria for epithelial and immune cells, or at least three criteria for stromal cells.（这过滤标准真够严格的）。

Selection of variable genes, dimensionality reduction and clustering

高变基因的选择，这个地方考虑到样本的差异

降维，这里用NMF替代了PCA，NMF was used to reduce the dimensionality of the full dataset to between 15 and 40 dimensions as the product of two non-negative matrices these NMF components were used to calculate the k-nearest neighbors (k-NN) graph (k = 21) based on a cosine similarity distance. This graph was clustered using stability optimizing graph clustering， to identify 7 top level cell type clusters（当然，这里的单细胞数据没有用我们常用的seurat软件，而是作者自己的方法）。

批次效应的去除，也采用的直接算法而不是软件，To minimize differences across samples due to technical reasons (e.g., 10x v2 versus 10x v3), gene expression measurements of individual genes were quantile normalized, separately among cells of each top-level cellular compartment, such that the expression CDFs for each gene matched across all batches. Next, the same dimensionality reduction by NMF and graph-clustering procedure was applied iteratively to the transcriptomes of each top-level cell type separately, resulting in a total of 88 cell clusters spanning distinct types or states。（毕竟降维是NMF，很难直接嵌套现有的批次去除的软件）。

值得注意的是，基于 PCA 的 louvain 聚类导致质量相似的细胞子集定义（数据未显示）。 然而，由于我们通过 NMF 从头发现基因表达程序，因此决定始终使用 NMF 而不是 PCA 来定义cluster。

Cluster connectivity，cluster之间的相关性，采用了PAGA。PAGA edge thresholds were selected by using the minimum edge weight of the corresponding minimum spanning tree for each k-NN graph。这里的PAGA主要是用来分析肿瘤的cluster和正常细胞的cluster之间的相关性。

NMF运用的第二个重点，Identification of gene expression programs by NMF。有关NMF识别programs在我之前的文章10X单细胞（10X空间转录组）数据分析之细胞等级（分数），分享过具体过程，大家可以回顾一下。我们来看看文献怎么识别的。

We used as input the weight components matrices (W matrices)（数学的知识大家自行查阅，NMF矩阵降维的知识） from an NMF procedure that was run on 50-200 subsampled gene x cell subsets，We excluded outlier components by sorting components by their cosine distance to the 20th nearest neighbor and excluding components with unusually high distance by an elbow-based criterion.（我们通过components按到第 20 个最近邻居的余弦距离对components进行排序来排除离群值components，并通过基于肘部的标准排除具有异常高距离的components。 这里需要关注的是，输入的矩阵已经是原始矩阵经过NMF降维过的矩阵，而寻找program在降维过的矩阵上进行再次分析）接下来，我们构建了一个 k-NN 图 (k = 30)，并使用稳定性优化图聚类（关于图聚类的知识，之前也分享了很多了，大家可以自行查阅）识别该图中高度相似components的cluster（注意这里是识别高度相似的components），具有指数变化的尺度参数（0.1 到 10）。每个cluster中的components被中值平均为一个 components，从而形成“共识 NMF”components的候选名单。 这些被用作所有细胞和高度可变基因的第二轮 NMF 的初始化components矩阵（二次NMF）。上述程序分别应用于每个top-level细胞群和来自正常通道的上皮细胞。 对于每种细胞类型，这产生了 8 个解决方案，8-48 个cluster对应于分辨率参数的不同选择。 对于每种细胞类型，基于平均cluster轮廓的检查、残差图的变化以及通过手动检查输出程序中的top-level基因来选择单个解决方案。

这个地方跟之前分享的NMF寻找有效的programs是一致的，值得大家好好借鉴。

每个components的top 150个基因作为一个expression programs，这个地方跟之前的研究一致，只是选取基因的数量上有不同。

接下来就是对选取的program的运用，Identification of shared gene programs in malignant epithelial cells。

为了识别来自多个个体患者的恶性上皮细胞共享的表达programs，首先将上述一致的 NMF programs分别应用于来自每个患者的恶性细胞（来自肿瘤通道并分类为如上所述的恶性细胞）。 对于每个患者，生成了一个单独的共识 NMF 表达programs集 (Wmatrix)，programs的数量根据残差图自动选择。 接下来，将类似的共识方法应用于所有每名患者共识 NMF program集（所有 W 矩阵，每个患者一个）以及一组 17 个正常上皮programs（如上所述识别 - 基因表达的识别）的组合列表 NMF programs），以便在单个组合 NMF 解决方案中捕获恶性和正常上皮programs。 在这一共识聚类programs完成后，包括一个或多个正常上皮programs的 NMF cluster被排除在外，并使用相应的正常 NMF programs代替它们。

Calculating NMF transcriptional program activity

In order to calculate the NMF program activity matrix (H), we used non-negative least-squares (NNLS，关于NNLS，大家可以参考文章最优化算法-非负最小二乘-NNLS), solving the following equation for the matrix H, H = argminH>0|X - WH|F , given X andW, where H is the ‘program activity’ matrix, k is the by cell matrix; X is the gene by cell expression matrix, and W is the gene by k NMF expression program matrix. W was restricted to at most top 100 weighted genes per NMF component.In this way we can calculate the activity values for any cell including cells not part of the original NMF procedure used to discover the program ‘‘dictionary’’。数学的知识非常重要，大家如果有能力，钻研一下。

测试共变 NMF 表达programs，这个地方非常难，有大牛的话，希望讨论一下

NMF的运用远不止如此，希望以后多多学习这方面的知识

生活很好，有你更好