Mathematical description of indirect effects

[ekernf01]

Dear Dr. Wang,

Dictys is a beautiful and thorough piece of work -- thank you for your continued efforts. I don't have a software issue to report, just a docs request. What's the mathematical formula for the total effects returned by dictys network indirect, and what scale are those effects on, or what experiment would be needed to observe the effect? I am not sure but I have some guesses.

• The function returns inverse of $B = (I - \beta)$ where $\beta$ contains the coefficients from the model, possibly after being scaled and/or regularized.
• To predict expression in wild type, use $B^{-1}\alpha$.
• To predict total effects of knockout, knockdown, or overexpression, use $B^{-1}\alpha$ and scale entry $j$ of $\alpha$ up or down if TF $j$ was targeted.
• To find the right value of $\alpha_j$, solve $x_j = B_j^{-1}\alpha$ where $x_j$ is the expression of gene $j$ post-intervention, $B^{-1}_j$ is row $j$ of $B^{-1}$, and $\alpha$ is known except for the $j$ th entry.

I have not studied stochastic processes much so I cannot tell at all whether this is equivalent to enforcing $X_j=0$ while simulating the differential equation model out to infinity.

Best regards. Eric Kernfeld

[lingfeiwang]

Hi Eric,

Thank you for taking the time into our theoretical derivations.

In this paper we are hoping to capture some of the kinetics that cannot be captured with traditional causal inference models, where stochastic process came into play. In this setting, we do not consider hard perturbations that simply set gene expression to specific values, like the 'do' operator. Such hard perturbations may be possible in social sciences like how the federal reserve sets the interest rate, but not really how things work in the physical world.

We considered three types of soft perturbations which we didn't elaborate in the paper.

• Perturbation in the basal trascriptional rate $\alpha$. Its first point of action in the system should be the corresponding gene, which can then propagate to other genes. Examples include CRISPRa/i.
• Perturbation in the degradation rate $\lambda$ (Section 3.4, Sup File 1), such as RNAi.
• Perturbation in the strength of this gene's regulation of all other genes, i.e. $\beta_{ji}$ for all $i$. One example is some CRISPR outcomes such as missense mutations for a TF.

I don't have the results at hand, but as far as I remember they are all equivalent and give the same result we used in the paper. We didn't consider complete KOs because they are almost identical with case 3 above because both give $xj\beta{ji}=0$, except having RNA level=0 which can lead to slight normalization differences.

Also, all the raw network parameters such as $\alpha$ and $B$ are already inferred before computing the total effect network. We don't need to find their values now. We are only considering infinitesimal perturbations on $\alpha$ (equivalent with other two cases above), which gave the theoretical result that $B^{(\infty)}=B^{-1}$. (Because the system is linear the effects can scale to finite perturbations.)

From here, there were several extra steps to consider for practical purposes:

• Starting from $B^{-1}$, we divided the perturbation outcome (logFC) on other genes v.s. on this gene, getting $\beta{ji}^{(\infty)}/\beta{ii}^{(\infty)}$. See https://github.com/pinellolab/dictys/blob/cb84353e77b011c1e8caabd307f82d29a37b7438/src/dictys/network.py#L890.
• We considered normalization changes because the perturbation can affect total mRNA count and the relevant level of other genes even if they are not indirectly regulated by the gene of interest. See https://github.com/pinellolab/dictys/blob/cb84353e77b011c1e8caabd307f82d29a37b7438/src/dictys/network.py#L890.
• We also needed to apply stronger regularization if the solution is divergent for Lyapunov equation, because CUDA doesn't yet provide a function to solve it and we had to use ad hoc approximations. See https://github.com/pinellolab/dictys/blob/cb84353e77b011c1e8caabd307f82d29a37b7438/src/dictys/network.py#L941.
• Due to the expression-dependent variance estimation bias in scRNA-seq, perturbation outcome scale can be distorted by regulator and target gene expression, so we further (inversely) scaled $\beta{ji}^{(\infty)}/\beta{ii}^{(\infty)}$ by the unexplained expression variance of genes $i$ and $j$. See https://github.com/pinellolab/dictys/blob/cb84353e77b011c1e8caabd307f82d29a37b7438/src/dictys/network.py#L969.

We only used this final output as total effect network, which indicates the relative logFC of target gene expression caused by 1 relative logFC of regulator gene expression. Relative logFC is defined as the mean log expression level change divided by the standard deviation of the stochastic noise level in log expression.

I hope that answers your question and don't hesitate to follow up.

Best, Lingfei

[ekernf01]

Thank you very much. This is extremely helpful. I do not fully understand the details yet, but I think I can proceed knowing just the last part of your response: relative FC in target per relative FC in regulator.

[ekernf01]

A couple of follow up questions now that I have read more:

• In your second bullet point, would you still do this when comparing to a validation experiment where scran or TMM normalization was used?
• In the last point, can you explain more about how this bias works?
• Can dictys.network.normalize be applied on either direct-effect or total-effect edge weights as needed?

Thanks again!

[lingfeiwang]

Sure!

• TMM normalization assumes most genes are not DE but it's possible they are weakly DE but insignificant with say 3 replicates. I image the best comparison would be not doing TMM with your data, or implement TMM in Dictys. But skipping the second bullet point appears a reasonable alternative too.
• Several papers mentioned this bias including Fig 3b in my Normalisr paper https://www.nature.com/articles/s41467-021-26682-1. When expression variance estimation is biased, estimations of the logFC in DE or effect size in GRN would be biased. Without this normalization, we observed edge strength to be dependent on expresison level, which is consistent with expression-dependent variance estimation bias. The dependency was gone after this normalization.
• Yes. By default, the direct and total effect networks already included these normalizations. You can also extract intermediate files without the normalization or apply partial normalizations. See https://github.com/pinellolab/dictys/blob/cb84353e77b011c1e8caabd307f82d29a37b7438/src/dictys/scripts/makefiles/common.mk#L158 and https://github.com/pinellolab/dictys/blob/cb84353e77b011c1e8caabd307f82d29a37b7438/src/dictys/scripts/makefiles/common.mk#L164

[github-actions[bot]]

Stale issue message