# Written by: Reilly Pidgeon

# Purpose: Script to merge exported metadata from individual studies, then to merge the combined metadata to relative abundance and gene families from individual studies

  # Info on curatedMetagenomicData can be found here: 
  # Bioconductor info: https://bioconductor.riken.jp/packages/3.14/data/experiment/manuals/curatedMetagenomicData/man/curatedMetagenomicData.pdf
  # All available studies: https://waldronlab.io/curatedMetagenomicData/articles/available-studies.html
  # General: https://waldronlab.io/curatedMetagenomicData/articles/curatedMetagenomicData.html

# The following packages should be downloaded and loaded using the library() command
library(data.table)
library(dplyr)
library(tidyverse)
library(purrr)
library(ggplot2)
library(RColorBrewer)
library(ggpmisc)

#======================================== METADATA =====================================================

# Set the working directory to the location where the folders with the following names are located: metadata, ra_boi, gf_goi
# Each of these folders contain parsed data (in csv format) from individual studies in curatedMetagenomicData (searched for bacteria of interest and genes of interest)

setwd("C:/Users/Reilly/Desktop/")

setwd("metadata/")

# List metadata files in the metadata folder
md_files <- list.files(pattern = "*.csv")
md_files_list = lapply(md_files, read.csv)

# There are 93 studies total - the individual studies in the folder are named Study#_md.csv

# The first column contains numbers that will mess up the merging of metadata from all sudies, so it will be removed for each csv file
for (i in 1:93) {
  print(paste("Removing first column from", i))
  md_files_list[[i]] <- md_files_list[[i]][-1]
  # Since some samples from different studies have the same identifiers, a unique number should be added to the sample_id
  print(paste("Adding index number", i, "to id"))
  md_files_list[[i]]$id <- paste(md_files_list[[i]]$id,"-",i, sep = "") #This is needed to distinguish samples with the same id in different studies
}

merged_df_md <- md_files_list[[1]] # Here, a new dataframe is created containing the metadata from the first study - all metadata from other studies will be appended to this df

# Loop through the remaining studies and append to the dataframe above
for (i in 2:93) {
  print(i)
  temp <- md_files_list[[i]]
  merged_df_md <- dplyr::full_join(merged_df_md, temp) # The full_join function is needed to append metadata which may contain different columns (as not all metadata is uniform across studies - lots of NA values)
}

# Change all periods and underscores to dashes in "id" column to be able to easily merge tables
merged_df_md$id <- str_replace_all(merged_df_md$id,"\\.","-")
merged_df_md$id <- str_replace_all(merged_df_md$id,"_","-")

# Write a csv containing merged metadata from all studies
write.csv(merged_df_md, file = "merged_metadata.csv")


#==================================== RELATIVE ABUNDANCE =====================================

setwd("C:/Users/Reilly/Desktop/")

setwd("ra_boi/")

# List metadata files in the metadata folder
ra_files <- list.files(pattern = "*.csv")
ra_files_list = lapply(ra_files, read.csv)

# The first column contains numbers that will mess up the merging of metadata from all sudies, so it will be removed for each csv file
for (i in 1:93) {
  print(paste("Removing first column from", i))
  ra_files_list[[i]] <- ra_files_list[[i]][-1]
  print(paste("Adding index number", i, "to column names"))
  names(ra_files_list[[i]]) <- paste(names(ra_files_list[[i]]),"-",i, sep = "")
  # Replace "taxa-i" with "taxa"
  names(ra_files_list[[i]])[1] <- "taxa"
}

merged_df_ra <- ra_files_list[[1]] # Here, a new dataframe is created containing the relative abundance from the first study - all ra from other studies will be appended to this df

# Loop through the remaining studies and append to the dataframe above
for (i in 2:93) {
  print(i)
  if(class(ra_files_list[[i]]$taxa) == "logical") {
    print(paste("Skipped iteration", i))
    next
  } else {
  temp <- ra_files_list[[i]]
  merged_df_ra <- dplyr::full_join(merged_df_ra, temp, by = "taxa")
  }
}
# Transpose merged ra matrix
merged_df_ra_t <- t(merged_df_ra)
# Rename the columns to the first row
colnames(merged_df_ra_t) <- merged_df_ra_t[1,]
# Remove the first row
merged_df_ra_t <- merged_df_ra_t[-1,]
# Convert to dataframe
merged_df_ra_t <- as.data.frame(merged_df_ra_t)
# Make the first column the "id" to match to metadata
merged_df_ra_t <- tibble::rownames_to_column(merged_df_ra_t, "id")
# Convert NAs to 0
merged_df_ra_t[is.na(merged_df_ra_t)] <- 0

# Change all periods and underscores to dashes in "id" column to be able to easily merge tables
merged_df_ra_t$id <- str_replace_all(merged_df_ra_t$id,"\\.","-")
merged_df_ra_t$id <- str_replace_all(merged_df_ra_t$id,"_","-")

# Write a csv containing merged relative abundance from all studies
write.csv(merged_df_ra_t, file = "merged_relative_abundance_t.csv")

#========================================= GENE FAMILIES =====================================

setwd("C:/Users/Reilly/Desktop/")

setwd("gf_goi/")

# List metadata files in the metadata folder
gf_files <- list.files(pattern = "*.csv")
gf_files_list = lapply(gf_files, read.csv)

# The first column contains numbers that will mess up the merging of metadata from all sudies, so it will be removed for each csv file
for (i in 1:93) {
  print(paste("Removing first column from", i))
  gf_files_list[[i]] <- gf_files_list[[i]][-1]
  print(paste("Adding index number", i, "to column names"))
  names(gf_files_list[[i]]) <- paste(names(gf_files_list[[i]]),"-",i, sep = "")
  # Replace "gene_families-i" with "gene_families"
  names(gf_files_list[[i]])[1] <- "gene_families"
}

merged_df_gf <- gf_files_list[[1]] # Here, a new dataframe is created containing the gene families from the first study - all gf from other studies will be appended to this df

# Loop through the remaining studies and append to the dataframe above
for (i in 2:93) {
  print(i)
  if(class(gf_files_list[[i]]$gene_families) == "logical") {
    print(paste("Skipped iteration", i))
    next
  } else {
    temp <- gf_files_list[[i]]
    merged_df_gf <- dplyr::full_join(merged_df_gf, temp, by = "gene_families")
  }
}
# Transpose merged ra matrix
merged_df_gf_t <- t(merged_df_gf)
# Rename the columns to the first row
colnames(merged_df_gf_t) <- merged_df_gf_t[1,]
# Remove the first row
merged_df_gf_t <- merged_df_gf_t[-1,]
# Convert to dataframe
merged_df_gf_t <- as.data.frame(merged_df_gf_t)
# Make the first column the "id" to match to metadata
merged_df_gf_t <- tibble::rownames_to_column(merged_df_gf_t, "id")
# Convert NAs to 0
merged_df_gf_t[is.na(merged_df_gf_t)] <- 0

# Change all periods and underscores to dashes in "id" column to be able to easily merge tables
merged_df_gf_t$id <- str_replace_all(merged_df_gf_t$id,"\\.","-")
merged_df_gf_t$id <- str_replace_all(merged_df_gf_t$id,"_","-")

# Write a csv containing merged gene families from all studies
write.csv(merged_df_gf_t, file = "merged_gene_families_t.csv")


# ================ MERGE RELATIVE ABUNDANCE AND GENE FAMILIES WITH METADATA ==========================

setwd("C:/Users/Reilly/Desktop/")

# Merge data tables based on the second argument (y) in this case
ra_gf_with_metadata <- dplyr::left_join(merged_df_md, merged_df_ra_t, by = "id") 
ra_gf_with_metadata <- dplyr::left_join(ra_gf_with_metadata, merged_df_gf_t, by = "id")

write.csv(ra_gf_with_metadata, file = "asv_table_boi_ra_with_metadata.csv")
print("Relative abundance and gene families with metadata saved")


# =================================== GRAPHING =====================================

# Filter all studies to only keep stool samples (exclude skin, milk, vagina, oralcavity, nasalcavity)
cMD_stool <- ra_gf_with_metadata
cMD_stool <- filter(cMD_stool, body_site == "stool")
# Save a csv containing only stool samples
write.csv(cMD_stool, file = "cMD_stool.csv")

#------------------------- Check if metabolizers are present ----------------------------------------------------------
    
      
    # Create a metabolizer_present column that returns TRUE if any of the relative adundance values for E. aspara, E. bolteae, or E. citroniae are above 0
    cMD_stool$metabolizer_present <- ifelse(cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Lachnospiraceae|g__Lachnoclostridium|s__Clostridium_asparagiforme` > 0 |
                                              cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Lachnospiraceae|g__Lachnoclostridium|s__Clostridium_bolteae` > 0 |
                                              cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Clostridiaceae|g__Clostridium|s__Clostridium_bolteae_CAG_59` > 0 |
                                              cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Lachnospiraceae|g__Lachnoclostridium|s__Clostridium_citroniae` > 0 |
                                              cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Clostridiales_unclassified|g__Clostridiales_unclassified|s__Clostridiales_bacterium_1_7_47FAA` > 0, TRUE, FALSE)
    summary(cMD_stool$metabolizer_present)
    
    # Create a new dataframe with id, study_name, metabolizer_present (TRUE or FALSE), and the realtive abundance of each species
    metabolizer_summary <- data.frame(id = cMD_stool$id, study_name = cMD_stool$study_name, metabolizer_present = cMD_stool$metabolizer_present, 
                                      E_asparagiformis = cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Lachnospiraceae|g__Lachnoclostridium|s__Clostridium_asparagiforme`, 
                                      E_bolteae = cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Lachnospiraceae|g__Lachnoclostridium|s__Clostridium_bolteae`,
                                      E_bolteae_CAG59 = cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Clostridiaceae|g__Clostridium|s__Clostridium_bolteae_CAG_59`,
                                      E_citroniae = cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Lachnospiraceae|g__Lachnoclostridium|s__Clostridium_citroniae`,
                                      L_pacaense = cMD_stool$`k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Clostridiales_unclassified|g__Clostridiales_unclassified|s__Clostridiales_bacterium_1_7_47FAA`)
    
    # Convert NA values in the genes_present column to FALSE to enable accurate counting of prevalence
    metabolizer_summary$metabolizer_present[is.na(metabolizer_summary$metabolizer_present)] <- FALSE
    
    # Count the frequency of TRUE and FALSE for metabolizer_present in each study
    metabolizer_summary_counts <- as.data.frame(xtabs(~ metabolizer_present + study_name, data = metabolizer_summary))
    
    metabolizer_summary_percentage <- metabolizer_summary %>%
      group_by(study_name) %>%
      summarise(percentage = (sum(metabolizer_present == TRUE)/sum(metabolizer_present == TRUE, metabolizer_present == FALSE, is.na(metabolizer_present))*100), num_subjects = n()) # Returns sum of TRUE values
      
    # Plot the number of samples in each study
    ggplot(metabolizer_summary, aes(x=study_name)) +
      geom_bar() +  theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1))
    
    # Plot absolute counts of metabolizer_present in each study with counts labeled on bar chart
    ggplot(metabolizer_summary, aes(x = study_name, fill = metabolizer_present)) + geom_bar() + theme_minimal() + theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) + geom_text(aes(label = ..count..), stat = "count", angle = 0, size = 2, vjust = 0, hjust=.5, position = position_dodge())
    # Plot relative counts of metabolizer_present in each study
    ggplot(metabolizer_summary, aes(x = study_name, fill = metabolizer_present)) + geom_bar(position = "fill") + theme_minimal() + theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) + ylab("proportion")
    
    ggplot(metabolizer_summary_percentage, aes(x = study_name, y = percentage)) + geom_point() + theme_bw() + theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1))
    #ggp + annotate(geom = "table", x = 9, y = 3, label = list(metabolizer_summary_percentage))
    
    write.csv(metabolizer_summary, file = "RP01-79_20230619_metabolizer_summary.csv")
    
#------------------------- Check if genes are present ----------------------------------------------------------
    
    
    # Create a metabolizer_present column that returns TRUE if any of the relative adundance values for E. aspara, E. bolteae, or E. citroniae are above 0
    cMD_stool$genes_present <- ifelse(
                                  cMD_stool$`UniRef90_C5EGQ4`>0 |
                                  cMD_stool$`UniRef90_A8RZR5`>0 |
                                  cMD_stool$`UniRef90_C5EGQ4|unclassified`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|unclassified`>0 |
                                  cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_asparagiforme`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_asparagiforme`>0 |  
                                  cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_bolteae`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_bolteae`>0 |
                                  cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_citroniae`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_citroniae`>0 |  
                                  cMD_stool$`UniRef90_C5EGQ4|g__Clostridiales_unclassified.s__Clostridiales_bacterium_1_7_47FAA`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|g__Clostridiales_unclassified.s__Clostridiales_bacterium_1_7_47FAA`>0 |
                                  cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_clostridioforme`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_clostridioforme`>0 |
                                  cMD_stool$`UniRef90_C5EGQ4|g__Hungatella.s__Hungatella_hathewayi`>0 |
                                  cMD_stool$`UniRef90_A8RZR5|g__Hungatella.s__Hungatella_hathewayi`>0,
                                  TRUE, FALSE)
    summary(cMD_stool$genes_present)
    
    # Create a new dataframe with id, study_name, metabolizer_present (TRUE or FALSE), and the realtive abundance of each species
    genes_summary <- data.frame(id = cMD_stool$id, study_name = cMD_stool$study_name, genes_present = cMD_stool$genes_present, number_reads = cMD_stool$number_reads, million_reads = cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_RPK = cMD_stool$`UniRef90_C5EGQ4`,
                                UniRef90_A8RZR5_RPK = cMD_stool$`UniRef90_A8RZR5`,
                                UniRef90_C5EGQ4_unclassified_RPK = cMD_stool$`UniRef90_C5EGQ4|unclassified`,
                                UniRef90_A8RZR5_unclassified_RPK = cMD_stool$`UniRef90_A8RZR5|unclassified`,
                                UniRef90_C5EGQ4_s__Clostridium_asparagiforme_RPK = cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_asparagiforme`,
                                UniRef90_A8RZR5_s__Clostridium_asparagiforme_RPK = cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_asparagiforme`,
                                UniRef90_A8RZR5_s__Clostridium_bolteae_RPK = cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_bolteae`,
                                UniRef90_C5EGQ4_s__Clostridium_bolteae_RPK = cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_bolteae`,
                                UniRef90_C5EGQ4_s__Clostridium_citroniae_RPK = cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_citroniae`,
                                UniRef90_A8RZR5_s__Clostridium_citroniae_RPK = cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_citroniae`,
                                UniRef90_C5EGQ4_s__Hungatella_hathewayi_RPK = cMD_stool$`UniRef90_C5EGQ4|g__Hungatella.s__Hungatella_hathewayi`,
                                UniRef90_A8RZR5_s__Hungatella_hathewayi_RPK = cMD_stool$`UniRef90_A8RZR5|g__Hungatella.s__Hungatella_hathewayi`,
                                UniRef90_C5EGQ4_s__Clostridiales_bacterium_1_7_47FAA_RPK = cMD_stool$`UniRef90_C5EGQ4|g__Clostridiales_unclassified.s__Clostridiales_bacterium_1_7_47FAA`,
                                UniRef90_A8RZR5_s__Clostridiales_bacterium_1_7_47FAA_RPK = cMD_stool$`UniRef90_A8RZR5|g__Clostridiales_unclassified.s__Clostridiales_bacterium_1_7_47FAA`,
                                UniRef90_C5EGQ4_s__Clostridium_clostridioforme_RPK = cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_clostridioforme`,
                                UniRef90_A8RZR5_s__Clostridium_clostridioforme_RPK = cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_clostridioforme`,
                                
                                UniRef90_C5EGQ4_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_unclassified_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|unclassified`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_unclassified_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|unclassified`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_s__Clostridium_asparagiforme_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_asparagiforme`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_s__Clostridium_asparagiforme_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_asparagiforme`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_s__Clostridium_bolteae_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_bolteae`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_s__Clostridium_bolteae_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_bolteae`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_s__Clostridium_citroniae_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_citroniae`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_s__Clostridium_citroniae_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_citroniae`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_s__Hungatella_hathewayi_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|g__Hungatella.s__Hungatella_hathewayi`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_s__Hungatella_hathewayi_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|g__Hungatella.s__Hungatella_hathewayi`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_s__Clostridiales_bacterium_1_7_47FAA_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|g__Clostridiales_unclassified.s__Clostridiales_bacterium_1_7_47FAA`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_s__Clostridiales_bacterium_1_7_47FAA_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|g__Clostridiales_unclassified.s__Clostridiales_bacterium_1_7_47FAA`) / cMD_stool$number_reads/1000000,
                                UniRef90_C5EGQ4_s__Clostridium_clostridioforme_RPKM = as.numeric(cMD_stool$`UniRef90_C5EGQ4|g__Lachnoclostridium.s__Clostridium_clostridioforme`) / cMD_stool$number_reads/1000000,
                                UniRef90_A8RZR5_s__Clostridium_clostridioforme_RPKM = as.numeric(cMD_stool$`UniRef90_A8RZR5|g__Lachnoclostridium.s__Clostridium_clostridioforme`) / cMD_stool$number_reads/1000000
                                )
    
    # Convert NA values in the genes_present column to FALSE to enable accurate counting of prevalence
    genes_summary$genes_present[is.na(genes_summary$genes_present)] <- FALSE
    
    # Count the frequency of TRUE and FALSE for genes_present in each study
    genes_summary_counts <- as.data.frame(xtabs(~ genes_present + study_name, data = genes_summary))
    
    genes_summary_percentage <- genes_summary %>%
      group_by(study_name) %>%
      summarise(percentage = sum(genes_present == TRUE)/n()*100, num_subjects = n()) # Returns sum of TRUE values

#===================== PLOTS ===============================================
    
    # Plot the number of samples in each study
    ggplot(genes_summary, aes(x=study_name)) +
      geom_bar() +  theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1))
    
    # Plot absolute counts of genes_present in each study with counts labeled on bar chart
    ggplot(genes_summary, aes(x = study_name, fill = genes_present)) + geom_bar() + theme_minimal() + theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) + geom_text(aes(label = ..count..), stat = "count", angle = 0, size = 2, vjust = 0, hjust=.5, position = position_dodge())
    
    # Plot relative counts of genes_present in each study
    ggplot(genes_summary, aes(x = study_name, fill = genes_present)) + geom_bar(position = "fill") + theme_minimal() + theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) + ylab("proportion")
    
    # Plot the prevalence of metabolizers in each study
    ggplot(data = metabolizer_summary_percentage, aes(x=percentage, y=study_name)) +
      geom_bar(position = 'dodge', stat = 'identity', fill = "red2", color = "black") +
      xlim(0, 100) +
      theme_bw() +
      geom_text(aes(label=num_subjects), position = position_dodge2(width = 0.9), size = 2.5, vjust = 0.4, hjust = -0.25)
    
    # Plot the prevalence of at least one ucd gene from any source in each study
    ggplot(data = genes_summary_percentage, aes(x=percentage, y=study_name)) +
      geom_bar(position = 'dodge', stat = 'identity', fill = "grey", color = "black") +
      xlim(0, 100) +
      theme_bw() +
      geom_text(aes(label=num_subjects), position = position_dodge2(width = 0.9), size = 2.5, vjust = 0.4, hjust = -0.25)
    
    # Merge percentages (prevalence) of boi and goi into a single table to plot side-by-side columns
    metabolizer_summary_percentage$type <- "bacteria"
    genes_summary_percentage$type <- "gene"
    merged_percentages <- rbind(metabolizer_summary_percentage, genes_summary_percentage)
    
    ggplot(data = merged_percentages, aes(x=percentage, y=study_name, fill=type)) +
      geom_bar(position = position_dodge(), stat = 'identity') +
      xlim(0, 100) +
      theme_bw() +
      xlab("Prevalence (%)") +
      ylab("Study") +
      scale_fill_manual(values=c("grey", "red2"))
    
    # Plot a heatmap of the number of samples in each study to merge with other graphs
    ggplot(data = merged_percentages, aes(x=0.1, y=study_name, fill = num_subjects)) +
      geom_tile() +
      scale_fill_gradient(low = "white", high = "black")
    
    # Only look at studies with > 200 individuals
    merged_percentages_large_studies <- filter(merged_percentages, num_subjects >= 200)
    
    ggplot(data = merged_percentages_large_studies, aes(x=percentage, y=study_name, fill=type)) +
      geom_bar(position = position_dodge(), stat = 'identity') +
      xlim(0, 100) +
      theme_bw() +
      xlab("Prevalence (%)") +
      ylab("Study") +
      scale_fill_manual(values=c("grey", "red2"))
    
    # Plot a heatmap of the number of samples in each study to merge with other graphs
    ggplot(data = merged_percentages_large_studies, aes(x=0.1, y=study_name, fill = num_subjects)) +
      geom_tile() +
      scale_fill_gradient(low = "white", high = "black") +
      geom_text(aes(label=num_subjects), position = position_dodge2(width = 0.9), size = 2.5, vjust = 0.4, hjust = -0.25)
    
    setwd("C:/Users/Reilly/Desktop/")
    write.csv(merged_percentages_large_studies, file = "merged_percentages_large_studies.csv")    
    