Gene names from components in NMF analysis

Question

This question was also asked on Biostars

I am using Rcppml package in R for my NMF analysis, I have a matrix from single-cell analysis. I have cells scores and I wanted to know how to extract genes instead of the components.

i used the following codes.

t1.cpm<- LogNormalize(data = as.matrix(t1@assays$RNA@counts)) # Log converting the data form seurat obj

t1.model <- RcppML::nmf(as.matrix(t1.cpm), 30, verbose = F, seed = 1234)

t1.w <- t1.model$w
t1.d <- t1.model$d
t1.h <- t1.model$h

dim(t1.w)
rownames(t1.w)<- rownames(t1.cpm)
colnames(t1.w)<- paste0("component", 1:30)
rownames(t1.h)<- paste0("component", 1:30)

and my output is as follows

        T_1_TACCTATCAGATGGCA T_1_CCTATTATCAGGCAAG
component1             0.000000000            0.0000000000
component2             0.000000000            0.0023590277
component3             0.001512931            0.0000000000
component4             0.000000000            0.0021574559
component5             0.001175799            0.0020222144
component6             0.000000000            0.0004463346

It would be helpful if you can let me know how to extract Gene names instead of components.

Thanks, Dave.

---

Rcppml is a machine learning library.

https://github.com/zdebruine/RcppML

NMF is Non-negative matrix factorization which is used in bioinformatics and AI, which can examine the latent relationships in experimental data sets.

Answer 1

NMF of a gene expression matrix into $k$ components works to decompose the expression matrix into an activity matrix ($H$) and latent factor matrix ($W$).

$W$ describes the $k$ latent factors (often called components) by providing a score for how much each gene contributes to each component. $H$ describes how much each component contributes to the expression of each cell.

To extract the activity and component matrices from your NMF model fit using RcppML, use the following code:

# cells is a Seurat object
model <- RcppML::nmf(cells@assays$RNA@data, k = 5, verbose = F, seed = 1234)
W <- model$w # Components Matrix
H <- model$h # Activity Matrix

# Name genes/components for latent factor matrix
colnames(W) <- paste0("component", 1:k)
rownames(W) <- rownames(cells) # Gene names

# Name components/cells for activity matrix
colnames(H) <- colnames(cells) # Cell IDs
row names(H) <- paste0("component", 1:k)

> head(W)
              component1   component2   component3 component4   component5
MIR1302-2HG 0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
FAM138A     0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
OR4F5       0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
AL627309.1  2.927886e-07 9.248074e-07 7.961105e-07          0 2.802758e-06
AL627309.3  3.486284e-08 0.000000e+00 1.916799e-07          0 7.938089e-08
AL627309.2  0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00

> head(H)
                 cell_1       cell_2       cell_3       cell_4       cell_5
component1 0.0001370465 0.0002195795 0.0001628898 0.0002239623 0.0001484157
component2 0.0003604048 0.0002239471 0.0005688147 0.0005242678 0.0004821531
component3 0.0003948121 0.0001485325 0.0002885584 0.0002242272 0.0002865587
component4 0.0001391838 0.0002418495 0.0001361317 0.0001382889 0.0001558255
component5 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000

People will often use the top $N$ genes to describe a component. The following code gets the top $N$ genes for each component and stores them in a list

# This code assumes the above code block has been already run
require(dplyr) #used for %>% piping
N <- 100
top_genes <- list()
for (i in 1:dim(W)[2]) {
  name <- colnames(W)[i]
  top_genes[[name]] <- sort(W[, i], decreasing = T) %>% head(N) %>% names()
}