## RNAseq differential analysis

```{r}
diffexpress_by_strain = NULL
gsea_diff_express=NULL

for(strain in c("1728", "1734", "1736", "1740", "1744", "1747", "1751")) {
  counts = as.data.frame(counts_df) %>% 
    dplyr::select('DGY1657-1','DGY1657-2',"DGY1657-3", colnames(counts_df)[str_detect(colnames(counts_df), strain)])
  rownames(counts) = counts_df$gene
  coldata = data.frame(type = c("wt", "wt", "wt", rep("cnv", sum(str_detect(colnames(counts_df), strain)))), 
                       sample = c('DGY1657-1','DGY1657-2',"DGY1657-3", colnames(counts_df)[str_detect(colnames(counts_df), strain)]), 
                       row.names = colnames(counts))
  coldata$type = factor(coldata$type, levels = c('wt', 'cnv'))
  dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = coldata,
                              design = ~ type)
  dds <- DESeq(dds)
  res <- results(dds, alpha=0.05)
  res <- lfcShrink(dds, coef=2, res=res, type="apeglm")
  diffexpress_by_strain = bind_rows(diffexpress_by_strain, 
                                    as_tibble(res) %>% 
                                      mutate(gene = rownames(res), strain = strain))
  geneList = res$log2FoldChange
  gene.df <- bitr(rownames(res), fromType = "ORF",
        toType = "ENTREZID",
        OrgDb = org.Sc.sgd.db::org.Sc.sgd.db)
  names(geneList) = gene.df$ENTREZID
  geneList = sort(geneList, decreasing = TRUE)
  ego <- gseGO(geneList     = geneList,
             OrgDb        = org.Sc.sgd.db::org.Sc.sgd.db,
             keyType = "ENTREZID",
              ont          = "BP",
             by="fgsea",
              minGSSize    = 50,
              maxGSSize    = 400,
              pvalueCutoff = 0.05,
              pAdjustMethod="fdr",
              verbose      = FALSE,
             seed = 1)
  ego = clusterProfiler::simplify(ego)
  gsea_diff_express = bind_rows(gsea_diff_express, as_tibble(ego@result) %>% mutate(sample = strain))
}
#write out so this doesn't have to run again
write_csv(gsea_diff_express, paste0(data_dir,"/gsea_diff_expres.csv"))
write_csv(diffexpress_by_strain, paste0(data_dir,"/diffexpress_by_strain.csv"))
```

```{r}
## these should probably be in sup ####
gsea_diff_express=read_csv(paste0(data_dir,"/gsea_diff_expres.csv"))

diffexpress_by_strain=read_csv(paste0(data_dir,"/diffexpress_by_strain.csv")) %>%
  mutate(strain = as.character(strain)) %>%
  left_join(strain_names)

write_csv(diffexpress_by_strain, paste0(fig_dir,"/SupplementaryTable6.csv"))
```

# Significantly differently expressed genes

Big heatmap of Log2FC for significantly different genes
This is any gene that was differently expressed in any strain
```{r}
# get significant genes
s_genes = diffexpress_by_strain %>% filter(padj < 0.05) %>%  
  filter(abs(log2FoldChange)>1.5) %>%
  select(gene) %>% distinct() %>% as_vector()

# subset diff express
diff_express_sig = diffexpress_by_strain %>%
  filter(gene %in% s_genes) %>%
  select(log2FoldChange, gene, strain) %>%
  pivot_wider(names_from = strain, values_from = log2FoldChange)

mat = as.matrix(diff_express_sig[,2:8])
rownames(mat) = diff_express_sig$gene

require(igraph)

cut.knn <- function(x, nn=100) {c <- x[order(x)][min(nn, length(x))]; x[x > c] <- NA; return(x)}

# Distance matrix & 100-NN graph adjacency matrix (inverse distance weight)
pbc.dist <- dist(mat)

pbc.cut.dist <- as.matrix(pbc.dist)
pbc.cut.dist[upper.tri(pbc.cut.dist, diag = T)] <- NA
pbc.cut.dist.cut <- apply(pbc.cut.dist, 1, cut.knn, nn=200)

pbc.cut.dist.cut <- 1 / pbc.cut.dist.cut
pbc.cut.dist.cut[is.na(pbc.cut.dist.cut)] <- 0
pbc.graph <- igraph::graph.adjacency(pbc.cut.dist.cut, mode = "undirected", weighted = TRUE, diag = TRUE)

# Leiden clustering Traag et al. 2019 From Louvain to Leiden: guaranteeing well-connected communities
pbc.cluster <- cluster_leiden(pbc.graph, resolution_parameter = 0.1)
pbc.gene.cluster <- data.frame(Gene=rownames(mat), Cluster=0)
for (i in 1:length(pbc.cluster)) {pbc.gene.cluster[pbc.gene.cluster$Gene %in% pbc.cluster[[i]], "Cluster"] <- i}

# Heirarchial clustering
gene.y.hclust <- hclust(pbc.dist, 'complete')
gene.y.order <- as.character(gene.y.hclust$labels[gene.y.hclust$order])
gene.y.order <- tibble(gene = gene.y.order, order = 1:length(gene.y.order))

rownames(pbc.gene.cluster) = pbc.gene.cluster$Gene
pbc.gene.cluster = subset(pbc.gene.cluster, select = c("Cluster") )
pbc.gene.cluster$Cluster = as.character(pbc.gene.cluster$Cluster)
t=rowAnnotation(df = pbc.gene.cluster, 
                col = list(Cluster=c("1"="#7fc97f",
                                     "2" ="#beaed4",
                                     "3"="#fdc086",
                                     "4"="#ffff99",
                                     "5"="#386cb0",
                                     "6"="#f0027f",
                                     "7"="#bf5b17",
                                     "8"="#666666")))

Heatmap(mat,
        name = "log2FC",
        show_column_dend = FALSE,
        show_row_dend = FALSE,
        show_row_names = FALSE,
        cluster_rows = gene.y.hclust) + t

# find enriched go terms in each cluster
require(GOstats)
require(org.Sc.sgd.db)

#for gene universe use expressed genes
expressed_genes=diffexpress_by_strain %>% 
  filter(!is.na(log2FoldChange)) %>%
  select(gene) %>%
  distinct()
expressed_genes =  bitr(expressed_genes$gene, fromType = "ORF",
                        toType = "ENTREZID",
                        OrgDb = org.Sc.sgd.db::org.Sc.sgd.db)

cluster_genes <- bitr(rownames(pbc.gene.cluster), fromType = "ORF",
                      toType = "ENTREZID",
                      OrgDb = org.Sc.sgd.db::org.Sc.sgd.db)


enriched.go.terms <- list()
for (i in 1:max(pbc.gene.cluster$Cluster)) {
  params <- new("GOHyperGParams",
                geneIds =cluster_genes[pbc.gene.cluster$Cluster == i, "ORF"],
                universeGeneIds = expressed_genes$ORF,
                annotation="org.Sc.sgd.db",
                ontology="BP",
                pvalueCutoff=0.0001,
                conditional=T, # do T, see if smaller list
                testDirection="over")
  enriched.go.terms[[i]] <- hyperGTest(params)
}
summary(enriched.go.terms[[5]])
```