In some cases it might be biologically meaningful (ha) to merge certain pathways nodes (e.g. protein complexes), as well as computationally helpful by summing low count vectors.

See section 2.2.2 and figure 2 in https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-13-20.

Meeting Notes:

1. better understand difference between separate and joint models (they are actually equivalent: https://books.google.ch/books?id=zyjWBgAAQBAJ&pg=PA137&lpg=PA137&dq=regression+for+to+separate+data+sets+indicator&source=bl&ots=OZiF7M0ShS&sig=ACfU3U3f5mni4Zj7-xY-RdvMsw8eVssHTQ&hl=de&sa=X&ved=2ahUKEwjahYC16e7oAhXM-KQKHR6fCdcQ6AEwAXoECAsQLw#v=onepage&q=regression%20for%20to%20separate%20data%20sets%20indicator&f=false)
2. try likelihood ratio test model (for single edge)
   • with delta vs without delta (with only one delta?)
3. log link function may be a better idea
4. try partial correlation with NB assumption
5. benchmark: set 90% of ground truth DCEs to 0 (makes setting more biologically relevant, maybe AUC can be used)
6. where to get DAGs from (KEGG, ...)
7. simulations: sampling beta and subtracting minimum biases beta (?)

graph1 <- create_random_DAG(n=p, prob=.05, lB=1)
Warning message:
In runif(length(negedges), min = lB[1], max = lB[2]) : NAs produced

In documentation it says lB is lower bound and uB is upper bound. I suppose you want to use that rather than lB[1], lB[2]

The simulations do probably not suffer much/at all from this. However, for real data we cannot just normalize the data, because we need counts for the glm.

How to deal with this problem?

Ideas:

• transform pathway to DAG (how?), how to validate
• adapt method (dynamical bayesian network??)

How to proceed, if samples in tissue < genes in pathway? Subsampling of genes? Prefiltering of genes, based on expression or position in pathway?

Is the somewhat complicated DAG creation in dce::create_random_DAG more or less similar to this much simpler implementation:

node_num <- 10
edge_prob <- .9
eff_min <- .2
eff_max <- 1.4

tmp <- matrix(rbinom(node_num * node_num, 1, edge_prob), node_num, node_num)
tmp[lower.tri(tmp)] <- 0
tmp[tmp != 0] <- runif(sum(tmp != 0), min = eff_min, max = eff_max)

dce::plot_network(tmp)

Computing Causal Effects using glm requires us to use a link function. Common choices are identity or log.

Cancer resources

• TCGA: https://portal.gdc.cancer.gov/
• TCGA barcode: https://docs.gdc.cancer.gov/Encyclopedia/pages/TCGA_Barcode/
• TCGA metadata: https://docs.cancergenomicscloud.org/docs/tcga-metadata/

Pathway resources

• KEGG: https://www.genome.jp/kegg-bin/network_viewer/

Modeling

• https://www.biorxiv.org/content/10.1101/599811v1.full#disp-formula-1

The following code generates warnings and an error if the value flip is applied.

A <- rnbinom(100, size=100, mu=1000)
B <- rnbinom(100, size=100, mu=0.1*A)

# value flip
A[1] <- 20
B[1] <- 2000

glm.nb(B ~ A, link="identity")

The error is no valid set of coefficients has been found: please supply starting values.

A subset of the warnings:

1: In log(y/mu) : NaNs produced
2: step size truncated due to divergence
3: In log(y/mu) : NaNs produced
4: step size truncated due to divergence
5: glm.fit: algorithm did not converge
6: In log(pmax(1, y)/mu) : NaNs produced
7: In log((y + .Theta)/(mu + .Theta)) : NaNs produced
[..]

Ground truth:
A -> B - > C

Perturbed:
A -> B - > C; A -> C

dce(A,C) is the same in both settings. However, the ground truth, would not compute a dce(A,C), since there is no edge. Is the dce(A,C) a false positive or not?

Idea 1

beta > 0 describes the relative change. 0.5 corresponds to halving and 2 to doubling the expression levels.
This is problematic because it requires a transformation of causal effects which is non-trivial (but possibly somehow doable?).

Idea 2

beta can be both positive and negative. Counts are propagated by multiplying beta with mean-standardized counts and adding noise.
This is problematic because standardizing might introduce artefacts and can lead to mu < 0 (which yields NaN counts).

beta <- -1.2

set.seed(42)
A.nb <- rnbinom(1000, size=10, mu=10)

B.nb <- beta * A.nb + rnbinom(1000, size=10, mu=10) # leads to negative counts
B.nb <- rnbinom(1000, size=10, mu=mean(A.nb) + beta * A.nb) # leads to negative mu, thus NA counts
B.nb <- rnbinom(1000, size=10, mu=10) + beta * scale(A.nb, scale=FALSE) # leads to negative counts
B.nb <- rnbinom(1000, size=10, mu=mean(A.nb) + beta * scale(A.nb, scale=FALSE)) # leads to negative mu, thus NA counts

MASS::glm.nb(B.nb ~ A.nb, link="identity")

Idea 3

Use a mean function for mu of rnbinom. This requires an appropriate link function during the regression.

beta <- -1.2

set.seed(42)
A.nb <- rnbinom(1000, size=10, mu=10)

B.nb <- rnbinom(1000, size=10, mu=exp(log(10) + beta * (A.nb - mean(A.nb)))) # link function keeps mu positive, exp can lead to extreme values

MASS::glm.nb(B.nb ~ A.nb, link="log")
glm(B.nb ~ A.nb, family=MASS::negative.binomial(theta=10, link="log"))
glm2::glm2(B.nb ~ A.nb, family=MASS::negative.binomial(theta=10, link="log"))

Consider the following graph:

A -> B -> C
|         ^
|_________|

Then the edge weight of A->C is not the total causal effect from A to C, but rather this total CE minus the CE of A->B->C

These tools shall be benchmarked next to ours.

Overview papers:

• https://academic.oup.com/bib/advance-article/doi/10.1093/bib/bbz158/5722384
• https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-3146-1

Real (e.g. KEGG) pathways may be missing true nodes and edges. Alternatively, they might have false ones.

We should investigate their effect using benchmarks.
The choice of valid adjustment set might influence this.

TPM seem to work as well or even better than raw counts with library size correction within the regression. Log normalization reduces accuracy.

benchmark_f1-score.pdf
benchmark_roc-auc.pdf

atm: negative binomial is hard coded

If theta and mu are constant for all genes in rnbinom, the actual theta propagated by the causal effects in he graph will mostly decrease. We can control for that by setting a lower mu and or higher theta for downstream genes.

I.e. produce p-value based on the dces for a pathway P.

Idea:
permutation test as in permute the node labels of the data and compute the differential causal effects. What is the statistic? absolute mean, absmax, abssum, ... dce?

Understand why KEGGgraph's gene splitting is wrong.
E.g. 3x HRAS

A few approaches from literature

Working with DE vs non-DE genes

• Choose real pathways (from GO) and randomly select some. Union of contained genes (+ noise) is labeled as significant. Then do FET for all pathways. (https://academic.oup.com/bioinformatics/article/22/13/1600/193669)

• Use real pathways (from MSigDB), randomly select one, and label its members as DE. (https://www.biorxiv.org/content/10.1101/659920v1)

Working with gene expression values

• Sample from normal distribution. Gene sets of varying differential expression and correlation patterns are simulated. (https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-10-47)

• Sample expression values from normal distribution. (https://www.nature.com/articles/s41540-017-0007-2)

• Sample fold changes from normal/uniform distribution. Then draw normalized estimated expression levels from normal distribution. (https://academic.oup.com/bioinformatics/article/25/2/211/218259)

• Log-expression values sampled from multivariate normal. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3458527/)

• Start from real expression data. Standardize data. Add signal to subsets. (https://link.springer.com/article/10.1186/s12859-019-3146-1)

• Draw expression from normal distribution. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5939912/)

Other simulation papers

• Draw expression level from gamma distribution. (https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btz752/5584234)

• Draw expression levels from negative binomial. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4635655/)

• Model interaction kinetics. Uses underlying network. (https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-7-43)

• Use stochastic simulations with delay. Uses underlying network. (https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-017-1731-8)

• Draw from negative binomial. (https://bioconductor.org/packages/release/bioc/vignettes/compcodeR/inst/doc/compcodeR.pdf)