how to merge metaprograms with genenmf

[mhvbyrne]

Hi

I have two metaprograms which are similar and I wish to merge them. I don't want to discard them like the Yerly code does.

I want to create a new combined metaprogram which is re-ranked and re-weighted, so that I can use non negative least squares to derive a weighted signature score.

Please could you advise on how I can merge them appropriately with your package?

[mhvbyrne]

I've done some digging and I don’t think it is currently possible to merge in this way because I can’t find a full gene-by-metaprogram weight matrix in the saved object output from getmetaprograms.

metaprograms.genes.weights appears to contain truncated named vectors for each MP, derived from markers.consensus, rather than full consensus weights across the whole gene list across all MPs. From the code, it looks like markers.consensus is filtered by weight.explained, min.confidence, and max.genes before being stored.

This means that downstream merging is hard to do correctly, because genes missing from an MP’s stored vector may have been dropped during truncation rather than having true weight 0.

Would it be possible for getMetaPrograms() to optionally store a full consensus weight vector for each metaprogram, in addition to the current truncated representation? That would make post hoc merging and other downstream analyses much easier.

Please let me know if there's an alternative option.

[mhvbyrne]

What I have done to merge is take the union of genes in MP1 and MP2, and I have taken the average of their weights. The issue is because of the truncation it is likely there are some artificial zeros. e.g. for genes A B C D, MP1 has A B C present; MP2 has B C D present.

gene name...MP_A...MP_B
gene W.........0.5........0
gene X..........0.3........0.4
gene Y..........0.2........0.4
gene Z..........0...........0.2

The other issue with this approach is it will be down weighted by a factor relative to how many programs the gene is not in; e.g. merge MP_A; MP_B: MP_C, and gene has a weight/factor of MPA=0.3; MPB=0.1; MPC= score not present because of truncation. It would then get a score of (0.3+0.1+0)/3.

Whereas, if there is a merge of MPA, MPB, MPC, MPD, where MPD is also truncated; then it would be (0.3+0.1+0+0)/4

[mass-a]

Hi Matthew,
you are right that there isn't an implementation for merging and reweighing two meta-programs. Usually one would expect that two similar MPs should cluster next to each other, and changing the resolution of the hierarchical clustering (ask for one fewer MP) may naturally merge them. But I understand it may not always neatly fall into place like that.

However, there is a relatively simple way to extract the full gene x MP matrix. If you run the code by setting:

geneNMF.metaprograms <- getMetaPrograms(geneNMF.programs,
                                        weight.explained = 1,
                                        max.genes = 10^5,
                                        nMP=10)

You are telling the method to output all genes with non-zero coefficients (weight.explained=100%). You can then fill out all remaining genes for that specific MP as zeros, without any truncation. Does this help?

[mhvbyrne]

Thank you for replying so quickly.

Am I correcting in thinking the work flow could be as follows:

1. Run getMetaPrograms with standard parameters (e.g weight.explained 0.5; specificity.weight=5; min.confidence=0.5)
2. Select the metaprograms I wish to merge (e.g. MP1, MP2), based on QC, % gene overlap, GSEA overlap.
3. Repeat getMetaPrograms with weight.explained increased to 1 and max.genes to 10^5; all other parameters remain the same as step 1.
4. Then union merge MP1 and MP2 from step 4 output, add in zeros where genes are missing, take a mean of the weights of each gene.

Will the labels for MP1 and MP2 match in step 1 and 3? Or do I need to perform any additional checks/processing?

[mass-a]

Hi, you don't even have to run it twice, just steps 3 and 4.
The weight.explained and max.genes are a post-processing step, to give shorter MPs with the most informative genes. The MPs are calculated on the full gene vectors and will be identical regardless of the value of weigh.explained (just truncated to the top genes, and with gene weights rescaled to sum 1). You can verify this by running getMetaPrograms with the two settings.

[mhvbyrne]

Thank you for the advice, I have already run step 1 and 2 on a number of different cell compartments and annotated them with GSEA. I performed step 3 as you recommended however I can't figure out how to map this back onto the the output from step 2 I tried both programs.clusters and the MP names but I believe both of these change when you alter the settings. I am reluctant to repeat the QC and GSEA again on the whole output as it is a lot of work. I am planning to just use a union merge with the truncated outputs from step 2 for now.

If a future version could save the whole weighted output as well as the truncated output that would be ideal as then you could do the standard geneNMF workflow using the truncated outputs which is a real strength of the package, and perform QC on this, and then access the full weighted output if needed later easily.

[mass-a]

Thanks, that's a good suggestion. We could easily add an additional slot to the MP object with the full gene list and corresponding weights. Will keep that in mind.
Cheers
-m