# Clear the R workspace rm(list = ls()) # Optionally, clear the console as well cat("\014") # This will clear the console in most R environments # Install required packages if not already installed if (!require("readr")) install.packages("readr") if (!require("dplyr")) install.packages("dplyr") if (!require("tidyr")) install.packages("tidyr") if (!require("ggplot2")) install.packages("ggplot2") # Load necessary libraries library(readr) library(dplyr) library(tidyr) library(ggplot2) # Load the ASV table asv_table <- read.csv(file='asv_table.csv') # Remove non-numeric columns (e.g., identifiers) for sequencing depth calculation numeric_asv_table <- asv_table %>% select(where(is.numeric)) # Calculate column sums cs <- colSums(numeric_asv_table) # Compute mean and median sequencing depth mean_depth <- mean(cs) median_depth <- median(cs) # Display results mean_depth median_depth # Rename the first column to ASV_ID colnames(asv_table)[1] <- "ASV_ID" # Define the columns of interest columns_of_interest <- c( "NIOZ153.005.006", "NIOZ153.007.008", "NIOZ153.009.010", "NIOZ153.011.012", "NIOZ153.013.014", "NIOZ153.015.016", "NIOZ153.019.020", "NIOZ153.021.022", "NIOZ153.023.024", "NIOZ153.025.026", "NIOZ153.029.030", "NIOZ153.031.032", "NIOZ153.033.034", "NIOZ153.035.036", "NIOZ153.037.038", "NIOZ153.039.040", "NIOZ153.041.042", "NIOZ153.043.044", "NIOZ153.045.046", "NIOZ153.047.048", "NIOZ153.049.050", "NIOZ153.051.052", "NIOZ153.053.054", "NIOZ153.055.056", "NIOZ153.057.058", "NIOZ153.059.060", "NIOZ153.061.062", "NIOZ309.001.250", "NIOZ309.003.252", "NIOZ309.005.254", "NIOZ309.007.256", "NIOZ309.009.002", "NIOZ309.011.004", "NIOZ309.013.006", "NIOZ309.015.008", "NIOZ309.017.010", "NIOZ309.019.012", "NIOZ309.021.014", "NIOZ309.023.016", "NIOZ309.025.018", "NIOZ309.027.020", "NIOZ309.029.022", "NIOZ309.031.024", "NIOZ309.033.026", "NIOZ309.035.028", "NIOZ309.037.030", "NIOZ309.039.032", "NIOZ309.041.034", "NIOZ309.043.036", "NIOZ309.049.042", "NIOZ309.051.044", "NIOZ309.053.046", "NIOZ309.055.048", "NIOZ153.065.066", "NIOZ153.069.070", "NIOZ153.073.074", "NIOZ153.075.076", "NIOZ153.077.078", "NIOZ153.081.082", "NIOZ153.085.086", "NIOZ153.089.090", "NIOZ153.093.094", "NIOZ153.097.098", "NIOZ153.101.102", "NIOZ153.105.106", "NIOZ153.191.192", "NIOZ153.113.114", "NIOZ153.067.068", "NIOZ153.071.072", "NIOZ153.079.080", "NIOZ153.083.084", "NIOZ153.087.088", "NIOZ153.091.092", "NIOZ153.095.096", "NIOZ153.099.100", "NIOZ153.103.104", "NIOZ153.107.108", "NIOZ153.111.112", "NIOZ153.115.116", "NIOZ153.117.118", "NIOZ153.119.120", "NIOZ153.121.122", "NIOZ153.123.124", "NIOZ153.125.126", "NIOZ153.129.130", "NIOZ153.133.134", "NIOZ153.137.138", "NIOZ153.141.142", "NIOZ153.145.146", "NIOZ153.149.150", "NIOZ153.157.158", "NIOZ153.161.162", "NIOZ153.165.166", "NIOZ153.169.170", "NIOZ153.173.174", "NIOZ153.127.128", "NIOZ153.131.132", "NIOZ153.135.136", "NIOZ153.139.140", "NIOZ153.143.144", "NIOZ153.151.152", "NIOZ153.159.160", "NIOZ153.163.164", "NIOZ153.167.168", "NIOZ153.171.172", "NIOZ153.175.176", "NIOZ153.181.182", "NIOZ153.183.184", "NIOZ153.185.186", "NIOZ153.187.188", "NIOZ309.057.050", "NIOZ309.059.052", "NIOZ309.061.054", "NIOZ309.063.056", "NIOZ309.065.058", "NIOZ309.067.060", "NIOZ309.069.062", "NIOZ309.071.064", "NIOZ309.073.066", "NIOZ309.075.068", "NIOZ309.077.070", "NIOZ309.079.072", "NIOZ309.085.078", "NIOZ309.087.080", "NIOZ309.089.082", "NIOZ309.095.088", "NIOZ309.099.092", "NIOZ309.103.096", "NIOZ309.107.100", "NIOZ309.111.104", "NIOZ309.115.108", "NIOZ309.119.112", "NIOZ309.123.116", "NIOZ309.127.120", "NIOZ309.131.124", "NIOZ309.135.128", "NIOZ309.163.156", "NIOZ309.097.090", "NIOZ309.101.094", "NIOZ309.105.098", "NIOZ309.109.102", "NIOZ309.113.106", "NIOZ309.117.110", "NIOZ309.121.114", "NIOZ309.125.118", "NIOZ309.129.122", "NIOZ309.133.126", "NIOZ309.137.130", "NIOZ309.165.158", "NIOZ309.171.164", "NIOZ309.173.166", "NIOZ309.175.168", "NIOZ309.177.170", "NIOZ319.001.240", "NIOZ319.003.242", "NIOZ319.005.244", "NIOZ319.007.246", "NIOZ319.009.248", "NIOZ319.011.250", "NIOZ319.013.252", "NIOZ319.015.254", "NIOZ319.017.256", "NIOZ319.019.002", "NIOZ319.021.004", "NIOZ319.023.006", "NIOZ319.025.008", "NIOZ319.027.010", "NIOZ319.029.012", "NIOZ319.031.014", "NIOZ319.033.016", "NIOZ319.035.018", "NIOZ319.037.020", "NIOZ319.039.022", "NIOZ319.043.026", "NIOZ319.045.028", "NIOZ319.047.030", "NIOZ319.049.032", "NIOZ319.051.034", "NIOZ319.053.036", "NIOZ319.055.038", "NIOZ319.057.040", "NIOZ319.059.042", "NIOZ319.085.068", "NIOZ319.087.070", "NIOZ319.089.072", "NIOZ319.091.074", "NIOZ319.093.076", "NIOZ319.095.078", "NIOZ319.097.080", "NIOZ319.099.082", "NIOZ319.101.084", "NIOZ319.103.086", "NIOZ319.105.088", "NIOZ319.107.090", "NIOZ319.109.092", "NIOZ319.113.096", "NIOZ319.115.098", "NIOZ319.117.100", "NIOZ319.119.102", "NIOZ319.121.104", "NIOZ319.123.106", "NIOZ319.125.108", "NIOZ319.127.110", "NIOZ319.129.112", "NIOZ319.131.114", "NIOZ319.133.116", "NIOZ319.135.118", "NIOZ319.137.120", "NIOZ319.139.122", "NIOZ319.141.124", "NIOZ319.143.126", "NIOZ319.145.128", "NIOZ319.147.130", "NIOZ319.149.132", "NIOZ319.173.156", "NIOZ319.175.158", "NIOZ319.177.160", "NIOZ319.179.162", "NIOZ319.181.164", "NIOZ319.183.166", "NIOZ319.185.168", "NIOZ319.187.170", "NIOZ319.189.172", "NIOZ319.191.174", "NIOZ319.193.176", "NIOZ319.195.178", "NIOZ319.197.180", "NIOZ319.199.182", "NIOZ319.201.184", "NIOZ319.203.186", "NIOZ319.205.188", "NIOZ319.207.190", "NIOZ319.209.192", "NIOZ319.211.194", "NIOZ319.213.196", "NIOZ319.229.212", "NIOZ321.001.232", "NIOZ319.239.222", "NIOZ319.241.224", "NIOZ319.243.226", "NIOZ321.003.234", "NIOZ321.005.236", "NIOZ321.007.238", "NIOZ321.009.240", "NIOZ321.011.242", "NIOZ321.013.244", "NIOZ321.015.246", "NIOZ321.017.248", "NIOZ321.019.250", "NIOZ321.021.252", "NIOZ321.023.254", "NIOZ321.025.256", "NIOZ321.027.002", "NIOZ321.029.004", "NIOZ321.031.006", "NIOZ321.033.008", "NIOZ321.035.010", "NIOZ321.037.012", "NIOZ321.039.014", "NIOZ321.041.016", "NIOZ321.043.018", "NIOZ321.057.032", "NIOZ321.059.034", "NIOZ321.061.036", "NIOZ321.063.038", "NIOZ321.065.040", "NIOZ321.067.042", "NIOZ321.069.044", "NIOZ321.071.046", "NIOZ321.073.048", "NIOZ321.075.050", "NIOZ321.077.052", "NIOZ321.079.054", "NIOZ321.081.056", "NIOZ321.083.058", "NIOZ321.085.060", "NIOZ321.087.062", "NIOZ321.089.064", "NIOZ321.091.066", "NIOZ321.093.068", "NIOZ321.095.070", "NIOZ321.097.072", "NIOZ321.099.074", "NIOZ321.101.076", "NIOZ321.103.078", "NIOZ321.121.096", "NIOZ321.123.098", "NIOZ321.125.100", "NIOZ321.127.102", "NIOZ321.129.104", "NIOZ321.131.106", "NIOZ321.133.108", "NIOZ321.135.110", "NIOZ321.137.112", "NIOZ321.139.114", "NIOZ321.141.116", "NIOZ321.143.118", "NIOZ321.145.120", "NIOZ321.147.122", "NIOZ321.149.124", "NIOZ321.151.126", "NIOZ321.153.128", "NIOZ321.155.130", "NIOZ321.157.132", "NIOZ321.159.134", "NIOZ321.161.136", "NIOZ321.163.138", "NIOZ321.165.140", "NIOZ321.167.142", "NIOZ321.181.156", "NIOZ321.183.158", "NIOZ321.185.160", "NIOZ321.187.162", "NIOZ321.189.164", "NIOZ321.191.166", "NIOZ321.193.168", "NIOZ321.195.170", "NIOZ321.197.172" ) # Filter out columns_of_interest to only include those that exist in asv_table available_columns <- columns_of_interest[columns_of_interest %in% colnames(asv_table)] # Select the ASV_ID and the available columns selected_data <- asv_table %>% select(ASV_ID, all_of(available_columns)) # Display the selected data to confirm head(selected_data) # Calculate the relative abundance (RA) for each selected column in available_columns selected_data <- selected_data %>% mutate(across(all_of(available_columns), ~ . / sum(., na.rm = TRUE), .names = "{.col}_RA")) # Kick out singletons: Keep only rows where the sum of numeric columns is greater than 1 no_singleton <- selected_data %>% filter(rowSums(select(., where(is.numeric))) > 1) # Load the taxonomy data, assuming it has the ASV_ID column and phylogenetic details taxonomy_data <- read_csv(file='taxotable.csv') # Rename the first column to "ASV_ID" if name is lacking colnames(taxonomy_data)[1] <- "ASV_ID" # Join the taxonomy data with the selected data based on ASV_ID selected_data_with_phylogeny <- no_singleton %>% left_join(taxonomy_data, by = "ASV_ID") head(selected_data_with_phylogeny[, (ncol(selected_data_with_phylogeny) - 8):ncol(selected_data_with_phylogeny)]) # Define lists of MOB genera and families mob_genera <- c( "Crenothrix", "ET-SHO", "IheB2-23", "Incertae Sedis", "Marine Methylotrophic Group 2", "Methylicorpusculum", "Methylobacter", "Methyloglobulus", "Methylomarinum", "Methylomicrobium", "Methylomonas", "Methyloprofundus", "Methylosarcina", "Methylosoma", "Methylosphaera", "Methylotuvimicrobium", "Methylovulum", "Milano-WF1B-03", "Milano-WF1B-42", "pItb-vmat-59", "pLW-20", "Candidatus Methylospira", "Incertae Sedis", "Methylocaldum", "Methylococcus", "Methylogaea", "Methylomagnum", "Methyloparacoccus", "Methyloterricola", "Methylohalobius", "Methylomarinovum", "Methylothermus", "Methylosinus", "Methylocystis", "Methylocella", "Methylocapsa", "Methyloferula", "Methyloceanibacter", "Methylacidiphilum", "Methylacidimicrobium", "Candidatus Methylomirabilis" ) mob_families <- c("Methylomonadaceae", "Methylococcacea", "Methylohalobiaceae", "Methylacidiphilaceae", "Methylomirabilaceae") # Add the MOB column: 'MOB' if genus or family matches, else NA merged_data <- selected_data_with_phylogeny %>% mutate(MOB = ifelse(Genus %in% mob_genera | Family %in% mob_families, "MOB", NA)) # Replace NA in Genus with "Family-NA" if MOB is "MOB" and Family is not NA merged_data <- merged_data %>% mutate(Genus = ifelse(is.na(Genus) & MOB == "MOB" & !is.na(Family), paste0(Family, "-NA"), Genus)) # Filter out rows where there is a “MOB” entry in the MOB column # and select only columns with `_RA` suffix and essential metadata columns mob_data <- merged_data %>% filter(MOB == "MOB") %>% select(ASV_ID, Genus, Family, MOB, matches("^NIOZ.*_RA$")) # Add the 'type_mob' column to the mob_data dataframe mob_data <- mob_data %>% mutate( type_mob = case_when( Family %in% c("Methylomonadaceae", "Methylococcaceae", "Methylohalobiaceae") ~ "I", Family %in% c("Beijerinckiaceae", "Methyloligellaceae") ~ "II", Family %in% c("Methylacidiphilaceae") ~ "III", Family %in% c("Methylomirabilaceae") ~ "N-AOM", TRUE ~ NA_character_ # Assign NA for non-matching families within MOB data ) ) head(mob_data[, (ncol(mob_data) - 5):ncol(mob_data)]) # Reorder the columns to place 'type_mob' right after 'MOB' mob_data <- mob_data %>% relocate(type_mob, .after = MOB) # Display the first few rows of the modified mob_data dataframe head(mob_data) # Group by Genus and retain Family, MOB, and type_mob while summing NIOZ columns merged_mob_data <- mob_data %>% group_by(Genus, Family, MOB, type_mob) %>% summarise(across(starts_with("NIOZ"), sum, na.rm = TRUE), .groups = 'drop') # Aggregate by summing NIOZ columns # Remove rows with no entries merged_mob_data <- merged_mob_data %>% filter( rowSums(select(., where(is.numeric)), na.rm = TRUE) > 0,) # Retain rows with non-zero numerical entries # View a preview of the resulting dataset head(merged_mob_data)