sva_simulation.Rmd

---
title: "SVA Simulation"
author: "Christopher Lo"
date: "9/6/2022"
output:
  pdf_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r echo=F, message=F, warning=F}
library(sva)
library(tidyverse)
library(pheatmap)
library(DESeq2)
theme_set(theme_bw())
```


```{r echo=F}
#useful functions
sum_of_span_residue = function(X, Y) {
  X = as.matrix(X)
  Y = as.matrix(Y)
  stopifnot(nrow(X) == nrow(Y))
  #Z-score X and Y by columns
  X = scale(X)
  Y = scale(Y)
  mod1 = lm.fit(X, Y)
  mod1_res = norm(mod1$fitted.values - Y, type = "F") / (nrow(Y) * ncol(Y))
  mod2 = lm.fit(Y, X)
  mod2_res = norm(mod2$fitted.values - X, type = "F") / (nrow(X) * ncol(X))
  return(mod1_res + mod2_res)
}

myPairs = function(df1, df2, bottom = "", left = "") {
  stopifnot(nrow(df1) == nrow(df2))
  plots = list()
  k = 1
  for(i in 1:ncol(df1)) {i
    for (j in 1:ncol(df2)) {
      correlation = round(cor(df1[, i], df2[, j]), 2)
      df = data.frame(x = df1[, i], y = df2[, j])
      plots[[k]] = ggplot(df, aes(x, y)) + geom_point(size = .8) + ggtitle(correlation) +
        theme(axis.text.x=element_blank(),
              axis.text.y=element_blank(),
              axis.ticks=element_blank(),
              legend.position="none",
              axis.title.x=element_blank(),
              axis.title.y=element_blank())
      k = k + 1
    }
  }
  top = paste0("SV Spanning Score: ", round(sum_of_span_residue(df1, df2), 4))
  print(marrangeGrob(plots, nrow = ncol(df1), ncol = ncol(df2), top = top, bottom = bottom, left = left))
}
```

Simulation studies inspired from
*"A general framework for multiple testing dependence"* (Leek et al. 2008)

# Simulation Set-Up, one single experiment

We generate $X$ from the following model: $$X = BS + \Gamma G + U$$ 

We have $m = 1000$ genes (tests), $n = 20$ samples, and $r = 2$ latent variables.

Sampling noise: $U_{m,n} \sim N(0, 1)$.

The design matrix $S$ is 10 cases and 10 controls: $S_{1, n} = 1$ for $n=1:20$. Then, $S_{2, n} = 0$ for $n = 1:10$, $S_{2, n} = 1$ for $n = 11:20$.

Control effect for all genes: $b_{m,1} \sim N(0, 1), m=1:1000$

Case effect for DE genes $m=1:300$: $b_{m, 2} \sim N(3, 1)$

Case effect for Non-DE genes $m=301:1000$: $b_{m, 2} = 0$

Latent design matrix  $G$: $G_{r, n} \sim Bernoulli(.4), n=1:10$. 

$G_{r, n} \sim Bernoulli(.6), n=11:20$, where $r = 1, 2$. 

Latent effect 1: $\Gamma_{m, 1} = 0, m=1:300$,  

$\Gamma_{m, 1} \sim N(1, 1), m=301:1000$. (Positive signal overlaps with Non-DE genes, will lead to FPs if not corrected)

Latent effect 2: $\Gamma_{m, 2} \sim N(-1, 1), m=1:300$, $\Gamma_{m, 2} = 0, m=301:1000$. (Negative signal overlaps with DE genes, will lead to FNs if not corrected)

Therefore, for every gene, whether it is DE or not, it will be affected by one of the two latent variables 

**To ask/consider: **

-   Currently we have negative expression due to Latent effect 2.

-   Should we normalize before running analysis?


```{r echo=F}
set.seed(2023)
m = 1000 #number of genes (tests)
n = 20 #number of samples
r = 2
DE_effect = 3
Latent_effect = 1
nSims = 100

#for(i in 1:nSims) {
  U = matrix(rnorm(m * n, mean = 0, sd = 1), nrow = m)
  
  b1 = rnorm(m, mean = 0, sd = 1) #intercept
  b2 = mapply(mu = c(rep(DE_effect, 300), rep(0, m - 300)),
              sigma = c(rep(1, 300), rep(0, m - 300)),
              function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  B = cbind(b1, b2)
  S = matrix(c(rep(1, n), rep(0, 10), rep(1, 10)), byrow = T, ncol = n)
  
  gamma1 = mapply(mu = c(rep(0, 300), rep(Latent_effect, m - 300)),
                  sigma = c(rep(0, 300), rep(1, m - 300)),
                  function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  gamma2 = mapply(mu = c(rep(-Latent_effect, 300), rep(0, m - 300)),
                  sigma = c(rep(1, 300), rep(0, m - 300)),
                  function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  Gamma = cbind(gamma1, gamma2)
  G = mapply(p = c(rep(.4, 10), rep(.6, 10)), #no need for intercept for kernel G.
             function(p) as.numeric(rbernoulli(2, p)))
  
  X = B %*% S + Gamma %*% G + U
  Xnull = B %*% S + U
  
  
#}
```

```{r echo=F}
case_control = data.frame(samples = c(rep(0, 10), rep(1, 10)))
DE_genes = data.frame(genes = c(rep(1, 300), rep(0, m - 300)))
colnames(Xnull) = case_control$samples
rownames(Xnull) = DE_genes$genes
pheatmap(Xnull, annotation_col = case_control, annotation_row = DE_genes)
colnames(X) = case_control$samples
rownames(X) = DE_genes$genes
pheatmap(X, annotation_col = case_control, annotation_row = DE_genes)
```


-   What is the correlation between true latent variables and primary variables?

Primary case/control vs. latent 1: `r cor(S[2 ,], G[1 ,])`

Primary case/control vs. latent 2: `r cor(S[2 ,], G[2 ,])`

Primary design matrix vs. latent design matrix span residual (not sure): `r sum_of_span_residue(t(S)[, 2], t(G))`


# Estimate the number of SVs:

```{r echo=T, warning=F, message=F}
n.sv = num.sv(X, t(S), method = "be")
cat("Number of SVs: ", n.sv, "\n")

pca = prcomp(t(X))
variance = pca$sdev^2 / sum(pca$sdev^2)
qplot(c(1:length(variance)), variance) + geom_line() + geom_point() + 
  geom_hline(yintercept=1/ncol(X), linetype = "dashed") +
  xlab("Principal Component") + ylab("Variance Explained") + ggtitle(paste0("Number of SVs: ", n.sv)) + ylim(0, 1)
```

# Estimate SVs, primary variable coefficients, and SV coefficients

-   Are latent variables are spanned by the estimated SVs?

-   Are the estimated coefficients similar to true coefficients?

-   Is the null p-value distribution uniform?

-   Do the ranks of top genes match?


```{r echo=T, warning=F, message=F}
nullMod = t(S)[, 1]
svobj = sva(X, t(S), nullMod, n.sv = n.sv)
```

Inferred SV vs. latent design matrix span residual (not sure): `r sum_of_span_residue(t(S)[, 2], svobj$sv)`.

Latent 1 vs. SV 1: `r cor(G[1 ,], svobj$sv[, 1])`

Latent 2 vs. SV 2: `r cor(G[2 ,], svobj$sv[, 2])`


```{r echo=T, warning=F, message=F}

nullmodsv = cbind(nullMod, svobj$sv)
modsv = cbind(t(S), svobj$sv)
fitsv = lm.fit(modsv, t(X))

#visually look at predicted coefficients
plot_df = data.frame(b1 = B[, 1], 
                     b2 = B[, 2], 
                     b1_hat = fitsv$coefficients[1 ,],
                     b2_hat = fitsv$coefficients[2 ,],
                     b2_labels = c(rep("alt", 300), rep("null", m - 300)),
                     gamma1 = Gamma[, 1],
                     gamma1_hat = fitsv$coefficients[3 ,],
                     gamma1_labels = c(rep("alt", 300), rep("null", m - 300)),
                     gamma2 = Gamma[, 2],
                     gamma2_hat = fitsv$coefficients[4 ,],
                     gamma2_labels = c(rep("alt", 300), rep("null", m - 300)))

ggplot(plot_df, aes(b1, b1_hat)) + geom_point() + labs(x = "True b_m1", y = "Est. b_m1")
ggplot(plot_df, aes(b2, b2_hat)) + geom_point() + facet_wrap(~b2_labels) + labs(x = "True b_m2", y = "Est. b_m2")
ggplot(plot_df, aes(gamma1, gamma1_hat)) + geom_point() + facet_wrap(~gamma1_labels) + labs(x = "True gamma_m1", y = "Est. gamma_m1")
ggplot(plot_df, aes(gamma2, gamma2_hat)) + geom_point() + facet_wrap(~gamma2_labels) + labs(x = "True gamma_m2", y = "Est. gamma_m2")
```

Not sure what's going on here yet regarding p-values and ranking.

```{r echo=F}
#just compute p-value of b_12 = 0 using F statistics. 
pValuesSV = f.pvalue(X, modsv, nullmodsv)
pValuesNoSV = f.pvalue(X, t(S), as.matrix(nullMod))

#double check with existing function:
pValuesSV2 = rep(NA, 1000)
fstat = rep(NA, 1000)
for(i in 1:1000) { 
  dat = data.frame(x = X[i ,], pv = S[2 ,], sv1 = svobj$sv[, 1], sv2 = svobj$sv[, 2])
  dat_nullmod = lm(x ~ sv1 + sv2, dat)
  dat_mod = lm(x ~ pv + sv1 + sv2, dat)
  an = anova(dat_nullmod, dat_mod)
  pValuesSV2[i] = an$`Pr(>F)`[2]
  fstat[i] = an$F[2]
}


plot(1:length(pValuesNoSV[301:1000])/(length(pValuesNoSV[301:1000])+1),sort(pValuesNoSV[301:1000]), main = "No SV control, No signal", xlab = "P-values", ylab = "Uniform(0, 1)", xlim = c(0, 1), ylim = c(0, 1))
abline(a = 0, b = 1)

plot(1:length(pValuesSV[301:1000])/(length(pValuesSV[301:1000])+1),sort(pValuesSV[301:1000]), main = "SV controlled, No signal", xlab = "P-values", ylab = "Uniform(0, 1)", xlim = c(0, 1), ylim = c(0, 1))
abline(a = 0, b = 1)

plot(1:length(pValuesNoSV[1:300])/(length(pValuesNoSV[1:300])+1),sort(pValuesNoSV[1:300]), main = "No SV control, Signal", xlab = "P-values", ylab = "Uniform(0, 1)", xlim = c(0, 1), ylim = c(0, 1))
abline(a = 0, b = 1)

plot(1:length(pValuesSV[1:300])/(length(pValuesSV[1:300])+1),sort(pValuesSV[1:300]), main = "SV controlled, Signal", xlab = "P-values", ylab = "Uniform(0, 1)", xlim = c(0, 1), ylim = c(0, 1))
abline(a = 0, b = 1)


FStat_order = order(fstat, decreasing = T)
true_order = order(abs(b2), decreasing = T)

qplot(FStat_order, true_order) + geom_smooth()
```


## "Knobs to turn" in estimating the number of SVs

$\Gamma_{m, 1}$: If strong effect relative to $b_{m, 1}$ (fixed), then this will generate noise on control samples, leading to false positives. 

$\Gamma_{m, 2}$: If strong effect relative to $b_{m, 2}$ (fixed), then this will generate noise on case samples, leading to false negatives.

Our certainty of $\Gamma$ to effect case or control samples depends on "the percentage of row space of $S$ explained by $G$". We appropx that by looking at $cor(G_r, S_2), r = 1, 2$. We probably can fix this value for now.



### Knob Speculation, within one experiment

| $\Gamma_{m, 1}$ 	| $\Gamma_{m, 2}$ 	| $cor(G_r, S_2)$ 	| **DE**           	| **Scree plot**     	|
|-----------------	|-----------------	|------------------	|------------------	|--------------------	|
| strong          	| weak            	| strong           	| more FPs         	| more even PCA      	|
| weak            	| strong          	| strong           	| more FNs         	| more even PCA      	|
| weak            	| weak            	| strong           	| neutral          	| more dominated PCA 	|
| strong          	| strong          	| strong           	| more FPs and FNs 	| more even PCA      	|


# Simulation with multiple experiments

$$
  X = \left[\begin{array}{ c | c }
    B_1 & B_2 \\
  \end{array}\right]
  \left[\begin{array}{ c | c }
    S_1 & 0 \\
    \hline
    0 & S_2
  \end{array}\right] +
  \left[\begin{array}{ c | c }
    \Gamma_1 & \Gamma_2 \\
  \end{array}\right]
  \left[\begin{array}{ c | c }
    G_1 & 0 \\
    \hline
    0 & G_2
  \end{array}\right] + U
$$
$$S_1 = 
\begin{bmatrix}
1 & ... & 1 & ... &  1\\
0 & ... & 1 & ... & 1
\end{bmatrix}$$

$$S_2 = 
\begin{bmatrix}
1 & ... & 1 & ... &  1\\
0 & ... & 1 & ... & 1
\end{bmatrix}$$

where study 2's effect relative to study 1 is modeled as $B_{21} \sim N(\alpha, 1)$
where $B_1$ and $B_2$ have the same set of DE genes with equal effects. 

We vary SV effects $\gamma$ via $\Gamma_{1:2, m, 1} = 0, m=1:300$,  $\Gamma_{1:2, m, 1} \sim N(\gamma, 1), m=301:1000$. 


```{r, echo=F, message=F, warning=F, results='hide'}
set.seed(2023)
m = 1000 #number of genes (tests)
n = 20 #number of samples
r = 2 #number of latent variables per studies
s = 2 #number of studies

DE_effect = 6
simulation = expand.grid(alpha = c(.75, 1.5, seq(3, 12, 3)), #effect of the second experiment relative to the first. 
                         gamma = c(.75, 1.5, seq(3, 12, 3))) #latent variable effect
simulation$n.sv = NA
simulation$num_PC_signif = NA
simulation$var_explained_PC1 = NA
simulation$sum_of_span_residue1_vs_truth = NA
simulation$sum_of_span_residue2_vs_truth = NA

forced_analysis = expand.grid(alpha = c(.75, 1.5, seq(3, 12, 3)), 
                              gamma = c(.75, 1.5, seq(3, 12, 3)),
                              forced_SV = 1:6) 
forced_analysis$sum_of_span_residue_joint_forced_vs_truth = NA

X_all = data.frame()

for(i in 1:nrow(simulation)) {
  alpha = simulation$alpha[i]
  gamma = simulation$gamma[i]
  
  #generate data
  U = matrix(rnorm(n * s * m, mean = 0, sd = 1), nrow = m)
  
  
  b11 = rnorm(m, mean = 0, sd = 1) #baseline effect in study 1
  b12 = mapply(mu = c(rep(DE_effect, 300), rep(0, m - 300)),
               sigma = c(rep(1, 300), rep(0, m - 300)),
               function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  B1 = cbind(b11, b12)
  b21 = rnorm(m, mean = alpha, sd = 1) #baseline effect in study 2
  b22 = mapply(mu = c(rep(DE_effect, 300), rep(0, m - 300)),
               sigma = c(rep(1, 300), rep(0, m - 300)),
               function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  B2 = cbind(b21, b22)
  B = cbind(B1, B2)
  
  S1 = matrix(c(rep(1, n), rep(0, n), rep(0, n/2), rep(1, n/2), rep(0, n)), byrow = T, ncol = s*n)
  S2 = matrix(c(rep(0, n), rep(1, n), rep(0, n), rep(0, n/2), rep(1, n/2)), byrow = T, ncol = s*n)
  S = rbind(S1, S2)
  
  
  gamma11 = mapply(mu = c(rep(0, 300), rep(gamma, m - 300)),
                   sigma = c(rep(0, 300), rep(1, m - 300)),
                   function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  gamma12 = mapply(mu = c(rep(-gamma, 300), rep(0, m - 300)),
                   sigma = c(rep(1, 300), rep(0, m - 300)),
                   function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  Gamma1 = cbind(gamma11, gamma12)
  gamma21 = mapply(mu = c(rep(0, 300), rep(gamma, m - 300)),
                   sigma = c(rep(0, 300), rep(1, m - 300)),
                   function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  gamma22 = mapply(mu = c(rep(-gamma, 300), rep(0, m - 300)),
                   sigma = c(rep(1, 300), rep(0, m - 300)),
                   function(mu, sigma) rnorm(1, mean = mu, sd = sigma))
  Gamma2 = cbind(gamma21, gamma22)
  Gamma = cbind(Gamma1, Gamma2)
  
  G_zeros = matrix(rep(0, n*2), byrow = T, ncol = n)
  G1 = cbind(mapply(p = c(rep(.4, 10), rep(.6, 10)),
                    function(p) rbinom(2, 1, p)), G_zeros)
  G2 = cbind(G_zeros, mapply(p = c(rep(.4, 10), rep(.6, 10)),
                             function(p) rbinom(2, 1, p)))
  G = rbind(G1, G2)
  
  X = B %*% S + Gamma %*% G + U
  Xnull = B %*% S + U
  X1 = B1 %*% S1 + Gamma1 %*% G1 + U
  X1 = X1[, 1:n]
  X2 = B2 %*% S2 + Gamma2 %*% G2 + U
  X2 = X2[, (n+1):(2*n)]
  
  X_long = data.frame(X) %>% 
            mutate(gene = 1:nrow(X)) %>%
            pivot_longer(!gene, names_to = "sample", values_to = "expression")
  X_long$sample = as.numeric(gsub("X", "", as.character(X_long$sample)))
  X_long$alpha = alpha
  X_long$gamma = gamma
  X_all = rbind(X_all, X_long)

  #Normalize
  # X[X < 0] = 0
  # dds <- DESeqDataSetFromMatrix(countData = X, colData = t(S_merged), design = ~t(S_merged)[, 2])
  # dds <- estimateSizeFactors(dds)
  # geneCountsNormalized <- counts(dds, normalized = TRUE)
  
  #Inference
  #PCA
  pca = prcomp(t(X))
  variance = pca$sdev^2 / sum(pca$sdev^2)
  simulation$num_PC_signif[i] = length(which(variance > 1/ncol(X)))
  simulation$var_explained_PC1[i] = variance[1]
  #SV
  #single study
  x1_n.sv = num.sv(X1, t(S1[, 1:n]), method = "be")
  x1_svobj = sva(X1, t(S1[, 1:n]), t(S1[, 1:n])[, 1], n.sv = x1_n.sv)
  x2_n.sv = num.sv(X2, t(S2[, (n+1):(2*n)]), method = "be")
  x2_svobj = sva(X2, t(S2[, (n+1):(2*n)]), t(S2[, (n+1):(2*n)])[, 1], n.sv = x2_n.sv)
  #Individual study spanning residue
  simulation$sum_of_span_residue1_vs_truth[i] = sum_of_span_residue(x1_svobj$sv, t(G1)[1:n ,])
  simulation$sum_of_span_residue2_vs_truth[i] = sum_of_span_residue(x2_svobj$sv, t(G2)[(n+1):(2*n) ,])
  #joint study
  simulation$n.sv[i] = num.sv(X, t(S), method = "be")
  #Joint study spanning residue
  idx = which(forced_analysis$alpha == simulation$alpha[i] & forced_analysis$gamma == simulation$gamma[i])
  for(j in idx) {
    cat("\n", forced_analysis$forced_SV[j], "\n")
    joint_svobj_forced = tryCatch({
        sva(X, t(S), t(S)[, c(1, 3)], n.sv = forced_analysis$forced_SV[j])
    }, error = function(e) {
        joint_svobj_forced = NA
    }, finally = {})
    if(!is.na(joint_svobj_forced)) {
      forced_analysis$sum_of_span_residue_joint_forced_vs_truth[j] = sum_of_span_residue(joint_svobj_forced$sv, t(G))
    }
  }

  # simulation$sum_of_span_residue1_vs_joint[i] = sum_of_span_residue(x1_svobj$sv, joint_svobj$sv[1:n, ])
  # simulation$sum_of_span_residue2_vs_joint[i] = sum_of_span_residue(x2_svobj$sv, joint_svobj$sv[(n+1):(2*n), ])
  # simulation$sum_of_span_residue1_vs_joint_forced[i] = sum_of_span_residue(x1_svobj$sv, joint_svobj_forced$sv[1:n ,])
  # simulation$sum_of_span_residue2_vs_joint_forced[i] = sum_of_span_residue(x2_svobj$sv, joint_svobj_forced$sv[(n+1):(2*n) ,])
  # 
  
  # myPairs(x1_svobj$sv, joint_svobj$sv[1:n ,], left = "Study 1", bottom = "Study 1 + 2")
  # myPairs(x2_svobj$sv, joint_svobj$sv[(n+1):(2*n) ,], left = "Study 2", bottom = "Study 1 + 2")
  # myPairs(x1_svobj$sv, joint_svobj_forced$sv[1:n ,], left = "Study 1", bottom = "Study 1 + 2 forced")
  # myPairs(x2_svobj$sv, joint_svobj_forced$sv[(n+1):(2*n) ,], left = "Study 2", bottom = "Study 1 + 2 forced")
  # 
}

```



```{r echo=F}
ggplot(simulation, aes(x = as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill =  n.sv)) + geom_tile() + 
  labs(x = "Noise:Signal of study effect = alpha/DE_effect",
       y = "Noise:Signal of latent variables = gamma/DE_effect") +
   ggtitle("Expected SV: 4")

ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = num_PC_signif)) + geom_tile() + 
  labs(x = "Noise:Signal of study effect = alpha/DE_effect",
       y = "Noise:Signal of latent variables = gamma/DE_effect",
       fill = "num_PC_signif") + ggtitle("Expected signif PC: 5")

ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = sum_of_span_residue1_vs_truth)) + geom_tile() + 
  labs(x = "Noise:Signal of study effect = alpha/DE_effect",
       y = "Noise:Signal of latent variables = gamma/DE_effect",
       fill = "Study 1\nSpan residule") +
  scale_fill_continuous(limits = c(0, .3))

ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = sum_of_span_residue2_vs_truth)) + geom_tile() + 
  labs(x = "Noise:Signal of study effect = alpha/DE_effect",
       y = "Noise:Signal of latent variables = gamma/DE_effect",
      fill = "Study 2\nSpan residule") +
  scale_fill_continuous(limits = c(0, .3))



# 
# ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = sum_of_span_residue1_vs_joint)) + geom_tile() + 
#   labs(x = "Noise:Signal of study effect = alpha/DE_effect",
#        y = "Noise:Signal of latent variables = gamma/DE_effect") +
#   scale_fill_continuous(limits = c(0, .1))
# 
# ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = sum_of_span_residue1_vs_joint_forced)) + geom_tile() + 
#   labs(x = "Noise:Signal of study effect = alpha/DE_effect",
#        y = "Noise:Signal of latent variables = gamma/DE_effect") +
#   scale_fill_continuous(limits = c(0, .1))
# 
# ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = sum_of_span_residue2_vs_joint)) + geom_tile() + 
#   labs(x = "Noise:Signal of study effect = alpha/DE_effect",
#        y = "Noise:Signal of latent variables = gamma/DE_effect") +
#   scale_fill_continuous(limits = c(0, .1))
# 
# ggplot(simulation, aes(x =  as.factor(alpha/DE_effect), y = as.factor(gamma/DE_effect), fill = sum_of_span_residue2_vs_joint_forced)) + geom_tile() + 
#   labs(x = "Noise:Signal of study effect = alpha/DE_effect",
#        y = "Noise:Signal of latent variables = gamma/DE_effect") +
#   scale_fill_continuous(limits = c(0, .1))
   

```


```{r echo=F, fig.width=14, fig.height=18}
ggplot(forced_analysis, aes(forced_SV, sum_of_span_residue_joint_forced_vs_truth)) + geom_point() + geom_line() + facet_grid(alpha/DE_effect ~ gamma/DE_effect)


# X = X_all %>% filter(alpha == 6 & gamma == 6) %>% select(-c("alpha", "gamma")) %>% pivot_wider(names_from = sample, values_from = expression) %>% select(-gene) %>% as.matrix()
# pheatmap(X)
# 
# ggplot(X, aes(x = sample, y = gene, fill = expression)) +  geom_tile()
# 
# ggplot(X_all, aes(x = sample, y = gene, fill = expression)) +  geom_tile() + facet_grid(alpha/DE_effect ~ gamma/DE_effect)
```