MT_tranplant_effect-revised-final.Rmd

---
title: "Microbiome transplants - transplant effect"
author: "Laura Brettell"
date: "06/10/2023"
output: html_document
keep_md: yes
---

# About

An experiment was carried out to determine how mosquitoes respond to a microbiome transplant. Transplants were performed using microbiomes from donors belonging to the same, or different mosquito species, collected from different environments (lab vs field). The recipient was always Aedes aegypti GALV laboratory-reared mosquitoes, reared in axenic conditions. A conventionally reared treatment was also included.  

The manuscript is entitled "Aedes aegypti gut transcriptomes respond differently to microbiome transplants from field-caught or laboratory-reared mosquitoes."

This half of the analysis address whether there were conserved responses to a transplant amongst recipients.



# Input required packages

```{r, eval=TRUE, message=FALSE}

#if (!requireNamespace("BiocManager", quietly = TRUE))
#        install.packages("BiocManager")
#BiocManager::install("DESeq2")


#if (!requireNamespace("BiocManager", quietly=TRUE))
#  install.packages("BiocManager")
#BiocManager::install("ComplexHeatmap")

#install.packages("UpSetR")
#install.packages("pheatmap")

```

```{r, eval=TRUE, message=FALSE}

library(DESeq2)
library(ggplot2)
library(UpSetR)
library(dplyr) # for %>%
library(tibble) # for tibble::rownames_to_column(df, "geneID")
library(ComplexHeatmap) # for extracting lists from upsetR and alternative heatmap plotting
library(ggrepel) # when overlaying text on volcano plots this makes then not overlap
library(readr) # for read_tsv
library(cowplot) # for plotting side by side
library(circlize) # for setting up colours in complexheatmaps
library(scales) # for sorting GOplots
```



# Input data and filter

These data are in the form of a raw counts file generated using HiSat2 and an associated metadata file. DEseq2 normalises to account for differences in library size and to reduce the effect of genes with zero or near-zero counts across all the data, so no prior normalisation is required. 

While it is not necessary to pre-filter low count genes before running the DESeq2 functions, in doing so we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. 

Here we are doing one initial filtering step, where we will keep only genes/rows which show >=10 reads in every replicate of at least one treatment.

```{r}
data <- read.csv("raw_counts.csv", row.names = 1) #import the count data

#Subset and generate vectors of logical arguments corresponding to whether the gene (row) has greater than 10 counts in all replicates for each condition 
data.AegF.L_GB_p <- apply(data[, 1:3], MARGIN = 1, function(x) all(x >= 10))
data.AegL.L_GB_p <- apply(data[, 4:6], MARGIN = 1, function(x) all(x >= 10))
data.AetF.AegL_GB_p <- apply(data[, 7:9], MARGIN = 1, function(x) all(x >= 10))
data.AetL.AegL_GB_p <- apply(data[, 10:12], MARGIN = 1, function(x) all(x >= 10))
data.Ag.Aeg_GB_p <- apply(data[, 13:15], MARGIN = 1, function(x) all(x >= 10))
data.CONV_p <- apply(data[, 16:18], MARGIN = 1, function(x) all(x >= 10))
data.Cu.Aeg_GB.P <- apply(data[, 19:21], MARGIN = 1, function(x) all(x >= 10))

 
#Merge the different logical vectors together
test.logical <- as.data.frame(rbind(data.AegF.L_GB_p, data.AegL.L_GB_p, data.AetF.AegL_GB_p, data.AetL.AegL_GB_p, data.Ag.Aeg_GB_p, data.CONV_p, data.Cu.Aeg_GB.P))

# we want this transposed - 
# we want a dataframe with 7 columns, one column for each condition
test.logical.2 <- as.data.frame(t(test.logical))

#At this point, at least one column in the row must be TRUE, which will indicate that all replicates in that condition have a count greater than 10
#We then get another vector of logical arguments, checking for the above
test.logical.3 <- rowSums(test.logical.2) >=1

#Finally, filter the original dataset with this final set of logicals
data.filtered <- data[test.logical.3, ]
# View(data.filtered)
# reduced from 19804 50 9723 genes

```


```{r}

coldata <- read.csv("coldata_expanded_meta.csv", row.names = 1) #import sample data
coldata <- coldata[,c("Microbiome.swap", "Direction", "microbiome.origin", "donor.env")]
coldata$Microbiome.swap <- factor(coldata$Microbiome.swap)
coldata$Direction <- factor(coldata$Direction)
coldata$microbiome.origin <- factor(coldata$microbiome.origin)
coldata$donor.env <- factor(coldata$donor.env)
head(coldata)

```

# Differential Expression analysis

## preparing the DESeq2 object 

Now we have the required data, we need to put it into an object of the required format, a DESeqDataSet (dds)


```{r}
dds <- DESeqDataSetFromMatrix(countData = data.filtered,
                              colData = coldata,
                              design = ~ Direction)

dds
```

Now to prepare for the Differential expression analysis, we ensure the reference level is set appropriately ie is the baseline to which all others will be compared (in our case the conventional treatment). This is all still essentially metadata at this point, no comparisons have been made yet. 

We will look at the Differential Expression of each of the recipient treatments individually compared to the Conventional control

```{r}
dds$Direction <- factor(dds$Direction, levels = c("AegL.L", "AegF.L", "AetF.AegL", "AetL.AegL", "Ang.Aeg", "Cu.Aeg", "Conventional"))

```


## running the differential expression analysis


```{r}
DE.dds <- DESeq(dds)

```
Now have a quick look at the results

```{r}

res.all <- results(DE.dds)
res.all
# by default this just shows the last comparison that R performed (alphabetically). 

```


## Extracting the results of interest and filtering

Firstly, extract all the results of interest ie every pairwise comparison to Conventional and apply cutoffs to remove those results which either didn't show much of a change, or which had a low probability of being correctly called. For this I will retain only the geneIDs in each pairwise Differential expression result with an adjusted p value of < 0.05 and a log2 fold change of >= 1.5.


```{r}

# extract the results of the first pairwise comparison as a dataframe
res.all.AegL.L.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AegL.L", "Conventional")))
# 9723 geneIDs

# filter by all those geneIDs with a log2 fold change of > 1.5 or <- 1.5
res.all.AegL.L.vs.conv <- res.all.AegL.L.vs.conv[abs(res.all.AegL.L.vs.conv$log2FoldChange) >= 1.5, ]
# 1042 geneIDs

# then filter by padj < 0.05, remove NAs (introduced during the initial filtering), and pull out the geneIDs
gene.IDs.AegL.L.vs.conv <- rownames(na.omit(res.all.AegL.L.vs.conv[res.all.AegL.L.vs.conv$padj < 0.05, ]))
# now is 730 geneIDs with na's omitted 

summary(gene.IDs.AegL.L.vs.conv)
```
Then repeat this for the rest of the pairwise comparisons (not the most elegant way of doing this!)

```{r}
res.all.AegF.L.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AegF.L", "Conventional")))
res.all.AegF.L.vs.conv <- res.all.AegF.L.vs.conv[abs(res.all.AegF.L.vs.conv$log2FoldChange)  >= 1.5, ]
gene.IDs.AegF.L.vs.conv <- rownames(na.omit(res.all.AegF.L.vs.conv[res.all.AegF.L.vs.conv$padj < 0.05, ]))
# 588 genes

res.all.AetF.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AetF.AegL", "Conventional")))
res.all.AetF.vs.conv <- res.all.AetF.vs.conv[abs(res.all.AetF.vs.conv$log2FoldChange)  >= 1.5, ]
gene.IDs.AetF.L.vs.conv <- rownames(na.omit(res.all.AetF.vs.conv[res.all.AetF.vs.conv$padj < 0.05, ]))  
# 195 genes

res.all.AetL.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AetL.AegL", "Conventional")))
res.all.AetL.vs.conv <- res.all.AetL.vs.conv[abs(res.all.AetL.vs.conv$log2FoldChange)  >= 1.5, ]
gene.IDs.AetL.L.vs.conv <- rownames(na.omit(res.all.AetL.vs.conv[res.all.AetL.vs.conv$padj < 0.05, ])) 
# 912 genes

res.all.Ang.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "Ang.Aeg", "Conventional")))
res.all.Ang.vs.conv <- res.all.Ang.vs.conv[abs(res.all.Ang.vs.conv$log2FoldChange)  >= 1.5, ]
gene.IDs.Ang.vs.conv <- rownames(na.omit(res.all.Ang.vs.conv[res.all.Ang.vs.conv$padj < 0.05, ]))
# 605 genes

res.all.Cu.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "Cu.Aeg", "Conventional")))
res.all.Cu.vs.conv <- res.all.Cu.vs.conv[abs(res.all.Cu.vs.conv$log2FoldChange)  >= 1.5, ]
gene.IDs.Cu.vs.conv <- rownames(na.omit(res.all.Cu.vs.conv[res.all.Cu.vs.conv$padj < 0.05, ]))
# 600 genes


# and stitch the lists
all.comps.to.conv <- list(gene.IDs.AegL.L.vs.conv, gene.IDs.AegF.L.vs.conv, gene.IDs.AetF.L.vs.conv, gene.IDs.AetL.L.vs.conv, gene.IDs.Ang.vs.conv, gene.IDs.Cu.vs.conv)


names(all.comps.to.conv) <- c("gene.IDs.AegL.L.vs.conv", "gene.IDs.AegF.L.vs.conv", "gene.IDs.AetF.L.vs.conv", "gene.IDs.AetL.L.vs.conv", "gene.IDs.Ang.vs.conv", "gene.IDs.Cu.vs.conv")
```


Then compile all geneIDs of interest ie all those genes that pass the threshold in at least one pairwise comparison. These lists are then concatenated and duplicates removed to obtain a full geneID list of interest. 


```{r}
# first to create a vector containing all the genes to be used 

all.gene.IDs <- gene.IDs.AegL.L.vs.conv
all.gene.IDs <- append(all.gene.IDs, gene.IDs.AegF.L.vs.conv)
all.gene.IDs <- append(all.gene.IDs, gene.IDs.AetF.L.vs.conv)
all.gene.IDs <- append(all.gene.IDs, gene.IDs.AetL.L.vs.conv)
all.gene.IDs <- append(all.gene.IDs, gene.IDs.Ang.vs.conv)
all.gene.IDs <- append(all.gene.IDs, gene.IDs.Cu.vs.conv)

# this vector contains 3630 genes,some of which are duplicates (ie feature in multiple lists), so need removing. 

all.gene.IDs.dupl.removed <- all.gene.IDs[!duplicated(all.gene.IDs)]
summary(all.gene.IDs.dupl.removed)
# now there are 1680 unique genes for inclusion in the heatmap 

```


# Upset plot to visualise which DE genes are common to multpile comparisons

Now to just consider which genes were identified in each comparison and whether the same genes were identified in multiple comparisons


```{r first_upset_plot}
# duplicating so i can try to make the names of the sets nicer for the plots, but still keeping the original just in case
all.comps.to.conv.renamed <- all.comps.to.conv
names(all.comps.to.conv.renamed) # to check the order

# make names nicer
names(all.comps.to.conv.renamed) <- c("Ae. aegypti (lab)", "Ae. aegypti (field)", "Ae. taeniorhynchus (field)", "Ae. taeniorhynchus (lab)", "An. gambiae (lab)", "Cx. tarsalis (lab)")

upsetplot1 <- UpSetR::upset(fromList(all.comps.to.conv.renamed), 
                order.by = "freq", 
                nsets = 38,
                queries = list(list(query = intersects, params = list(
                  "Ae. aegypti (lab)", 
                  "Ae. aegypti (field)", 
                  "Ae. taeniorhynchus (field)", 
                  "Ae. taeniorhynchus (lab)", 
                  "An. gambiae (lab)", 
                  "Cx. tarsalis (lab)"), 
                  color = "lightseagreen", active = T)))

upsetplot1

```

There are lots of genes which are only on one, or a subset of comparisons, but there are 71 (in teal) which are common to all swaps. 


# investigating the 'swap effect' in more depth

Now to look at the 71 genes that were identified as being enriched in every swap relative to the conventional treatment. We don't know yet if they are affected in the same direction.


So, firstly to extract the list of 71 geneIDs to assign annotations         
        
```{r}
# first, create a matrix of counts from the upsetR plot

comb_mat_to_conv <- make_comb_mat(all.comps.to.conv) # this is from complexheatmap package
comb_mat_to_conv

```


```{r}
# then we want to extract the genes in our intersect of interest by calling them using binary codes
#Pulls out the set containing all items ie all 1's
shared_genes_to_conv <- extract_comb(comb_mat_to_conv, "111111")

```


# Are these genes showing the same direction of regulation in the different comaprisons? 

This requires extracting all the fold change results from each of the pairwise comparisons for the 71 geneIDs of interest, then collating and creating a heatmap


# extracting and collating results

```{r}

# extract the result
res.full.AegF.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AegF.L", "Conventional")))

# just take the log2foldchange column
foldchange.AegF <- res.full.AegF.vs.conv  %>% select(log2FoldChange)

# rename the column to the sample name (otherwise when I bind, all samples will be indistinguishable)
colnames(foldchange.AegF) <- c("AegF")

# give the geneIDs their own column, to be able to merge later
foldchange.AegF <- tibble::rownames_to_column(foldchange.AegF, "geneID")

# now for the other pairwise comps

res.full.AegL.L.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AegL.L", "Conventional")))
res.full.AetF.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AetF.AegL", "Conventional")))
res.full.AetL.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "AetL.AegL", "Conventional")))
res.full.Ang.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "Ang.Aeg", "Conventional")))
res.full.Cu.vs.conv <- as.data.frame(results(DE.dds, contrast = c("Direction", "Cu.Aeg", "Conventional")))

foldchange.AegL <- res.full.AegL.L.vs.conv  %>% select(log2FoldChange)
foldchange.AetF <- res.full.AetF.vs.conv  %>% select(log2FoldChange)
foldchange.AetL <- res.full.AetL.vs.conv  %>% select(log2FoldChange)
foldchange.Ang <- res.full.Ang.vs.conv  %>% select(log2FoldChange)
foldchange.Cu <- res.full.Cu.vs.conv  %>% select(log2FoldChange)

colnames(foldchange.AegL) <- c("AegL")
colnames(foldchange.AetF) <- c("AetF")
colnames(foldchange.AetL) <- c("AetL")
colnames(foldchange.Ang) <- c("Ang")
colnames(foldchange.Cu) <- c("Cu")

foldchange.AegL <- tibble::rownames_to_column(foldchange.AegL, "geneID")
foldchange.AetF <- tibble::rownames_to_column(foldchange.AetF, "geneID")
foldchange.AetL <- tibble::rownames_to_column(foldchange.AetL, "geneID")
foldchange.Ang <- tibble::rownames_to_column(foldchange.Ang, "geneID")
foldchange.Cu <- tibble::rownames_to_column(foldchange.Cu, "geneID")

bound <-merge(foldchange.AegF, foldchange.AegL, by="geneID")
bound <-merge(bound, foldchange.AetF, by="geneID")
bound <-merge(bound, foldchange.AetL, by="geneID")
bound <-merge(bound, foldchange.Ang, by="geneID")
bound <-merge(bound, foldchange.Cu, by="geneID")

head(bound)

# now to re-set the row names as geneIDs
rownames(bound) <- bound$geneID

# remove the redundant column 'geneID'
bound <- bound[c(2:7)]

# and subset the results to keep only the genes of interest
bound.shared.to.conv <- bound[shared_genes_to_conv, ]
```


# make a heatmap of up/down regulation (log2fold change) from conventional


```{r}

col_fun3 = colorRamp2(c(-20, -10.001, -10, -1.5, -1.499, 1.499, 1.5, 10, 10.001, 20), c( "#2166AC", "#2166AC", "#2166AC", "#D1E5F0","snow3", "snow3", "#FDDBC7", "#B2182B", "#B2182B", "#B2182B"))
col_fun3(seq(-10, 10))

heatmap_shared_to_conv_no_rownames <- Heatmap(bound.shared.to.conv, name = "log2foldchange", col = col_fun3, show_row_names=F)

heatmap_shared_to_conv_no_rownames

```

and write the log2fold change dataframe to csv for supp

```{r}
# write.csv(bound.shared.to.conv, "bound.shared.to.conv.csv")

```


All 71 geneIDs identified as deferentially expressed in all comparisons to conventional have consistent directions of change. So, make a separate version of this heatmap with rownames, then extract the up and down regulated genenames and assign Vectorbase product descriptions and kegg metabolic pathways.

```{r}

heatmap_shared_to_conv_with_rownames <- Heatmap(bound.shared.to.conv, name = "log2foldchange", col = col_fun3, show_row_names=T, row_names_gp = gpar(fontsize = 3))

heatmap_shared_to_conv_with_rownames

```


# Pull together all transcripts in Upsetplot1 and identify their direction of change

The heatmap of the 71 genes identified in all comparisons showed consistency in the direction of change. Does this pattern hold more generally?


```{r}

# I previously extracted the results from each pairwise comparison to conv
# Here I will take the log2foldchange column AND the padj column
foldchange.padj.AegF <- res.full.AegF.vs.conv  %>% select(log2FoldChange, padj)

# rename the column to the sample name 
colnames(foldchange.padj.AegF) <- c("AegF", "padjAegF")

# give the geneIDs their own column, to be able to merge later
foldchange.padj.AegF <- tibble::rownames_to_column(foldchange.padj.AegF, "geneID")

# if padj >0.05, replace the value in the fold change column (now AegF) with NA (later rows will remain for non-signif values - if the value was signif in another comparison, so these need replacing with NAs). Doing this in 3 parts is not efficient but works.
foldchange.padj.AegF <- foldchange.padj.AegF %>% mutate(AegF = ifelse(padjAegF >= 0.05, "NA", AegF))

```

now for the other pairwise comps...

```{r}
foldchange.padj.AetF <- res.full.AetF.vs.conv  %>% select(log2FoldChange, padj)
foldchange.padj.AetL <- res.full.AetL.vs.conv  %>% select(log2FoldChange, padj)
foldchange.padj.AegL <- res.full.AegL.L.vs.conv  %>% select(log2FoldChange, padj)
foldchange.padj.Ang <- res.full.Ang.vs.conv  %>% select(log2FoldChange, padj)
foldchange.padj.Cu <- res.full.Cu.vs.conv  %>% select(log2FoldChange, padj)

colnames(foldchange.padj.AetF) <- c("AetF", "padjAetF")
colnames(foldchange.padj.AetL) <- c("AetL", "padjAetL")
colnames(foldchange.padj.AegL) <- c("AegL", "padjAegL")
colnames(foldchange.padj.Ang) <- c("Ang", "padjAng")
colnames(foldchange.padj.Cu) <- c("Cu", "padjCu")

foldchange.padj.AetF <- tibble::rownames_to_column(foldchange.padj.AetF, "geneID")
foldchange.padj.AetL <- tibble::rownames_to_column(foldchange.padj.AetL, "geneID")
foldchange.padj.AegL <- tibble::rownames_to_column(foldchange.padj.AegL, "geneID")
foldchange.padj.Ang <- tibble::rownames_to_column(foldchange.padj.Ang, "geneID")
foldchange.padj.Cu <- tibble::rownames_to_column(foldchange.padj.Cu, "geneID")

foldchange.padj.AetF <- foldchange.padj.AetF %>% mutate(AetF = ifelse(padjAetF >= 0.05, "NA", AetF))
foldchange.padj.AetL <- foldchange.padj.AetL %>% mutate(AetL = ifelse(padjAetL >= 0.05, "NA", AetL))
foldchange.padj.AegL <- foldchange.padj.AegL %>% mutate(AegL = ifelse(padjAegL >= 0.05, "NA", AegL))
foldchange.padj.Ang <- foldchange.padj.Ang %>% mutate(Ang = ifelse(padjAng >= 0.05, "NA", Ang))
foldchange.padj.Cu <- foldchange.padj.Cu %>% mutate(Cu = ifelse(padjCu >= 0.05, "NA", Cu))


#then combine
bound.padj <-merge(foldchange.padj.AegF, foldchange.padj.AegL, by="geneID")
bound.padj <-merge(bound.padj, foldchange.padj.AetF, by="geneID")
bound.padj <-merge(bound.padj, foldchange.padj.AetL, by="geneID")
bound.padj <-merge(bound.padj, foldchange.padj.Ang, by="geneID")
bound.padj <-merge(bound.padj, foldchange.padj.Cu, by="geneID")

head(bound.padj)

# now to re-set the row names as geneIDs
rownames(bound.padj) <- bound.padj$geneID

# remove the redundant column 'geneID' and the padj columns
bound.padj <- bound.padj[c(2,4,6,8,10,12)]

# and subset the results to keep only the genes of interest
bound.padj <- bound.padj[all.gene.IDs.dupl.removed, ]

# write a csv of these results for supp
#write.csv(bound.padj, "all_comps_to_conv.csv")
```


```{r}

bound.padj.no.nas <- bound.padj

bound.padj.no.nas <- data.frame(apply(bound.padj.no.nas, 2, as.numeric), check.names=, row.names = rownames(bound.padj.no.nas))  # Convert all variable types to numeric, keeping row names as they are

bound.padj.no.nas[is.na(bound.padj.no.nas)] <- 0 # change all NAs to 0
```

```{r}

heatmap_1680_genes_to_conv <- Heatmap(bound.padj.no.nas, name = "log2fold change", col = col_fun3, show_row_names=F)

heatmap_1680_genes_to_conv
```


Make a version with row names (GeneIDs) to use for GO enrichment analysis

```{r}
heatmap_1680_genes_to_conv_rownames <- Heatmap(bound.padj.no.nas, name = "mat", col = col_fun3, row_names_gp = gpar(fontsize = 0.5))
heatmap_1680_genes_to_conv_rownames

```

from the pdf of this plot, then take the row names for all those genes upregulated, downregulated, or a mix
'fold_changes_compared_to_conventional.xlsx'
Then use those genelists as inputs to Vectorbase GO enrichment analysis



## generating and plotting GO enrichment data

Here, we assign GO terms to the genes that were consistently enhanced or suppressed in at least one transplant group across the dataset as a whole (ie the 'fold_changes_compared_to_conventional.xlsx') and used this as input to Vectorbase Gene Ontology enrichment analysis (https://vectorbase.org/vectorbase/app) tool (once for the enhanced and once for the suppressed lists). Now to import the results and plot


```{r}
# the lists of terms upregulated compared to conventional - ie when a swap is done
# biological process
go_enhanced_BP <- read_tsv("hiddenGoEnrichmentResult_up_to_conv_BP.tsv")

# cutoff those with a bonferoni adjusted p value of > 0.05 
go_enhanced_BP <- go_enhanced_BP[go_enhanced_BP$Bonferroni <= 0.05, ]
#make a new column for fold change for later merging with the enhanced list
go_enhanced_BP$foldchange2 <- go_enhanced_BP$`Fold enrichment`
# make a new column to specify which ontology category
go_enhanced_BP$ontology <- "Biological_process"



go_enhanced_MF <- read_tsv("hiddenGoEnrichmentResult_up_to_conv_MF.tsv")
go_enhanced_MF <- go_enhanced_MF[go_enhanced_MF$Bonferroni <= 0.05, ]
go_enhanced_MF$foldchange2 <- go_enhanced_MF$`Fold enrichment`
go_enhanced_MF$ontology <- "Molecular_function"


go_enhanced_CC <- read_tsv("hiddenGoEnrichmentResult_up_to_conv_CC.tsv")
go_enhanced_CC <- go_enhanced_CC[go_enhanced_CC$Bonferroni <= 0.05, ]
go_enhanced_CC$foldchange2 <- go_enhanced_CC$`Fold enrichment`
go_enhanced_CC$ontology <- "Cellular_component"


go_enhanced_all <- rbind(go_enhanced_BP, go_enhanced_MF, go_enhanced_CC)

```


```{r}
# the list of terms downregulated/suppressed compared to conventional - ie when a swap is done

go_suppressed_BP <- read_tsv("hiddenGoEnrichmentResult_down_to_conv_BP.tsv")
go_suppressed_BP <- go_suppressed_BP[go_suppressed_BP$Bonferroni <= 0.05, ]
go_suppressed_BP$foldchange2 <- go_suppressed_BP$`Fold enrichment`*(-1)
# this time I changed the fold enrichment to a negative to denote the suppression/downregulation
go_suppressed_BP$ontology <- "Biological_process"

go_suppressed_MF <- read_tsv("hiddenGoEnrichmentResult_down_to_conv_MF.tsv")
go_suppressed_MF <- go_suppressed_MF[go_suppressed_MF$Bonferroni <= 0.05, ]
go_suppressed_MF$foldchange2 <- go_suppressed_MF$`Fold enrichment`*(-1)
go_suppressed_MF$ontology <- "Molecular_function"


go_suppressed_CC <- read_tsv("hiddenGoEnrichmentResult_down_to_conv_CC.tsv")
go_suppressed_CC <- go_suppressed_CC[go_suppressed_CC$Bonferroni <= 0.05, ]
go_suppressed_CC$foldchange2 <- go_suppressed_CC$`Fold enrichment`*(-1)
go_suppressed_CC$ontology <- "Cellular_component"


go_suppressed_all <- rbind(go_suppressed_BP, go_suppressed_MF, go_suppressed_CC)

```

```{r}
# the list of terms suppressed compared to conventional - ie when a swap is done

go_all <- rbind(go_enhanced_all, go_suppressed_all)
go_all$neg_log10_bonf <- -log(go_all$Bonferroni)

# write.csv(go_all, "GOresults_mtransplant_effect_all.csv")

```

## plot

```{r}
# make an extra column to say whether there is a positive or negative change
go_all$direction_change <- go_all$foldchange2

go_all <- go_all %>% 
  mutate(direction_change = if_else(direction_change > 0, "pos", "neg"))

#plot
GOplot2 <- ggplot(data = go_all, aes(x = reorder(Name, neg_log10_bonf), y = foldchange2, fill = direction_change)) + 
        geom_col() + 
        scale_fill_manual(values = c("pos" = "#B2182B", "neg" = "#2166AC")) +
        scale_x_discrete(labels = label_wrap(55)) + 
       facet_grid(ontology~., switch = "y", scales = "free_y", space = "free") +

#Below flips the x- and y-axis
        coord_flip() + 
        theme_light() +
        theme(axis.title = element_text(size = rel(2)), axis.text = element_text(size = rel(0.9)), axis.title.y = element_blank(), axis.title.x = element_blank(), legend.position = "none") #+ 
      #  scale_y_continuous(labels = scales::percent, limits = c(0,1)) + 
      #  labs(fill = "-Log10(P-value)")

GOplot2

```

```{r}

# combine plots with 1 row and 2 columns
par(mfcol = c(1,2))

heatmap_1680_genes_to_conv

GOplot2

```