This is the 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 we look at the abundance of MAGs and to summarize the data into the file MAGs_summary_Oslofjord.tsv.

Linux stuff

Things we did before:

Trimming adapters and low quality bases (fastp)
Normalization (BBTools 38.61b)
Co-assembly (Megahit 1.2.9)
Mapping back to the co-assembly (Bowtie v2.3.5.1)
Binning (MetaBAT2 v2.15, BinSanity v0.5.3)
Completness and contamination (CheckM v1.1.3)
Dereplication (DASTool v1.1.2)

Which is described in detail on: Suarez et al 2022 (https://doi.org/10.1111/1462-2920.16006)

Taxonomy

I estimate taxonomy with GTDB-Tk v2.1.0

Note that in Suarez et al 2022, the GTDB R202 taxonomy was used. Here taxonomy is estimated using R207v2 taxonomy.

gtdbtk classify_wf --cpus 20 --extension fa \
--genome_dir MAGs_Oslofjord --out_dir gtdbtk_results \
--pplacer_cpus 10

I get two files with the taxonomy classification:gtdbtk.bac120.summary.tsv and gtdbtk.ar122.summary.tsv. There is some useful info, but here I only care about the first two columns.

cat gtdbtk.bac120.summary.tsv | grep -v "user_genome" | \
awk -v OFS='\t' '{print $1,$2}' > gtdbtk.taxonomy_edit.tsv
# same as above, but append the archaea
cat gtdbtk.ar122.summary.tsv | grep -v "user_genome" | \
awk -v OFS='\t' '{print $1,$2}' >> gtdbtk.taxonomy_edit.tsv
#Col names for the taxonomy file
echo -e "bin\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies"> MAGs.taxonomy.txt
#Separate classification column and append
sed 's/\;/\t/g' gtdbtk.taxonomy_edit.tsv >> MAGs.taxonomy.txt
rm gtdbtk.taxonomy_edit.tsv

Stuff we need in R

Now we do some R magic

We load some libraries I’ll need

library(tidyverse)

Data

We load the data into R

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

#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)

CheckM data

We load the results from checkM v1.1.3

checkm <- read.csv(file="data_R/CheckM.txt", row.names=NULL, sep="\t")

CheckM was done in Suarez et al 2022 prior to dereplication and assigning the final names to the MAGs. Therefore in the CheckM results the bins ID don’t mach the MAGs names. This is one example of the bins ID

head(checkm$Bin.Id)

## [1] "BinSanityLC-kmean-bin_0-bin_5"               
## [2] "BinSanityLC-kmean-bin_1-bin_4_sub"           
## [3] "BinSanityLC-kmean-bin_1-bin_5"               
## [4] "BinSanityLC-kmean-bin_11-bin_0-refined_0_sub"
## [5] "BinSanityLC-kmean-bin_11-bin_6-refined_1_sub"
## [6] "BinSanityLC-kmean-bin_13-bin_0-refined_1_sub"

Fortunately we have a file matching the bins IDs with the MAGs names

checkm_names <- read.csv(file="data_R/oldnames.txt", row.names=NULL, sep=" ", header=FALSE)
colnames(checkm_names) <- c("MAG", "Bin.Id")
head(checkm_names)

And we put it all together

checkm_names$MAG <- substr(checkm_names$MAG,1,nchar(checkm_names$MAG)-3)

checkm <- checkm %>%
  inner_join(checkm_names , by="Bin.Id")  %>%
  rename(contigs = X..genomes) %>%
  select(MAG, Completeness, Contamination)

Additional data

We also include some stats about the MAGs and the NCBI accessions.

#stats from bbtools
data_stats <- read_tsv(file = "data_R/stats.txt", col_select = c("n_contigs", "contig_bp", "filename")) %>%
  mutate(MAG = str_sub(str_replace(filename, '.*\\/', ''), end = -4)) %>%
  select(-filename)  %>%
  rename(bp_MAG = contig_bp)

#NCBI accessions
#I got this from https://www.ncbi.nlm.nih.gov/bioproject/755678
data_NCBI <- read_delim(file="data_R/PRJNA755678_AssemblyDetails.txt", 
                        comment = "#", 
                        col_names = c("Assembly" ,"Level", "WGS", "BioSample", "MAG", "Taxonomy")) %>%
  select(Assembly, BioSample, MAG)

Export the data

This is how the MAGs_summary_Oslofjord.tsv in Zenodo was made.

data_export <- rel_abundance %>%
  inner_join(checkm, by="MAG") %>%
  inner_join(data_stats, by="MAG") %>%
  inner_join(taxonomy, by="MAG") %>%
  inner_join(data_NCBI, by="MAG")

colnames(data_export)[2:9] <- as.character(metadata$ID)
#write_tsv(data_export, file = "MAGs_summary_Oslofjord.tsv")

Figure - Summary abundance

We do some data wrangling with the relative abundance and metadata for making a heatmap

#Join relative abundance with metadata
big.sum <- rel_abundance %>%
  pivot_longer(!MAG,
             names_to = "ID",
             values_to = "abundance") %>%
  inner_join(metadata, by="ID")

Heathmap

We select the top20 MAGs (by median abundance) and add some taxonomic information

#filter the top 20
xr <- big.sum %>%
  filter(MAG != "unmapped") %>%
  select(MAG, ID, abundance) %>%
  group_by(ID) %>%
  distinct(MAG, .keep_all = TRUE) %>%
  group_by(MAG)%>%
  summarise(av=median(abundance)) %>%
  arrange(desc(av)) %>%
  top_n(20) #top 20

## Selecting by av

#data for the plot
xr <-  big.sum %>%
  semi_join(xr, by = "MAG") %>% #filter
  inner_join(taxonomy, by = "MAG") %>%
  group_by(Phylum, Sample) %>%
  mutate(Species = na_if(Species, "s__")) %>%  #Use NA if unclassified
  mutate(Genus = na_if(Genus, "g__")) %>%
  mutate(Family = na_if(Family, "f__")) %>%
  mutate(Order = na_if(Order, "o__"))

#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)
#Chose the lowest taxonomic rank (chose the lowest that is not NA)
xr$S2 <- xr$Species %>% coalesce(xr$Species, xr$Genus, xr$Family, xr$Order, xr$Class)
xr$S2 <- paste(xr$Phylum, xr$S2, xr$MAG, sep = ", ")

#Order of the y-axis. By median abundance
x = tapply(xr$abundance, xr$S2, function(x) median(x))
x = sort(-x, TRUE)
xr$S2 = factor(as.character(xr$S2), levels=names(x))
xr$ID = factor(as.character(xr$ID), levels=sort(unique(xr$ID)))

And now we make the plot

ggplot(xr, aes(S2, ID)) +
    facet_grid(cols = vars(Location),
             scales = "free_x") +
  geom_tile(aes(fill=abundance), colour="white") +
  theme_bw() +
  coord_flip() +
   theme(
     legend.title = element_blank(),
     panel.border = 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),
      strip.background = element_blank(),
        strip.text = element_text(colour = "black", size=7, face="bold"),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        panel.grid = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.margin=margin(0,0,0,0)) +
  labs(y="Samples", fill="Relative abundance (%)") + 
  scale_fill_viridis_c(trans="sqrt")

#ggsave("Graphs/MAG_heathmap.svg", width = 95, height = 90, units = "mm", dpi = 600, limitsize = FALSE)

Boxplot

Find the top5 Phyla by median abundance

#filter the top5 phyla
xr <- big.sum %>%
  inner_join(taxonomy, by = "MAG") %>%
  filter(MAG != "unmapped") %>%
  group_by(ID, Phylum) %>%
  summarise(av=sum(abundance)) %>%
  group_by(Phylum)%>%
  summarise(av=median(av)) %>%
  arrange(desc(av)) %>%
  top_n(5)

## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
## Selecting by av

#data for the plot
#Abundance of the top MAGs
xr <-  big.sum %>%
  inner_join(taxonomy, by = "MAG") %>%
  semi_join(xr, by = "Phylum") %>% #filter
  group_by(ID, Phylum) %>%
  summarise(av=sum(abundance))

## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.