getINITModel2

[fgualdr]

Hi,

I'm running through the tutorial (https://sysbiochalmers.github.io/Human-GEM-guide/gem_extraction/) to run getINITModel2.
It all works smoothly and nice. Anyhow I need to trick the setup in order to use the tool for specific cells and is not clear to me how should I design the inputs to the function to get per cell per tissue models.

Thanks a lot

Francesco

[JonathanRob]

Hi @fgualdr, I converted this issue into a discussion, since it's more of a question than a software problem. Also, we think it's a nice way to try out the new "Discussion" feature on GitHub.

Regarding your question, could you clarify a bit more what you would like to do with tINIT? Are you wanting to make a cell-type (sub-tissue) specific GEM, or many GEMs from many RNA-Seq profiles, or something else?

[fgualdr]

Hi @jonathanRO, thanks for the reply and sorry you are right is a question.
What I am trying to model is several independent RNAseq.. are the same cells but in a time course post stimulation.
What I would like to use is not really indipendendent extraction but using all the data together.
Given that a gene across the time course I am doing might show its maximum or its minimum in other time point I would like to get a per-time model but relative to the other conditions... and as it seemed that some of the features in the getINITModel2 would fit this purpose.
Thanks a lot

[JonathanRob]

Ok, I think I understand your question.

It sounds like you are interested in the more traditional approach of tINIT, which scores the genes based on their expression level relative to other samples rather than a single threshold. The Human-GEM guide walks through an example where an expression threshold of 1 TPM is used to weight the genes, where genes with an expression level higher than 1 TPM will receive a positive weight (to stay in the model), and genes with an expression lower than 1 TPM will receive a negative weight (to be removed from the model).

If you would instead prefer to score your genes based on their expression relative to other samples (i.e., other time points in your case), then the arrayData structure passed to the getINITModel2 function should not contain a threshold field (or it should be empty: []). Instead, your arrayData structure should contain a matrix of expression levels (in the levels field), where rows correspond to genes and columns correspond to "tissues", which in your case are different samples (time points). If you leave the threshold field empty, then the tINIT function will instead score each gene based on its expression level in the specified tissue (time point) relative to its mean expression across all tissues (time points).

As yet another alternative, if you would rather score the genes relative to a single time point (e.g., prior to stimulation or at time zero) instead of relative to the mean value across all time points, then you can set arrayData.threshold equal to the vector of gene expression values at time zero.

[fgualdr]

Excellent thank you so much!
Extremely useful explanation.
If I supply z-scores to the scoreComplexModel would it still behave normally?
We have performed in parallel proteomics and RNAseq as we know there are discrepancy between RNA and protein expression. To combine both into the reaction score is there a trick?
Thanks again

Francesco

[JonathanRob]

No problem!

If I supply z-scores to the scoreComplexModel would it still behave normally?

By default, tINIT will score genes (via the scoreComplexModel function) using the formula: 5*log(level/threshold), where level is the value for each gene provided in the arrayData.levels field, and threshold is either a single value or vector of mean values as described in my previous comment. The algorithm by default also restricts the gene scores to a range of [-5, 10], so for example if a gene level is zero and threshold is 1, its score: 5*log(0/1) = -Inf would be adjusted to -5.

This is a long way of saying that z-scores would be problematic because arrayData.levels (or arrayData.threshold) cannot be negative, since it will cause problems for the log transformation. If you wanted to pass a Z-score into the algorithm without having to modify the code itself, you could instead just set arrayData.levels equal to the exponentiated Z-scores (exp(Z)) so they revert back to Z-scores within the algorithm. Not the most elegant solution, unfortunately.

We have performed in parallel proteomics and RNAseq as we know there are discrepancy between RNA and protein expression. To combine both into the reaction score is there a trick?

In theory, yes, but in practice, it's complicated. tINIT does have this feature built-in, where proteomic data is supplied as the hpaData input and transcriptomic as arrayData. It will then use the proteomic data when available, and only use the transcriptomic data as a backup to fill in for genes that have missing values in the proteomic data. Furthermore, the gene scores calculated from protein data can range from [-8, 20], meaning they are given more weight than genes with transcriptome-derived scores ([-5, 10]).

The reason it's complicated is because, as the name suggests, the hpaData input was designed specifically for data retrieved from the Human Protein Atlas (HPA) and its format is rather specific and non-intuitive.

So as I see it, you have a few options here:

1. Familiarize yourself with the hpaData structure format by looking at the parseHPA and scoreComplexModel function code, and then format your protein data in the same way (probably a major pain in the ass).
2. Pre-combine your protein and RNA-Seq data into arrayData.levels by setting the values equal to the protein Z-scores, and then replacing any missing values with the RNA-Seq Z-scores (effectively what the tINIT algorithm would do if you were able to get the hpaData input to work)
3. Pre-combine your protein and RNA-Seq data into arrayData.levels using some other approach of your choosing.

[manas-kohli]

Hi @JonathanRob, I wanted to just ask a few follow-ups regarding this. I do know that things have slightly changed with ftINIT but I still had a few questions:

1. The threshold used in lots of your papers for the TPM for bulk data is 1.0. Is there a reason that this threshold is used/was it selected after you guys determined it was generally optimal for most datasets? The reason I'm asking is that I have reproducibility issues when I run tINIT on different cohorts mostly because of batch effect issues with the transcriptomics. I know in the ftINIT paper you guys performed batch correction on your transcriptomics data to account for this but because the threshold is 1.0, it leads to some issues to do with things like comparing metabolic function and gene essentiality if I build GEMMs per sample/patient.
2. If I understand the above formula for the inclusion and exclusion of tINIT correctly if you don't have a strict threshold, it should assign half the samples a negative weight because log(level/threshold) will be negative when level < threshold (which is the mean). Are these too many reactions in general to exclude given half the samples will have a given reaction excluded?