This is additional R code for “Acidification from microbial oxidation of N, S, Fe and Mn as potential key mechanism for sprayed concrete deterioration in a subsea road tunnel”.

Here figures are made with the genome annotation of DRAM, FeGenie and DiSCo.

Libraries and data

Load R packages

suppressMessages(library(tidyverse))
suppressMessages(library(patchwork))

I load the data. Here we load the relative abundance, metadata, and GTDB taxonomy

#Some information about the samples.
metadata <- read.csv(file="data_R/metadata.txt", row.names=NULL, sep="\t") %>%
  mutate(Location = if_else(Location == "P", "P2", "T"))

#Relative abundance of MAGs
rel_abundance <- read.csv(file = "data_R/rel_abundance.txt", sep="\t", row.names=NULL)
#The sample ID in the publication is different that the sequencing ID
#We replace with the publication ID
colnames(rel_abundance)[2:9] <- metadata$ID

rel_abundance <- rel_abundance %>%
  rename(MAG = Genome) %>%
  pivot_longer(!MAG,
               names_to = "ID",
               values_to = "abundance")

#taxonomy estimated with GTDBTk
taxonomy <- read.csv(file="data_R/MAGs.taxonomy_R207.txt", row.names=NULL, sep="\t")
taxonomy <- taxonomy %>%
  rename(MAG = bin)

DRAM data

We also load the DRAM annotation. They are in two separate files The files are huge, so that make things easier I just import a few columns

#DRAM data
annotation_1 <- read_tsv(file="DRAM/annotations_1.tsv", 
                       col_select = c(1, fasta, kegg_id, kegg_hit))

print("List of MAGs:")
## [1] "List of MAGs:"
unique(annotation_1$fasta)
##   [1] "OFTM1"   "OFTM10"  "OFTM100" "OFTM101" "OFTM102" "OFTM103" "OFTM11" 
##   [8] "OFTM110" "OFTM114" "OFTM115" "OFTM119" "OFTM120" "OFTM124" "OFTM125"
##  [15] "OFTM131" "OFTM132" "OFTM134" "OFTM135" "OFTM136" "OFTM137" "OFTM14" 
##  [22] "OFTM141" "OFTM142" "OFTM144" "OFTM146" "OFTM147" "OFTM149" "OFTM15" 
##  [29] "OFTM151" "OFTM153" "OFTM156" "OFTM158" "OFTM16"  "OFTM160" "OFTM166"
##  [36] "OFTM168" "OFTM169" "OFTM17"  "OFTM170" "OFTM172" "OFTM173" "OFTM174"
##  [43] "OFTM175" "OFTM177" "OFTM180" "OFTM181" "OFTM182" "OFTM183" "OFTM188"
##  [50] "OFTM190" "OFTM191" "OFTM195" "OFTM196" "OFTM20"  "OFTM201" "OFTM205"
##  [57] "OFTM209" "OFTM211" "OFTM213" "OFTM214" "OFTM216" "OFTM217" "OFTM219"
##  [64] "OFTM22"  "OFTM220" "OFTM221" "OFTM222" "OFTM224" "OFTM226" "OFTM227"
##  [71] "OFTM228" "OFTM23"  "OFTM230" "OFTM234" "OFTM235" "OFTM236" "OFTM237"
##  [78] "OFTM238" "OFTM239" "OFTM244" "OFTM248" "OFTM249" "OFTM25"  "OFTM250"
##  [85] "OFTM251" "OFTM254" "OFTM255" "OFTM256" "OFTM257" "OFTM26"  "OFTM260"
##  [92] "OFTM261" "OFTM264" "OFTM265" "OFTM267" "OFTM268" "OFTM271" "OFTM274"
##  [99] "OFTM275" "OFTM276" "OFTM279" "OFTM28"  "OFTM283" "OFTM285" "OFTM286"
## [106] "OFTM288" "OFTM291" "OFTM293" "OFTM294" "OFTM295" "OFTM296" "OFTM297"
## [113] "OFTM298" "OFTM301" "OFTM303" "OFTM304" "OFTM305" "OFTM307" "OFTM308"
## [120] "OFTM31"  "OFTM310" "OFTM311" "OFTM312" "OFTM313" "OFTM316" "OFTM318"
## [127] "OFTM319" "OFTM320" "OFTM321" "OFTM322" "OFTM324" "OFTM33"  "OFTM331"
## [134] "OFTM333" "OFTM334" "OFTM335" "OFTM336" "OFTM337" "OFTM339" "OFTM34" 
## [141] "OFTM340" "OFTM341" "OFTM342" "OFTM345" "OFTM346" "OFTM347" "OFTM348"
## [148] "OFTM349" "OFTM350" "OFTM351" "OFTM352" "OFTM353" "OFTM354" "OFTM356"
## [155] "OFTM361" "OFTM363" "OFTM366" "OFTM367" "OFTM368" "OFTM370" "OFTM371"
## [162] "OFTM372" "OFTM377" "OFTM378" "OFTM379" "OFTM380" "OFTM381" "OFTM383"
## [169] "OFTM385" "OFTM386" "OFTM387" "OFTM388" "OFTM389" "OFTM39"  "OFTM392"
## [176] "OFTM393" "OFTM397" "OFTM399" "OFTM4"   "OFTM40"  "OFTM401" "OFTM42" 
## [183] "OFTM44"  "OFTM46"  "OFTM47"  "OFTM49"  "OFTM5"   "OFTM51"  "OFTM52" 
## [190] "OFTM55"  "OFTM6"   "OFTM61"  "OFTM64"  "OFTM65"  "OFTM7"   "OFTM71" 
## [197] "OFTM72"  "OFTM74"  "OFTM75"  "OFTM76"  "OFTM77"  "OFTM78"  "OFTM8"  
## [204] "OFTM81"  "OFTM82"  "OFTM83"  "OFTM87"  "OFTM89"  "OFTM9"   "OFTM91" 
## [211] "OFTM94"  "OFTM95"  "OFTM96"  "OFTM97"  "OFTM99"
annotation_2 <- read_tsv(file="DRAM/annotations_2.tsv", 
                       col_select = c(1, fasta, kegg_id, kegg_hit))

print("List of MAGs:")
## [1] "List of MAGs:"
unique(annotation_2$fasta)
##   [1] "OFTM104" "OFTM105" "OFTM106" "OFTM107" "OFTM108" "OFTM109" "OFTM111"
##   [8] "OFTM112" "OFTM113" "OFTM116" "OFTM117" "OFTM118" "OFTM12"  "OFTM121"
##  [15] "OFTM122" "OFTM123" "OFTM126" "OFTM127" "OFTM128" "OFTM129" "OFTM13" 
##  [22] "OFTM130" "OFTM133" "OFTM138" "OFTM139" "OFTM140" "OFTM143" "OFTM145"
##  [29] "OFTM148" "OFTM150" "OFTM152" "OFTM154" "OFTM155" "OFTM157" "OFTM159"
##  [36] "OFTM161" "OFTM162" "OFTM163" "OFTM164" "OFTM165" "OFTM167" "OFTM171"
##  [43] "OFTM176" "OFTM178" "OFTM179" "OFTM18"  "OFTM184" "OFTM185" "OFTM186"
##  [50] "OFTM187" "OFTM189" "OFTM19"  "OFTM192" "OFTM193" "OFTM194" "OFTM197"
##  [57] "OFTM198" "OFTM199" "OFTM2"   "OFTM200" "OFTM202" "OFTM203" "OFTM204"
##  [64] "OFTM206" "OFTM207" "OFTM208" "OFTM21"  "OFTM210" "OFTM212" "OFTM215"
##  [71] "OFTM218" "OFTM223" "OFTM225" "OFTM229" "OFTM231" "OFTM232" "OFTM233"
##  [78] "OFTM24"  "OFTM240" "OFTM241" "OFTM242" "OFTM243" "OFTM245" "OFTM246"
##  [85] "OFTM247" "OFTM252" "OFTM253" "OFTM258" "OFTM259" "OFTM262" "OFTM263"
##  [92] "OFTM266" "OFTM269" "OFTM27"  "OFTM270" "OFTM272" "OFTM273" "OFTM277"
##  [99] "OFTM278" "OFTM280" "OFTM281" "OFTM282" "OFTM284" "OFTM287" "OFTM289"
## [106] "OFTM29"  "OFTM290" "OFTM292" "OFTM299" "OFTM3"   "OFTM30"  "OFTM300"
## [113] "OFTM302" "OFTM306" "OFTM309" "OFTM314" "OFTM315" "OFTM317" "OFTM32" 
## [120] "OFTM323" "OFTM325" "OFTM326" "OFTM327" "OFTM328" "OFTM329" "OFTM330"
## [127] "OFTM332" "OFTM338" "OFTM343" "OFTM344" "OFTM35"  "OFTM355" "OFTM357"
## [134] "OFTM358" "OFTM359" "OFTM36"  "OFTM360" "OFTM362" "OFTM364" "OFTM365"
## [141] "OFTM369" "OFTM37"  "OFTM373" "OFTM374" "OFTM375" "OFTM376" "OFTM38" 
## [148] "OFTM382" "OFTM384" "OFTM390" "OFTM391" "OFTM394" "OFTM395" "OFTM396"
## [155] "OFTM398" "OFTM400" "OFTM41"  "OFTM43"  "OFTM45"  "OFTM48"  "OFTM50" 
## [162] "OFTM53"  "OFTM54"  "OFTM56"  "OFTM57"  "OFTM58"  "OFTM59"  "OFTM60" 
## [169] "OFTM62"  "OFTM63"  "OFTM66"  "OFTM67"  "OFTM68"  "OFTM69"  "OFTM70" 
## [176] "OFTM73"  "OFTM79"  "OFTM80"  "OFTM84"  "OFTM85"  "OFTM86"  "OFTM88" 
## [183] "OFTM90"  "OFTM92"  "OFTM93"  "OFTM98"

I do some data wrangling

#I only care about the ones with kegg_id
annotation_1 <- drop_na(annotation_1, kegg_id)
annotation_2 <- drop_na(annotation_2, kegg_id)

#put all together
annotation_big <- bind_rows(annotation_1, annotation_2)
rm(annotation_1, annotation_2)

#change column names
annotation_big <- annotation_big  %>%
  rename(gene = 1) %>%
  rename(MAG = fasta) 

head(annotation_big)

FeGenie data

We load the FeGenie data and do some data wrangling

data_iron <- read.csv(file="fegenie_out/FeGenie-geneSummary.csv", 
                      row.names=NULL, sep=",", header=TRUE, skip = 1) %>%
  select(1:4) %>%
  #we only include some categories
  filter(category %in% c("iron_oxidation", 
                         "possible_iron_oxidation_and_possible_iron_reduction",
                         "iron_reduction", 
                         "probable_iron_reduction")) %>%
  rename(gene_name = HMM) %>%
  rename(gene = orf) %>%
  #get the MAG name from the gene id
  separate(gene, c("MAG", NA, NA), sep="_", remove = FALSE)  %>%
  #some genes are renamed
  mutate(gene_name = case_when(gene_name == "MtrB_TIGR03509" ~ "MtrB",
                               gene_name == "MtrC_TIGR03507" ~ "MtrC",
    grepl("[:{:]pc", gene_name) ~ "pc",
                              grepl("[:{:]porin", gene_name) ~ "porin",
                              grepl("Cyc2_", gene_name) ~ "Cyc2",
                              grepl("[:{:]omc", gene_name) ~ "omc",
                            TRUE ~ as.character(gene_name))) %>%
  select(-genome.assembly)

head(data_iron)

DiSCo data

In R load the data from DiSCo v1.0.0

DiSCo_result <- read_tsv(file="Disco_DRAM/OFTMs_DRAM_genes.DiSCo.filtered.txt",
                         skip = 1,
                         col_names =  c("gene", "protein", "type"),
                         col_select = c(1, 2, 3)) %>%
  #simplify the gene names
  separate(protein, c("gene_name", NA, NA), sep="_", fill = "right")  %>%
  #get the MAGs name from the gene name
  separate(gene, c("MAG", NA, NA, NA), sep="_", fill = "right")

head(DiSCo_result)

Osmotic stress

How many MAGs have ectC (L-ectoine synthase), ectD (ectoine hydroxylase) or betB (betaine aldehyde dehydrogenase)?

annotation_big %>%
  mutate(gene_name = case_when(
    str_detect(kegg_id, "K06720") ~ "ectC",
    str_detect(kegg_id, "K10674") ~ "ectD",
    str_detect(kegg_id, "K00130") ~ "betB",
    TRUE ~ "others"))%>%
  filter(gene_name != "others") %>%
  distinct(MAG, gene_name) %>%
  group_by(gene_name) %>%
  summarise(number_of_MAGs = n())

Iron cycling

We use the FeGenie data

Some Iron-cycling stats

How many MAGs had a potential for iron oxidation?

data_iron %>%
  filter(category == "iron_oxidation") %>%
  distinct(MAG) %>%
  nrow()
## [1] 62

What was abundance of the MAGs with potential for iron oxidation?

data_iron %>%
  filter(category == "iron_oxidation") %>%
  distinct(MAG) %>%
  inner_join(rel_abundance, by="MAG")  %>%
  group_by(ID)  %>%
  summarise(my_sum = sum(abundance)) %>%
  summarise(min = min(my_sum),
            max = max(my_sum))

How many MAGs had a potential for iron reduction?

data_iron %>%
  filter(category == "iron_reduction") %>%
  distinct(MAG) %>%
  nrow()
## [1] 34

What was abundance of the MAGs with potential for iron reduction?

data_iron %>%
  filter(category == "iron_reduction") %>%
  distinct(MAG) %>%
  inner_join(rel_abundance, by="MAG")  %>%
  group_by(ID)  %>%
  summarise(my_sum = sum(abundance)) %>%
  summarise(min = min(my_sum),
            max = max(my_sum))

How many MAGs have key genes linked to iron cycling?

#we make simpler the Fegeneie dataset 
data_iron_simple <- data_iron %>%
  #we only show some genes
  filter(gene_name %in% c("Cyc1", "Cyc2", "MtoA", "MtrA", "MtrB", "MtrC")) %>%
  #for the plots I only care if a gene is a present in a MAG
  group_by(MAG, gene_name) %>%
  summarise(copies = n())
## `summarise()` has grouped output by 'MAG'. You can override using the `.groups`
## argument.
data_iron_simple %>%
  group_by(gene_name) %>%
  summarise(number_of_MAGs = n())

Let’s look a bit more detail on Cyc2. For each phyla, how many MAGs have Cyc2?

data_iron_simple %>%
  filter(gene_name == "Cyc2") %>%
  select(MAG)  %>%
  inner_join(taxonomy, by="MAG") %>%
  group_by(Phylum) %>%
  summarise(number_of_MAGs = n()) %>%
  arrange(desc(number_of_MAGs))

How many of the MAGs with Cyc2 are Zetaproteobacteria?

data_iron_simple %>%
  filter(gene_name == "Cyc2") %>%
  select(MAG)  %>%
  inner_join(taxonomy, by="MAG") %>%
  filter(Class == "c__Zetaproteobacteria") %>%
  nrow()
## [1] 5

More details on MtoA. For each phyla, how many MAGs have MtoA. ?

data_iron_simple %>%
  filter(gene_name == "MtoA") %>%
  select(MAG)  %>%
  inner_join(taxonomy, by="MAG") %>%
  group_by(Phylum) %>%
  summarise(number_of_MAGs = n()) %>%
  arrange(desc(number_of_MAGs))

How many have mtrABC?

data_iron_simple %>%
  filter(gene_name %in% c("MtrA", "MtrB", "MtrC")) %>%
  pivot_wider(names_from = gene_name, values_from = copies) %>%
  drop_na() %>%
  nrow()
## [1] 11

What is the taxonomy of the ones with mtrABC?

data_iron_simple %>%
  filter(gene_name %in% c("MtrA", "MtrB", "MtrC")) %>%
  pivot_wider(names_from = gene_name, values_from = copies) %>%
  drop_na() %>%
  inner_join(taxonomy, by="MAG")

Of these how many are in the family Woeseiaceae?

data_iron_simple %>%
  filter(gene_name %in% c("MtrA", "MtrB", "MtrC")) %>%
  pivot_wider(names_from = gene_name, values_from = copies) %>%
  drop_na() %>%
  inner_join(taxonomy, by="MAG") %>%
  filter(Family == "f__Woeseiaceae") %>%
  nrow()
## [1] 6

Iron figures

Number of MAGs with iron oxidation genes

#I add taxonomic data
plot_iron <- data_iron_simple %>%
  inner_join(taxonomy, by="MAG") %>%
   mutate(Phylum = case_when(
     Phylum == "p__Proteobacteria"~ "Proteobacteria",
     Phylum == "p__Nitrospinota"~ "Nitrospinota",
     Phylum == "p__Nitrospirota"~ "Nitrospirota",
     Phylum == "p__Planctomycetota" ~ "Planctomycetota",
    TRUE ~ "Other Phyla"))

#order of the colours in the plot
plot_iron$Phylum = factor(as.character(plot_iron$Phylum), 
                             levels=c("Nitrospinota",
                                      "Nitrospirota",
                                      "Planctomycetota",
                                      "Proteobacteria",
                                      "Other Phyla"
                            ))


my_palette  <- palette.colors(palette = "Tableau 10")
names(my_palette) <- NULL

ggplot(plot_iron, aes(gene_name)) +
  geom_bar(aes(fill=Phylum), colour="white", linewidth=0.2) +
   theme_bw() +
    theme(legend.title = element_blank(),
        axis.text = element_text(colour = "black", size=7),
        axis.title = element_text(colour = "black", size=7),
        legend.text = element_text(colour = "black", size=7),
        axis.title.y = element_text(colour = "black", size=7),
        panel.border = element_blank(),
        panel.grid = element_blank(),
        legend.position = "right",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        legend.margin=margin(0,0,0,0)) +
   scale_fill_manual(values=my_palette) +
  labs(y="Number of MAGs", x="Gene Name") +
  guides(fill = guide_legend(keywidth = 0.7, keyheight = 0.7))

#ggsave("Graphs/Gene_Phylum_iron.svg", width = 80, height = 45, units = "mm", dpi = 600)

Relative abundance of MAGs with iron oxidation genes.

plot_iron <- data_iron_simple %>%
  distinct(MAG, gene_name) %>% #just in case I have duplicates
  #Add relative abundance info
  inner_join(rel_abundance, by="MAG") %>%
  select(gene_name, MAG, ID, abundance) %>%
  group_by(ID, gene_name) %>%
  summarise(av = sum(abundance)) %>%
  inner_join(metadata, by="ID")
## Warning in inner_join(., rel_abundance, by = "MAG"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 57 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
ggplot(plot_iron, aes(gene_name, av, colour=Location, shape=as.character(Year))) +
  geom_jitter(height = 0, width = 0.05) +
  theme_bw() +
     theme(legend.title = element_text(colour = "black", size=8, face="bold"),
           panel.grid =  element_blank(),
           strip.background = element_blank(),
           strip.text = element_text(colour = "black", size=7, face="bold"),
        axis.text = element_text(colour = "black", size=7),
        axis.text.x = element_text(face="italic"),
        legend.text = element_text(colour = "black", size=7),
        axis.title = element_text(colour = "black", size=7),
        legend.position = "right",
        legend.margin=margin(0,0,0,0)) +
  scale_colour_viridis_d(begin = 0.2, end=0.8, name="Location", direction = -1) +
  scale_shape(name="Year") +
  labs(y="MAG abudance (%)", x="Gene name") +
  guides(x=guide_axis(angle = 90), 
         shape=guide_legend(keywidth = 0.6, keyheight = 0.6),
         colour=guide_legend(keywidth = 0.6, keyheight = 0.6))

#ggsave("Graphs/Gene_abundance_iron2.pdf", width = 70, height = 50, units = "mm", dpi = 600)

Nitrogen cyling

To make the nitrogen figures I use the DRAM data. I chose the genes with their kegg codes.

data_nitrogen <- annotation_big %>%
  mutate(gene_name = case_when(
    str_detect(kegg_id, "K10944") ~ "amoA-pmoA",
    #str_detect(kegg_id, "K10535") ~ "hao",
    str_detect(kegg_id, "K20932") ~ "hzsA",
    str_detect(kegg_id, "K03385") ~ "nrfA",
    str_detect(kegg_id, "K20935") ~ "hdh",
    str_detect(kegg_id, "K00370") ~ "narG-nxrA",
    str_detect(kegg_id, "K02567") ~ "napA",
    str_detect(kegg_id, "K00368") ~ "nirK",
    str_detect(kegg_id, "K15864") ~ "nirS",
    str_detect(kegg_id, "K04561") ~ "norB",
    str_detect(kegg_id, "K00376") ~ "nosZ",
    TRUE ~ "others")) %>%
  filter(gene_name != "others") 

head(data_nitrogen)

Some Nitrogen stats

How many MAGs have a N-cycling gene?

data_nitrogen %>%
  distinct(MAG, gene_name)  %>%
   group_by(gene_name) %>%
  summarise(number_of_MAGs = n())

How many MAGs have napA and/or narG/nxrA?

data_nitrogen %>%
  distinct(MAG, gene_name)  %>%
  filter(gene_name %in% c("napA", "narG-nxrA")) %>%
  distinct(MAG) %>%
  nrow()
## [1] 191

How many MAGs have nirK and/or nirS?

data_nitrogen %>%
  distinct(MAG, gene_name)  %>%
  filter(gene_name %in% c("nirK", "nirS")) %>%
  distinct(MAG) %>%
  nrow()
## [1] 117

Nitrogen Figures

Number of MAGs with Nitrogen genes.

data_plot <- data_nitrogen %>%
  distinct(MAG, gene_name)  %>%
  group_by(gene_name) %>%
  summarise(av = n()) #count how many

plotA <- ggplot(data_plot, aes(gene_name, av,)) +
  geom_col() +
  theme_bw() +
     theme(legend.title = element_text(colour = "black", size=8, face="bold"),
           panel.grid =  element_blank(),
           strip.background = element_blank(),
           strip.text = element_text(colour = "black", size=7, face="bold"),
        axis.text = element_text(colour = "black", size=7),
        axis.text.x = element_blank(),
        legend.text = element_text(colour = "black", size=7),
        axis.title = element_text(colour = "black", size=7),
        legend.position = "right",
        axis.title.x = element_blank(),
        legend.margin=margin(0,0,0,0)) +
  scale_colour_viridis_d(begin = 0.2, end=0.8, name="Location", direction=-1) +
  scale_shape(name="Year") +
  labs(y="Number of MAGs", x="Gene name") +
  guides(x=guide_axis(angle = 90), 
         shape=guide_legend(keywidth = 0.6, keyheight = 0.6),
         colour=guide_legend(keywidth = 0.6, keyheight = 0.6))

plotA

Relative abundance of MAGs with Nitrogen genes

data_plot <- data_nitrogen %>% 
    distinct(MAG, gene_name) %>% 
  inner_join(rel_abundance, by="MAG") %>% #add relative abundance info
  select(gene_name, MAG, ID, abundance) %>%
  group_by(ID, gene_name) %>%
  summarise(av = sum(abundance)) %>%
  inner_join(metadata, by="ID")
## Warning in inner_join(., rel_abundance, by = "MAG"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 9 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
plotB <- ggplot(data_plot, aes(gene_name, av, colour=Location, shape=as.character(Year))) +
  geom_jitter(height = 0, width = 0.05) +
  theme_bw() +
     theme(legend.title = element_text(colour = "black", size=8, face="bold"),
           panel.grid =  element_blank(),
           strip.background = element_blank(),
           strip.text = element_text(colour = "black", size=7, face="bold"),
        axis.text = element_text(colour = "black", size=7),
        axis.text.x = element_text(face="italic"),
        legend.text = element_text(colour = "black", size=7),
        axis.title = element_text(colour = "black", size=7),
        legend.position = "right",
        legend.margin=margin(0,0,0,0)) +
  scale_colour_viridis_d(begin = 0.2, end=0.8, name="Location", direction=-1) +
  scale_shape(name="Year") +
  labs(y="MAG abudance (%)", x="Gene name") +
  guides(x=guide_axis(angle = 90), 
         shape=guide_legend(keywidth = 0.6, keyheight = 0.6),
         colour=guide_legend(keywidth = 0.6, keyheight = 0.6))

plotB

We put the two figures together

plotA / plotB + 
  plot_annotation(tag_levels = 'a') & 
  theme(plot.tag = element_text(size = 8))

#ggsave("Graphs/Gene_abundance_nitrogen.png", width = 90, height = 90, units = "mm", dpi = 600)

Sulfur figure

For the plots of Sulfur oxidation and reduction I use the DRAM data. I chose the genes with their kegg codes.

data_sulfur <- annotation_big %>%
  mutate(gene_name = case_when(
    str_detect(kegg_id, "K11180") ~ "dsrA",
    str_detect(kegg_id, "K11181") ~ "dsrB",
    str_detect(kegg_id, "K00394") ~ "aprA",
    str_detect(kegg_id, "K00395") ~ "aprB",
    str_detect(kegg_id, "K17222") ~ "SoxA",
    str_detect(kegg_id, "K17224") ~ "SoxB",
    str_detect(kegg_id, "K08352") ~ "phsA",
    str_detect(kegg_id, "K08358") ~ "ttrB",
    str_detect(kegg_id, "K21307") ~ "soeA",
    str_detect(kegg_id, "K05301") ~ "sorA",
    str_detect(kegg_id, "K17218") ~ "sqr",
    TRUE ~ "others")) %>%
  filter(gene_name != "others") 

head(data_sulfur)

Some sulfur stats

Number of MAGs with a sulfur cycling gene

data_sulfur  %>%
  distinct(MAG, gene_name)  %>%
   group_by(gene_name) %>%
  summarise(number_of_MAGs = n())

Abundance of MAGs for each gene

func_sulfur_ab <- function(var_gene_name) {
  qr <- data_sulfur  %>%
  filter(gene_name == var_gene_name) %>%
  distinct(MAG) %>%
  inner_join(rel_abundance, by="MAG")  %>%
  group_by(ID)  %>%
  summarise(my_sum = sum(abundance)) %>%
  summarise(min = min(my_sum),
            max = max(my_sum))
  return(qr)
}

sulfur_genes <- sort(unique(data_sulfur$gene_name))

sapply(sulfur_genes, func_sulfur_ab)
##     aprA     aprB     dsrA      dsrB      phsA     soeA     sorA      SoxA    
## min 3.871024 4.084233 0.3212479 0.3212479 1.939891 6.525422 0.4180097 2.461686
## max 16.03586 16.51732 5.260237  5.260237  8.99131  8.915912 8.5432    6.302103
##     SoxB     sqr      ttrB    
## min 2.469701 9.31766  4.009362
## max 6.355511 16.12243 12.24704

Figure sulfur

First I made a figure with the relative abundance of MAGs having sulfur oxidation or reduction genes.

data_plot <- data_sulfur  %>%
  distinct(MAG, gene_name) %>% #I only care about presence of a gene
  inner_join(rel_abundance, by="MAG") %>% #add rel abundance
  select(gene_name, MAG, ID, abundance) %>%
  group_by(ID, gene_name) %>%
  summarise(av = sum(abundance)) %>%
  inner_join(metadata, by="ID")
## Warning in inner_join(., rel_abundance, by = "MAG"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 9 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
#Order of the x-axis
data_plot$gene_name = factor(as.character(data_plot$gene_name), 
                             levels=c("dsrA",
                                      "dsrB",
                                      "aprA",
                                      "aprB",
                                      "soeA",
                                      "sorA",
                                      "SoxA",
                                      "SoxB",
                                      "sqr",
                                      "ttrB",
                                      "phsA"
                            ))

ggplot(data_plot, aes(gene_name, av, colour=Location, shape=as.character(Year))) +
  geom_jitter(height = 0, width = 0.05) +
  theme_bw() +
     theme(legend.title = element_text(colour = "black", size=8, face="bold"),
           panel.grid =  element_blank(),
           strip.background = element_blank(),
           strip.text = element_text(colour = "black", size=7, face="bold"),
        #text = element_text(colour = "black", size=7),
        axis.text = element_text(colour = "black", size=7),
        axis.text.x = element_text(face = "italic"),
        #axis.title = element_text(colour = "black", size=7),
        legend.text = element_text(colour = "black", size=7),
        #legend.position = "bottom",
        axis.title = element_text(colour = "black", size=7),
        legend.position = "right",
        #axis.title.x = element_blank(),
        legend.margin=margin(0,0,0,0)) +
  scale_colour_viridis_d(begin = 0.2, end=0.8, name="Location", direction=-1) +
  scale_shape(name="Year") +
  labs(y="MAG abudance (%)", x="Genes") +
  guides(x=guide_axis(angle = 90), 
         shape=guide_legend(keywidth = 0.6, keyheight = 0.6),
         colour=guide_legend(keywidth = 0.6, keyheight = 0.6))

#ggsave("Graphs/Gene_abundance_sulfur2.pdf", width = 70, height = 55, units = "mm", dpi = 600)

Now I make a plot with the number of MAGs having a gene

data_plot <- data_sulfur %>%
  distinct(MAG, gene_name) %>% #I only care about presence of a gene 
  inner_join(taxonomy, by="MAG") %>% #add taxonomy
   mutate(Phylum = case_when(
     Phylum == "p__Proteobacteria"~ "Proteobacteria",
     Phylum == "p__Nitrospirota"~ "Nitrospirota",
     Phylum == "p__Planctomycetota" ~ "Planctomycetota",
     Phylum == "p__Bacteroidota" ~ "Bacteroidota",
     Phylum == "p__Actinobacteriota" ~ "Actinobacteriota",
    TRUE ~ "Other Phyla"))

#Order of the x-axis
data_plot$gene_name = factor(as.character(data_plot$gene_name), 
                             levels=c("dsrA",
                                      "dsrB",
                                      "aprA",
                                      "aprB",
                                      "soeA",
                                      "sorA",
                                      "SoxA",
                                      "SoxB",
                                      "sqr",
                                      "ttrB",
                                      "phsA"
                            ))

data_plot$Phylum = factor(as.character(data_plot$Phylum), 
                             levels=c("Actinobacteriota",
                                      "Bacteroidota",
                                      "Nitrospirota",
                                      "Planctomycetota",
                                      "Proteobacteria",
                                      "Other Phyla"
                            ))

my_palette  <- palette.colors(palette = "Tableau 10")
names(my_palette) <- NULL

ggplot(data_plot, aes(gene_name)) +
  geom_bar(aes(fill=Phylum), colour="white", linewidth=0.2) +
   theme_bw() +
    theme(legend.title = element_blank(),
        axis.text = element_text(colour = "black", size=7),
        axis.title = element_text(colour = "black", size=7),
        legend.text = element_text(colour = "black", size=7),
        axis.title.y = element_text(colour = "black", size=7),
        panel.border = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        legend.margin=margin(0,0,0,0)) +
   scale_fill_manual(values=my_palette) +
  labs(y="Number of MAGs", x="Gene Name") +
  guides(fill = guide_legend(nrow = 3, keywidth = 0.6, keyheight = 0.6))

#ggsave("Graphs/Gene_Phylum_sulfur.svg", width = 70, height = 60, units = "mm", dpi = 600)

Dsr accessory genes

We want to make a plot showing the presence/absence of SOX genes and dsr accessory genes across MAGs. We also look at the presence/absence of soeA and sorA genes.

First we select the MAGs with sulfur oxidation genes from the DRAM data

#Only Sox, soeA or sorA
data_sulfur <- annotation_big %>%
  mutate(gene_name = case_when(
    str_detect(kegg_id, "K17222") ~ "SoxA",
    str_detect(kegg_id, "K17224") ~ "SoxB",
    str_detect(kegg_id, "K17225") ~ "SoxC",
    str_detect(kegg_id, "K22622") ~ "SoxD",
    str_detect(kegg_id, "K17223") ~ "SoxX",
    str_detect(kegg_id, "K17226") ~ "SoxY",
    str_detect(kegg_id, "K17227") ~ "SoxZ",
    str_detect(kegg_id, "K21307") ~ "soeA",
    str_detect(kegg_id, "K05301") ~ "sorA",
    TRUE ~ "others"))%>%
  filter(gene_name != "others") %>%
  select(MAG, gene_name) %>%
  distinct() %>%
  mutate(type = case_when(
    str_detect(gene_name, "Sox") ~ "SOX",
    gene_name == "soeA" ~ "soeA",
    gene_name == "sorA" ~ "soeA"))

To identity dsr accessory genes we use the DisCo results. For identification of dsr and apr genes, we use DisCo instead of DRAM. This is because in DisCo there is a classification of these genes as reductive or oxidative.

To keep the plot simple, we only show MAGs having aprAB or dsrAB, or the the ones with soxAB

#Create a filter of only MAGs with SoxA, SoxB
data_sulfur_filter <- data_sulfur %>%
  filter(gene_name %in% c("SoxA", "SoxB")) %>%
  distinct(MAG) %>%
  pull(MAG)

#Create a filter of MAGs where DisCo found AprA, AprB, DsrA or DsrB
DiSCo_result_filter <- DiSCo_result %>%
  filter(gene_name %in% c("AprA", "AprB", "DsrA", "DsrB")) %>%
  distinct(MAG) %>%
  pull(MAG)

large_filter <- union(DiSCo_result_filter, data_sulfur_filter)
xr <- DiSCo_result %>% #load the DisCo data
  bind_rows(data_sulfur) %>% #add the DRAM sulfur data
  filter(MAG %in% large_filter)  %>% #keep only MAGs with aprAB or dsrAB or soxAB
  left_join(taxonomy, by="MAG") %>% #add taxonomy
  #if unclasfied make a NA
  mutate(Species = na_if(Species, "s__")) %>%
  mutate(Genus = na_if(Genus, "g__")) %>%
  mutate(Family = na_if(Family, "f__")) %>%
  mutate(Order = na_if(Order, "o__")) %>%
  mutate(type = case_when(type == "oxidative" ~ "dsr oxidative",
                          type == "reductive" ~ "dsr reductive",
                          type == "unspecific" ~ "dsr unspecific",
                          TRUE ~ as.character(type)))
#Y-axis names
xr$Species <- as.character(xr$Species)
xr$Genus <- as.character(xr$Genus)
xr$Family <- as.character(xr$Family)
xr$Order <- as.character(xr$Order)
xr$Class <- as.character(xr$Class)
#Get the lowest classifed taxonomic rank
xr$S2 <- xr$Species %>% coalesce(xr$Species, xr$Genus, xr$Family, xr$Order, xr$Class)
xr$S2 <- paste(xr$S2, xr$MAG, sep = ", ")
xr$S3 <- paste(xr$Phylum, xr$Class, xr$S2, sep = ", ")

xr <- xr %>%
  arrange(S3) %>% 
  mutate(S2 = factor(S2, unique(S2)))

We make the plot

ggplot(xr, aes(gene_name, S2, fill=type)) +
  geom_tile(colour="white", linewidth=0.2) +
  theme_bw() +
     theme(legend.title = element_text(colour = "black", size=7, face="bold"),
        axis.text = element_text(colour = "black", size=7),
        legend.text = element_text(colour = "black", size=7),
        panel.border = element_blank(),
        axis.title = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.margin=margin(0,0,0,0)) +
  scale_fill_viridis_d(direction =1,
                       option = "D") +
  guides(fill=guide_legend(nrow = 2))

#ggsave("Graphs/DiSCo_v2.svg", width = 130, height = 175, units = "mm", dpi = 600, limitsize = FALSE)

Table

It might be useful to create a big table summarizing all the data

data_nitrogen_export <- data_nitrogen %>%
  group_by(MAG, gene_name) %>%
  summarise(copies = n())

data_sulfur_export <- data_sulfur %>%
  group_by(MAG, gene_name) %>%
  summarise(copies = n()) 

data_osmotic <- annotation_big %>%
  mutate(gene_name = case_when(
    str_detect(kegg_id, "K06720") ~ "etcC",
    str_detect(kegg_id, "K10674") ~ "etcD",
    str_detect(kegg_id, "K00130") ~ "betB",
    TRUE ~ "others"))%>%
  filter(gene_name != "others") %>%
  group_by(MAG, gene_name) %>%
  summarise(copies = n())

data_biogechemical <- bind_rows(data_osmotic, 
                                data_sulfur_export, 
                                data_nitrogen_export, 
                                data_iron_simple) %>%
  pivot_wider(names_from = gene_name, values_from = copies, values_fill = 0) 

#relative abundance
data_abundance <-rel_abundance %>% 
  select(MAG, abundance, ID) %>%
  pivot_wider(names_from = ID, values_from = abundance)

#All together
data_multiple  <-data_abundance %>% 
  inner_join(taxonomy, by="MAG") %>%
  inner_join(data_biogechemical, by="MAG")  

#write_tsv(data_multiple, "summary_biogeochemical.tsv")