---
title: "compare_by_experiment"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(DESeq2)
require(data.table)
library(pheatmap)
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"))
allGenes = rownames(count_per_gen)
```

## DESeq by experiment
divide the dataset into 6 subsets
run DESeq on each subset and save:
1. dds results in .RData file
2. normalized counts in .csv files

```{r}
library(rmarkdown)
# divide the dataset into 6 subsets

sb1 = list()
sb1$name = "YPD_to_Low_N.0"
sb1$idx = which(sumExp$experiment_type== "YPD_to_Low_N" & sumExp$time_from_experiment_start_min == "0")

sb2 = list()
sb2$name = "YPD_to_Low_N.60"
sb2$idx = which(sumExp$experiment_type== "YPD_to_Low_N" & sumExp$time_from_experiment_start_min == "60")

sb3 = list()
sb3$name = "YPD_to_Low_N.960"
sb3$idx = which(sumExp$experiment_type== "YPD_to_Low_N" & sumExp$time_from_experiment_start_min == "960")

sb4 = list()
sb4$name = "CellCycle"
sb4$idx = which(sumExp$experiment_type== "CellCycle")

sb5 = list()
sb5$name = "Growth_curve"
sb5$idx = which(sumExp$experiment_type== "Growth_curve")

sb6 = list()
sb6$name = "SC_to_No_Pi"
sb6$idx = which(sumExp$experiment_type== "SC_to_No_Pi")

sb7 = list()
sb7$name = "all_experiments"
sb7$idx = 1:dim(sumExp)[1]


sbs = list(sb1, sb2, sb3, sb4, sb5, sb6, sb7)

# for (i in 1:length(sbs)) {
#   render("~/R_DESeq_local/compare_by_experiment_report_template.Rmd",
#        output_file=paste0("~/R_DESeq_local/compare_by_experiment_results/", sbs[[i]]$name, ".html"),
#        params= sbs[[i]])
# }

```


# load dds per experiment and analyze

```{r}
ddsList = list()
genes.sp = matrix(0, length(allGenes), length(sbs)-1)
genes.cis = matrix(0, length(allGenes), length(sbs)-1)
genes.trans = matrix(0, length(allGenes), length(sbs)-1)
currAlpha = 0.05
currLfcTh = 1

for (i in 1:(length(sbs)-1)) {
  currName = sbs[[i]]$name
  currList = list()
  currList$name = currName
  load(paste0(homeDir, "R_DESeq/DESeq_out/dds_gb_", currName, ".RData"))
  currList$dds = dds
  load(paste0(homeDir, "R_DESeq/DESeq_out/dds_trans_", currName, ".RData"))
  currList$dds2 = dds2
  ddsList[[i]] = currList
  
  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)
  
  write.csv(as.data.frame(res_sp), file = paste0(homeDir, "R_DESeq/DESeq_out/res_sp", currName, ".csv"))
  write.csv(as.data.frame(res_cis), file = paste0(homeDir, "R_DESeq/DESeq_out/res_cis", currName, ".csv"))
  write.csv(as.data.frame(res_trans), file = paste0(homeDir, "R_DESeq/DESeq_out/res_trans", currName, ".csv"))
  
  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)
  
  genes.sp[, i] = is.element(allGenes, rownames(dds)[sp.Genes])
  genes.cis[, i] = is.element(allGenes, rownames(dds)[cis.Genes])
  genes.trans[, i] = is.element(allGenes, rownames(dds2)[trans.Genes])
}

```

save significant genes to csv 
```{r}
experimentNames = vector()
for (i in 1:(length(sbs)-1)) {
  experimentNames[i] = sbs[[i]]$name
}

rownames(genes.sp) = allGenes
rownames(genes.cis) = allGenes
rownames(genes.trans) = allGenes
colnames(genes.sp) = experimentNames
colnames(genes.cis) = experimentNames
colnames(genes.trans) = experimentNames
write.csv(genes.sp, file=paste0(homeDir, "/R_DESeq/DESeq_out/genes.sp.csv"))
write.csv(genes.cis, file=paste0(homeDir, "/R_DESeq/DESeq_out/genes.cis.csv"))
write.csv(genes.trans, file=paste0(homeDir, "/R_DESeq/DESeq_out/genes.trans.csv"))
```

how many genes are changing in how many experiments?
```{r}
par(mfrow = c(1,3))
xl = c(.5, 6)
yl = c(0, 1000)
hist(rowSums(genes.sp), xlim = xl, ylim = yl)
hist(rowSums(genes.cis), xlim = xl, ylim = yl)
hist(rowSums(genes.trans), xlim = xl, ylim = yl)

```
## how many genes change per experiment?
```{r}
perExp.sp = colSums(genes.sp)
perExp.cis = colSums(genes.cis)
perExp.trans = colSums(genes.trans)

d = data.frame(experimentNames = experimentNames, sp = perExp.sp, cis= perExp.cis, trans = perExp.trans)
d2 = melt(data = d, id.vars = "experimentNames")

ggplot(data = d2, aes(x=experimentNames, y = value, fill=variable )) + geom_bar(stat="identity", position=position_dodge()) + theme(axis.text.x = element_text(angle = 90))
```


trans diverging genes
```{r}
dds = ddsList[[1]]$dds
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)

x = matrix(0, length(allGenes), 2)
x[is.element(allGenes, rownames(dds)), 1] = res_sp$log2FoldChange
x[is.element(allGenes, rownames(dds)), 2] = res_cis$log2FoldChange
plot(x[,1], x[,2])

always.trans = which(rowSums(genes.trans)==6)
mostly.trans = which(rowSums(genes.trans) >= 3)
points(x[always.trans, 1], x[always.trans, 2], col="red")
points(x[mostly.trans, 1], x[mostly.trans, 2], col="green")

mostly.cis = which(rowSums(genes.cis) >= 3)
points(x[mostly.cis, 1], x[mostly.cis, 2], col="purple")
```

combine normalized counts from all experiments into one matrix
```{r}
Nsamples = 2166
normCounts = data.frame(matrix(NA, length(allGenes), Nsamples))
justCounts = data.frame(matrix(NA, length(allGenes), Nsamples))
th = 1
colnames(normCounts) = colnames(count_per_gen)
colnames(justCounts) = colnames(count_per_gen)
for (i in 1:(length(sbs)-1)) {
  ntd = normTransform(ddsList[[i]]$dds)
  currCounts = assay(ntd)
  currCounts2 = matrix(NA, length(allGenes), dim(currCounts)[2])
  currCounts2[is.element(allGenes, rownames(currCounts)), ] = currCounts
  normCounts[, is.element(colnames(normCounts), colnames(currCounts))] = currCounts2
  
  currCounts = counts(ddsList[[i]]$dds, normalized=FALSE)
  currCounts2 = matrix(NA, length(allGenes), dim(currCounts)[2])
  currCounts2[is.element(allGenes, rownames(currCounts)), ] = currCounts
  justCounts[, is.element(colnames(normCounts), colnames(currCounts))] = currCounts2
}

```

## comparable conditions
using log2 normalized counts

```{r}
comparable_samples = read_xlsx(paste0(homeDir, "Screen/read_organize_data/comparable_samples.xlsx"))
cis_switched_replicates = read_xlsx(paste0(homeDir, "gene_gene_corr/cistrans_cross_replicates/cis_switched_replicates.xlsx"), 
                                    col_names=FALSE)

currNames = colnames(normCounts)

exp_cer = normCounts[, which(is.element(currNames, comparable_samples$cer))]
exp_cer = exp_cer[, match(comparable_samples$cer, colnames(exp_cer))]

exp_par = normCounts[, which(is.element(currNames, comparable_samples$par))]
exp_par = exp_par[, match(comparable_samples$par, colnames(exp_par))]

nameshyc = gsub("hyb", "hyc", comparable_samples$hyb)
nameshyp = gsub("hyb", "hyp", comparable_samples$hyb)

exp_hyc = normCounts[, which(is.element(currNames, nameshyc))]
exp_hyc = exp_hyc[, match(nameshyc, colnames(exp_hyc))]

exp_hyp = normCounts[, which(is.element(currNames, nameshyp))]
exp_hyp = exp_hyp[, match(nameshyp, colnames(exp_hyp))]

nameshyc.switched = cis_switched_replicates$...1
nameshyp.switched = gsub("hyc", "hyp", cis_switched_replicates$...1)

exp_hyc_replicates = normCounts[, which(is.element(currNames, nameshyc.switched))]
exp_hyc_replicates = exp_hyc_replicates[, match(nameshyc.switched, colnames(exp_hyc_replicates))]

exp_hyp_replicates = normCounts[, which(is.element(currNames, nameshyp.switched))]
exp_hyp_replicates = exp_hyp_replicates[, match(nameshyp.switched, colnames(exp_hyp_replicates))]

sp = exp_cer - exp_par
cis = exp_hyc - exp_hyp
cis_replicates = exp_hyc_replicates - exp_hyp_replicates
trans_sameRep = sp - cis
trans = sp - cis_replicates
nC = dim(sp)[2]
corr_cross_rep = vector(mode="numeric", nC)
corr_same_rep = vector(mode="numeric", nC)
for (i in 1:nC) {
  corr_cross_rep[i] = cor(cis[, i], trans[, i], use="complete.obs")
  corr_same_rep[i] = cor(cis[, i], trans_sameRep[, i], use="complete.obs")
}
par(mfrow=c(1,2),mar=c(2,2,1,1))
hist(corr_same_rep)
hist(corr_cross_rep)
```

check if OK: visualize expression levels of TDH3
```{r}
currGene = "YGR192C" # TDH3
# currGene = "YOL151W"
currGene = which(rownames(count_per_gen) == currGene)
d = data.frame(cer = as.vector(t(exp_cer[currGene, ])), par = as.vector(t(exp_par[currGene, ])), 
               hyc = as.vector(t(exp_hyc[currGene, ])), hyp = as.vector(t(exp_hyp[currGene, ])))
pheatmap(d, cluster_cols = FALSE)

```


## how many genes change per condition?
```{r}
N = 1
sp.Genes = which(rowSums(genes.sp) >= N)
cis.Genes = which(rowSums(genes.cis) >= N)
trans.Genes = which(rowSums(genes.trans) >= N)
length(sp.Genes)
length(cis.Genes)
length(trans.Genes)
currLfcTh = 1


genes.per.cond = data.frame(between_species = apply(sp, 2, function(x) sum(abs(x[sp.Genes]) > currLfcTh, na.rm=TRUE)),
                            cis = apply(cis, 2, function(x) sum(abs(x[cis.Genes]) > currLfcTh, na.rm=TRUE)),
                            trans = apply(trans, 2, function(x) sum(abs(x[trans.Genes]) > currLfcTh, na.rm=TRUE)))


genes.per.cond.1 = melt(data = genes.per.cond)

genes.per.cond.1 %>%
  ggplot(aes(x = value, group = variable, fill = variable)) + geom_histogram(alpha=.5, position="identity", binwidth=50) +
  xlab("# differentially-expressed genes") + ylab("# conditions") + theme_bw()
# genes.per.cond.1 %>%
#   ggplot(aes(x = value, group = variable, fill = variable)) + geom_density(alpha=.5)

```

## in how many conditions a gene changes?
```{r}

conds.per.gene = data.frame(between_species = apply(sp, 1, function(x) sum(abs(x) > currLfcTh, na.rm=TRUE)),
                            cis = apply(cis, 1, function(x) sum(abs(x) > currLfcTh, na.rm=TRUE)),
                            trans = apply(trans, 1, function(x) sum(abs(x) > currLfcTh, na.rm=TRUE)))
conds.per.gene.2 = matrix(0, length(allGenes), 3)
conds.per.gene.2[sp.Genes, 1] = conds.per.gene[sp.Genes, 1]
conds.per.gene.2[cis.Genes, 2] = conds.per.gene[cis.Genes, 2]
conds.per.gene.2[trans.Genes, 3] = conds.per.gene[trans.Genes, 3]
conds.per.gene = data.frame(between_species = conds.per.gene.2[, 1], 
                            cis = conds.per.gene.2[, 2],
                            trans = conds.per.gene.2[,3])

conds.per.gene.1 = melt(data = conds.per.gene)
conds.per.gene.1 %>% filter(value > 0) %>%
  # ggplot(aes(x = value, group = variable, fill = variable)) + geom_histogram(alpha=.5, position="identity") +
  # xlab("# conditions") + ylab("# genes") 
  ggplot(aes(x = value, group = variable, fill = variable)) + geom_density(alpha=.5, position="identity") +
  xlab("# conditions") + ylab("# genes") 

```

```{r}


```


