---
params:
  name: name
  idx: idx
title: "`r params$name`"
output: html_document
---
## load data
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(DESeq2)
library("dplyr")
library("ggplot2")
library(ashr)
library(vsn)
library(readxl)
library(knitr)

homeDir = "//data.wexac.weizmann.ac.il/barkailab/gatk/"
count_per_gen = readRDS(paste0(homeDir, "R_DESeq/count_per_gen.rds"))
sumExp = readRDS(paste0(homeDir, "R_DESeq/sumExp.rds"))
```

current subset of the data: 
```{r}
currExp = params$name
idx = params$idx
print(currExp)
sumExp1 = sumExp[idx, ]
sumExp1 = droplevels(sumExp1)
count_per_gen1 = count_per_gen[, idx]
```

## compare gene expression between the four genotypes: b/w species & cis 
design = ~gb
```{r}
dds = DESeqDataSetFromMatrix(count_per_gen1, sumExp1, design = ~ gb)
```

minminal count = 5 reads. How many genes are included?
```{r}
c = counts(dds)
Ngenes = dim(c)[1]
sum_per_gb = matrix(nrow=Ngenes, ncol=4)
gb = levels(colData(dds)$gb)
for (i in 1:length(gb)) {
  sum_per_gb[, i] = rowSums(c[ ,which(colData(dds)$gb == gb[i])])
}

minCounts = 5
good_genes_idx = apply(sum_per_gb, 1, function(x) all(x>minCounts))  
Ngenes = sum(good_genes_idx)
dds = dds[good_genes_idx, ]
print(dim(dds))
```

DESeq
```{r}
dds = DESeq(dds)
```

write the normalized counts to a csv file
```{r}
ntd = normTransform(dds)
write.csv(assay(ntd), file=paste0(homeDir, "R_DESeq/DESeq_out/log2NormCounts_", currExp, ".csv"))
```

## trans effect
The trans effect is ratio of ratios (trans = (cer/par)/(hyc/hyp)). It is therfore the interaction term F:sp. 
We use Likelihood Ratio Test (LRT) and a reduced model that doesn't involve the interaction. See: https://support.bioconductor.org/p/61509/
```{r}
dds2 = DESeqDataSetFromMatrix(count_per_gen1, sumExp1, design = ~F + sp + F:sp)
dds2 = dds2[good_genes_idx, ]
dds2 = DESeq(dds2, test="LRT", reduced = ~F+sp)
```
save files 
```{r}
save(dds, file = paste0(homeDir, "R_DESeq/DESeq_out/dds_gb_", currExp, ".RData"))
save(dds2, file = paste0(homeDir, "R_DESeq/DESeq_out/dds_trans_", currExp, ".RData"))
```
## Results + log fold change shrink
```{r}
res_sp = lfcShrink(dds, contrast = c("gb", "cer", "par"), type = "ashr", quiet=TRUE)
res_cis = lfcShrink(dds, contrast = c("gb", "hyc", "hyp"), type = "ashr", quiet=TRUE)
res_trans = lfcShrink(dds2, coef="F1.sppar", type="ashr", quiet=TRUE)

currAlpha = 0.05
currLfcTh = 1

sp.Genes = which(res_sp$padj < currAlpha & abs(res_sp$log2FoldChange) > currLfcTh)
cis.Genes = which(res_cis$padj < currAlpha & abs(res_cis$log2FoldChange) > currLfcTh)
trans.Genes = which(res_trans$padj < currAlpha & abs(res_trans$log2FoldChange) > currLfcTh)
```

```{r}
printResults = function(res) {
  pv = which(res$padj < currAlpha)
  LFC.high = intersect(which(res$log2FoldChange > currLfcTh), pv)
  LFC.low = intersect(which(res$log2FoldChange < (-currLfcTh)), pv)
  writeLines(paste(paste0("out of ", Ngenes, " with nonzero total read count"), 
        paste0("adjusted p-value < ", currAlpha), 
        paste0("LFC > ", currLfcTh, " (cer) : ", length(LFC.high), ", ", 
               format(round(length(LFC.high)*100/Ngenes, 2), nsmall = 2), "%"),
        paste0("LFC < -", currLfcTh, " (par) : ", length(LFC.low), ", ", 
               format(round(length(LFC.low)*100/Ngenes, 2), nsmall = 2), "%"), sep="\n"))
}
```

## DESeq statistics
```{r}
par(mfrow = c(3,3), mar = c(2,2,1,1))
plotMA(res_sp)
plotMA(res_cis)
plotMA(res_trans)

plot(res_sp$baseMean, -log10(res_sp$padj), xlab = "mean counts", ylab = "-log10(p-value)", log="x")
plot(res_cis$baseMean, -log10(res_cis$padj), xlab = "mean counts", ylab = "-log10(p-value)", log="x")
plot(res_trans$baseMean, -log10(res_trans$padj), xlab = "mean counts", ylab = "-log10(p-value)", log="x")

plot(res_sp$log2FoldChange, -log10(res_sp$padj), xlab = "log2 FC", ylab = "-log10(p-value)")
abline(v = -currLfcTh); abline(v = currLfcTh)
plot(res_cis$log2FoldChange, -log10(res_cis$padj), xlab = "log2 FC", ylab = "-log10(p-value)")
abline(v = -currLfcTh); abline(v = currLfcTh)
plot(res_trans$log2FoldChange, -log10(res_trans$padj), xlab = "log2 FC", ylab = "-log10(p-value)")
abline(v = -currLfcTh); abline(v = currLfcTh)
```





## How many genes in each category?

b/w species:
```{r}
printResults(res_sp)
```

cis:
```{r}
printResults(res_cis)
```

trans:
```{r}
printResults(res_trans)
```

## fold change 0

b/w species:
```{r}
currLfcTh = 0
printResults(res_sp)
```

cis:
```{r}
printResults(res_cis)
```

trans:
```{r}
printResults(res_trans)
```



## b/w species vs. cis plot
``` {r}
x = res_sp$log2FoldChange
y = res_cis$log2FoldChange

plot(x, y)
points(x[cis.Genes], y[cis.Genes], col="red")
points(x[trans.Genes], y[trans.Genes], col="blue")
cis_plus_trans.Genes = intersect(cis.Genes, trans.Genes)
points(x[cis_plus_trans.Genes], y[cis_plus_trans.Genes], col="purple")
```

