Difference between gene_spectra_score.k_3.dt_0_1.txt and gene_spectra_tpm.k_3.dt_0_1.txt

[tujchl]

Hi,

I applied cNMF to analyze my scRNA-seq data, and I focused on two output files: gene_spectra_score.k_3.dt_0_1.txt and gene_spectra_tpm.k_3.dt_0_1.txt. According to the manual, these files should represent essentially the same information but on different scales (Z-score vs. TPM). However, I noticed a significant discrepancy between the top 50 genes from each file, with only 15 genes overlapping between them. I used the R code below to extract the top genes for each GEP. Could there be an issue with my approach?
##########################################
exprSet = read.table(file = "MC.Merge.gene_spectra_score.k_3.dt_0_1.txt", header = TRUE, sep = "\t")
OR
exprSet = read.table(file = "MC.Merge.gene_spectra_tpm.k_3.dt_0_1.txt", header = TRUE, sep = "\t")
exprSet = exprSet[,-1]
rownames(exprSet) = paste(rep("C", 3), seq(1, 3, 1), sep = "")
top_genes = apply(exprSet, 1, function(x){ names(sort(x, decreasing = TRUE))[1:50]})
###########################################
Additionally, I exported the counts file with only HVGs from Seurat before applying cNMF, assuming it would save computation time. Is this approach acceptable?

[OlivierBakker]

Hi,

Not an official answer and late to the party, but as I had the same issue and there is no answer: As I understand it the spectra_score is done on the per-gene centered and scaled TPM data, which means there is no difference between genes with different levels of expression, whereas the spectra_tpm is not normalized to gene expression level, which is why the ranking is so different. The issue with single cell is for low expression genes, z-scoring the TPM data may not be 100% correct, but it should work in a practical sense. So for gene ranking goals you want to use the spectra_score as far as I understand, but an offical conformation would be nice :)

[dylkot]

Hi, the two files should be interpreted differently.

• gene_spectra_score.k_3.dt_0_1.txt - should be used to rank marker genes. The top genes for each program in this file are the ones most strongly linked to the program relative to the other programs. It varies between negative values (more strongly associated with other programs) to positive values (most strongly associated with the program).
• gene_spectra_tpm.k_3.dt_0_1.txt - represents how much each gene would be expressed in units of TPM in a cell that only expressed that one program. It doesn't necessarily related directly to how good of a marker gene each gene is.

Yes, it is reasonable to run cNMF only on the HVGs. However, it is not necessary to only provide this in the counts matrix. cNMF computes HVGs internally and only runs on that. And then it actually relearns the program representation on all of the genes if they are available. So I actually recommend providing it the full counts matrix rather than just the HVGs. You can actually specify a specific list of HVGs if you want using the --genes-file argument and providing a file with 1 gene name per row.

I hope this helps