Processing metabarcoding results

1 Metabarcoding

See scripts/ for QIIME2 commands used to process raw sequences.

Below code is for analyzing all of the sequence data and compiling it for analysis in R.

1.1 Make manifest file

# | eval: false
library(tidyverse)

Make a list of files document so we can create a manifest file and import. I need to check total number of sequences per sample and the overall quality.

# | eval: false
seq_list <- read_delim(file = "scripts/list_of_files.txt", delim = "\t",col_names = F)
length(seq_list$X1)
length(unique(seq_list$X1))
fullpath <- "/scratch/group/hu-lab/data/tag-seq-data/GoM-2023-18S_2025-05/"

Parse file names for sample IDs.

# | eval: false
parsed <- seq_list %>%
  separate(X1, into = c("Order", "gom", "stn", "niskin", "template", "DNA", "REP", "RUN", "suffix", "l", "r"), remove = FALSE) %>%
  mutate(`sample-id` = case_when(
    RUN == "orig" ~ paste(gom, stn, niskin, template, REP, RUN, sep = "_"),
    (RUN != "orig") ~ paste(gom, stn, niskin, template, REP, sep = "_"))) %>%
  mutate(READ = case_when(
    (r == "R1" | r == "R2") ~ r,
    (l == "R1" | l == "R2") ~ l,
    (suffix == "R1" | suffix == "R2") ~ suffix
  )) %>%
  select(`sample-id`, READ, X1)  %>%
  pivot_wider(names_from = READ, values_from = X1)

# | eval: false
manifest_gom23 <-parsed %>%
  mutate(`forward-absolute-filepath` = paste(fullpath, R1, sep = ""),
         `reverse-absolute-filepath` = paste(fullpath, R2, sep = "")) %>%
  select(`sample-id`, `forward-absolute-filepath`, `reverse-absolute-filepath`)

Write output as a manifest file

# | eval: false
write.table(manifest_gom23, file = "manifest-gom-2023", quote=FALSE,col.names=TRUE,row.names=FALSE,sep="\t")

1.2 Compile QIIME2 output files

Import ASV table

# | eval: false
asvs <- read_delim(file = "input-data/gom-2023-18s-asv-table.tsv", delim = "\t", skip =1 ) %>% 
  select(FeatureID = `#OTU ID`, starts_with("GOM"))
# head(asvs)

Import taxonomy.

# | eval: false
tax <- read_delim("input-data/taxonomy.tsv", delim = "\t") %>% 
  select(FeatureID = `Feature ID`, Taxon)
# head(tax)

# | eval: false
asv_wtax_GoM23_062025 <- asvs %>%
  left_join(tax)

2 Sequence QC

Set up R & import starting data

# | eval: false
library(tidyverse);library(phyloseq); library(decontam)
library(compositions); library(patchwork); 
library(ggupset); library(gt)
library(plotly); library(viridis); library(vegan)

# | eval: false
seq_info <- read.csv("input-data/seq-metadata-dict.csv")
# head(seq_info)

2.1 Sequence quality control

Import to Phyloseq.

# | eval: false
tax_mat <- asv_wtax_GoM23_062025 %>% 
  select(FeatureID, Taxon) %>% 
  separate(Taxon, c("Domain", "Supergroup", 
                  "Division", "Subdivision","Class", "Order",
                  "Family", "Genus", "Species"), sep = ";", remove = FALSE) %>% 
  column_to_rownames(var = "FeatureID") %>% 
  as.matrix

asv_mat <- asv_wtax_GoM23_062025 %>% 
  select(FeatureID, starts_with("GOM")) %>% 
  column_to_rownames(var = "FeatureID") %>% 
  as.matrix

rownames(tax_mat) <- row.names(asv_mat)

2.2 Phyloseq integration

# | eval: false
ASV = otu_table(asv_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
gom_phylo = phyloseq(ASV, TAX)

# | eval: false
# sample_names(ASV)
seq_info <- seq_info %>% 
  column_to_rownames(var = "SAMPLES")

samplenames <- sample_data(seq_info)
gom_phylo_sample <- merge_phyloseq(gom_phylo, samplenames)
gom_phylo_sample

2.3 Decontam

# | eval: false
sample_data(gom_phylo_sample)$is.neg <- sample_data(gom_phylo_sample)$Sample_Ctrl  == "Control"

When “Control” appears in “Sample_or_Control column, this is a negative control” > 0.5 - this threshold will ID contaminants in all samples that are more prevalent in negative controls than in positive samples.

# | eval: false
# ID contaminants using Prevalence information
contam_prev <- isContaminant(gom_phylo_sample, 
                               method="prevalence", 
                               neg="is.neg", 
                               threshold = 0.5, normalize = TRUE) 
# ?isContaminant()
# Report number of ASVs IDed as contaminants
table(contam_prev$contaminant)

# | eval: false
list <- filter(contam_prev, contaminant == TRUE)
list_to_rm <- as.character(row.names(list))

2.4 Remove decontam ASVs and check stats

Total number of sequences and ASVs. 48,463,174 sequences 61,813 ASVs

766 ASVs to be removed.

# | eval: false
length(unique(asv_wtax_GoM23_062025$FeatureID))
sum(asv_mat)
length(list_to_rm)

# | eval: false
cleaned <- as.data.frame(asv_mat) %>% 
  rownames_to_column(var = "FeatureID") %>% 
  filter(!(FeatureID %in% list_to_rm)) %>% 
  column_to_rownames(var = "FeatureID") %>% 
  as.matrix
sum(cleaned)

After decontam 42,775,994 sequences 766 ASVs

Save cleaned ASV table files

# | eval: false
# glimpse(asv_wtax_GoM23_062025)

asv_wtax_wide_062025 <- asv_wtax_GoM23_062025 %>% 
  filter(!(FeatureID %in% list_to_rm)) %>% 
  select(FeatureID, Taxon, starts_with("GOM"))

asv_wtax_long_062025 <- asv_wtax_GoM23_062025 %>% 
  filter(!(FeatureID %in% list_to_rm)) %>% 
  pivot_longer(cols = -c(FeatureID, Taxon), names_to = "SAMPLES", values_to = "SEQUENCE_COUNT") %>% 
  filter(SEQUENCE_COUNT > 0) %>% 
  separate(Taxon, into = c("Domain", "Supergroup", "Division", "Subdivision", "Class", "Order", "Family", "Genus", "Species"), sep = ";", remove = FALSE)

# length(unique(asv_wtax_wide_062025$FeatureID))

2.5 Get sequence stats

# | eval: false
# unique(asv_wtax_long_062025$SAMPLES)
head(asv_wtax_long_062025)
asv_wtax_long_062025 %>% 
  group_by(SAMPLES) %>% 
    summarize(TOTAL_SEQ = sum(SEQUENCE_COUNT),
              TOTAL_ASV = n()) %>% 
  ggplot(aes(x = TOTAL_SEQ, y = TOTAL_ASV)) + 
  geom_point(shape = 21, color = "black", fill = "grey20") +
  theme_classic() +
  theme(axis.text = element_text(color = "black"),
      panel.grid.major = element_line(color = "grey90")) +
  labs(x = "Total sequences", y = "Total ASVs", title = "Distribtion of sequences and ASVs for all samples")

# | eval: false
# unique(asv_wtax_long_062025$SAMPLES)
# head(asv_wtax_long_062025)
asv_wtax_long_062025 %>% 
  group_by(SAMPLES, Domain) %>% 
    summarize(TOTAL_SEQ = sum(SEQUENCE_COUNT),
              TOTAL_ASV = n()) %>% 
  ggplot(aes(y = SAMPLES, x = TOTAL_ASV, fill = Domain)) + 
  geom_bar(stat = "identity", position = "stack", color = "black") +
  theme_classic() +
  theme(axis.text = element_text(color = "black"),
      panel.grid.major = element_line()) +
  labs(x = "Total sequences", y = "Total ASVs", title = "Distribtion of sequences and ASVs for all samples")

# | eval: false
summary(asv_wtax_long_062025$SEQUENCE_COUNT)
sample_count <- asv_wtax_long_062025 %>% 
  group_by(SAMPLES) %>% 
    summarize(TOTAL_SEQ = sum(SEQUENCE_COUNT),
              TOTAL_ASV = n()) %>% 
  # try removing some sequences
  filter(TOTAL_SEQ > 40000)

samples_too_low <- setdiff((unique(asv_wtax_long_062025$SAMPLES)), (unique(sample_count$SAMPLES)))

samples_too_low
head(asv_wtax_long_062025)
ctrl_samples <- unique((seq_info %>% filter(Sample_Ctrl == "Control"))$SAMPLES)

Remove control samples and samples with too few sequences.

# | eval: false
asv_long_cleaned <- asv_wtax_long_062025 %>% 
  # Remove controls
  filter(!SAMPLES %in% ctrl_samples) %>% 
  # remove samples with too few sequences
  filter(!SAMPLES %in% samples_too_low)

# | eval: false
length(unique(asv_long_cleaned$FeatureID))
sum(asv_long_cleaned$SEQUENCE_COUNT)

42,302,959 total sequences 60044 total ASV

# | eval: false
table_seq_stats <- asv_long_cleaned %>% 
  group_by(SAMPLES) %>% 
  summarise(total_asvs = n(),
            sequence_count = sum(SEQUENCE_COUNT)) %>% 
  separate(SAMPLES, into = c("rm", "STN", "NISKIN", "18s", "Excess"), remove = FALSE) %>% 
  # mutate(Replicate = case_when(
  #   Excess == "orig" ~ "1",
  #   TRUE ~ "2"
  # )) %>% 
  group_by(STN, NISKIN) %>% 
    summarise(Sequence_files = str_c(SAMPLES, collapse = ", "),
              Total_ASVs = sum(total_asvs),
              Total_sequences = sum(sequence_count))
  # pivot_wider(names_from = Replicate, values_from = c(total_asvs, sequence_count)) %>% 

table_seq_stats
# write.csv(table_seq_stats, file = "output-tables/seq-asv-stats.csv")

2.5.1 Save necessary files

# load("input-data/asv_wtax_qc_GoM23_102025.RData")
tmp <- asv_long_cleaned %>% 
  separate(SAMPLES, into = c("gom", "STATION", "NISKIN", "extra"), sep = "_") %>%
  mutate(Station = as.integer(str_remove(STATION, "S")),
         Niskin = as.integer(str_remove(NISKIN, "N"))) %>%
  unite(STN_NISKIN, STATION, NISKIN, sep = "_")
stn_niskin_seqs <- as.character(unique(tmp$STN_NISKIN))

# | eval: false
save(samples_too_low, asv_long_cleaned, stn_niskin_seqs, file =  "input-data/asv_wtax_qc_GoM23_102025.RData")

3 Curate metadata

Import & compile metadata

# Import asv data
load("input-data/asv_wtax_qc_GoM23_102025.RData")

seq_info <- read.csv("input-data/seq-metadata-dict.csv")

# Direct CTD data from shared GRAD cruise drive
master_metadata <- read.csv("input-data/Grad23_MasterSpreadsheet_01-21-26.xlsx - Master.csv", 
                            header = TRUE,
                            check.names = TRUE, na.strings = "--") %>% 
  select(-starts_with("X"), -Notes) %>% 
  separate_longer_delim(c(`NO...µmol.L.`, `PO...µmol.L.`, `SIL..µmol.L.`, `NO...µmol.L..1`, `NH...µmol.L.`), delim = " / ")
# View(master_metadata)

cellcounts <- read.csv("input-data/cellcounts-sept22.csv")

features <- read.csv("input-data/depth-features-02042026.csv")

manual_stn_info <- read.csv("input-data/manual_stn_classification.csv")

Add station ordering

offshore_on_shore_order <- c("S1", "S2", "S3", "S4", "S5","S9", "S8", "S7", "S11", "S12", "S14", "S15")

transect_labels <- c("transect1", "transect1", "transect1", "transect1", "transect1", "transect2", "transect2", "transect2", "transect3", "transect3", "transect3", "transect3")

Distance from shore

library(geosphere)
# Long, Lat order
outflow_lat <- 28.95
outflow_long <- 89.39
outflow <- c(outflow_long, outflow_lat)
outflow

Get sunrise and sunset - These calculations are based on equations provided by the National Oceanic & Atmospheric Administration (NOAA).

library(suntools)

sunriset(
  # Longitude, Latitude
  matrix(c(-86.990, 27.839), nrow = 1), # Add a negative to N longitude
  as.POSIXct("2023-07-30 19:47:00", tz = "US/Central"),
  direction='sunrise',
  # direction = 'sunset',
  POSIXct.out = TRUE
)
# Include added function
sunrise_set_newcol <- function(lat, lon, date, direction = c("sunrise", "sunset")) {
  lat_long <- cbind(lon, lat)
  sunriset(lat_long, date, direction = direction, POSIXct.out=TRUE)[,2]
}

# glimpse(master_metadata)
# unique(master_metadata$Station)
# names(master_metadata)

Add time of day in central time (local) and add in sunrise and sunset.

library(lubridate)

# Combine date and time in GMT.
metadata_allsamples <- master_metadata %>% 
  filter(Station != "Pete") %>% 
  filter(Station != "") %>% 
  filter(!is.na(Station)) %>%
  mutate(Station = as.integer(Station)) %>% 
  filter(Station != 63) %>%
  select(Station, Niskin, Depth = `CTD.Recorded.Depth..m.`, 
         Latitude = `NMEA.Latitude`, Longitude = `NMEA.Longitude`, 
         Pressure = `CTD.Pressure..db.`,
         Date, Time_UTC = `Time..UTC.`, 
         Temperature = `CTD.Temperature`, Salinity = `Corrected.CTD.Salinity`,
         Oxygen_CTD = `Corrected.CTD.Oxygen`, Oxygen_bottle = `Bottle.Oxygen..ml.l.`, 
         NO3 = `NO...µmol.L.`, PO4 = `PO...µmol.L.`, 
         SIL = `SIL..µmol.L.`, NO2 = `NO...µmol.L..1`, 
         NH4 = `NH...µmol.L.`,
         DIC = `Calculated..DIC..Apollo.`, pH = `Apollo.pH`,
         TA = `TA_Final`) %>%
  left_join(cellcounts %>% select(Station, Niskin, Prok_count = `prok.count`)) %>%
  # separate_longer_delim(c(NO3, PO4, SIL, NO2, NH4), delim = " / ") %>%
  pivot_longer(cols = c(Temperature, Salinity, Oxygen_CTD, NO3, PO4, SIL, NO2, NH4, Prok_count, DIC, TA, pH), names_to = "ENV_VARIABLE", values_to = "VALUE", values_drop_na = TRUE, values_transform = as.numeric) %>%
  group_by(Station, Niskin, Depth, Latitude, Longitude, Pressure, Date, Time_UTC, ENV_VARIABLE) %>%
    summarise(value = mean(VALUE)) %>% 
  ungroup() %>% 
  mutate(DATE_GMT = mdy(Date),
         TIME_GMT = format(strptime(Time_UTC, "%H:%M:%S"), "%H:%M:%S")) %>%
  # Use glue to combine
  mutate(datetime_string = str_glue("{DATE_GMT} {TIME_GMT}"),
    # Parse the combined string using ymd_hms() keep in GMT time.
    datetime = ymd_hms(datetime_string, tz = "GMT")) %>% 
  # Create new column with a new time zone
  mutate(datetime_cst = with_tz(datetime, tzone = "US/Central")) %>% 
  # Add sunrise
  mutate(SUNRISE = sunrise_set_newcol(Latitude, -Longitude, 
                                      datetime_cst, direction = "sunrise"),
         SUNSET = sunrise_set_newcol(Latitude, -Longitude,
                                      datetime_cst, direction = "sunset"),
         DAYTIME_INT = interval(ymd_hms(SUNRISE), ymd_hms(SUNSET)),
         DAY_NIGHT = case_when(
           ymd_hms(datetime_cst) %within% DAYTIME_INT ~ "Day",
           TRUE ~ "Night"
         )) %>% 
  mutate(HR_OF_DAY = hour(datetime_cst)) %>% 
  select(-DAYTIME_INT, -datetime_string, -datetime) %>% 
  mutate(stn = paste("S", Station, sep = "")) %>% 
  mutate(STN_ORDER = factor(stn, levels = offshore_on_shore_order),
         TRANSECT = factor(stn, levels = offshore_on_shore_order, labels = transect_labels),
         DIST_OUTFLOW = distHaversine(cbind(Longitude, Latitude), cbind(outflow_long, outflow_lat))) %>%
  mutate(DEPTH_BIN = case_when(
    Depth <= 5 ~ "< 5m",
    Depth > 5 & Depth <=25 ~ "5-25m",
    Depth > 25 & Depth <=60 ~ "25-60m",
    Depth > 60 & Depth <=270 ~ "60-270m",
    Depth > 270 & Depth <=400 ~ "270-400m",
    Depth > 400 & Depth <=800 ~ "400-800m",
    Depth > 800 & Depth <=1200 ~ "800-1200m",
    Depth > 1200 & Depth <=2000 ~ "1200-2000m",
    Depth > 2000 & Depth <=3000 ~ "2000-3000m",
  ))

Where HR_OF_DAY 19 is at sunset and 6 is at sunrise for this study.

# head(seq_info)
# stn_niskin_seqs

metadata_seqsamples <- metadata_allsamples %>% 
  mutate(STN_NISKIN = paste("S", Station, "_N", Niskin, sep = "")) %>% 
  filter(STN_NISKIN %in% stn_niskin_seqs) %>%
  filter(STN_NISKIN != "S12_N10")


metadata_allsamples <- metadata_allsamples %>% 
  mutate(STN_NISKIN = paste("S", Station, "_N", Niskin, sep = "")) %>% 
  filter(STN_NISKIN != "S12_N10")

Station information.

metadata_stn_info <- features %>% 
  select(Station, Date, Latitude, Longitude, Transect, OilRig_Count,
         MLD = MLD_BV, SAL_MAX = Salinity.maximum.depth..m., O2_MIN = Oxygen.minimum.depth..m., DCM = DCM.depth..m.,
         Depth_range = Depth_range..m., Surface_temp = Surface_temp..C, Sample_n) %>% 
  left_join(manual_stn_info)

4 Curate ASV data

# | eval: false
load(file = "input-data/asv_wtax_qc_GoM23_102025.RData", verbose = TRUE)

# head(asv_long_cleaned)
asv_long_cleaned_mod <- asv_long_cleaned %>%
  separate(SAMPLES, into = c("gom", "STATION", "NISKIN", "extra"), sep = "_") %>%
  mutate(Station = as.integer(str_remove(STATION, "S")),
         Niskin = as.integer(str_remove(NISKIN, "N"))) %>%
  unite(STN_NISKIN, STATION, NISKIN, sep = "_") %>%
  filter(STN_NISKIN != "S12_N10") %>%
  filter(STN_NISKIN != "S0_N0") %>%
  select(-gom, -extra) %>% 
  mutate(stn = paste("S", Station, sep = "")) 
  
head(asv_long_cleaned_mod)

4.1 Depth profiles across transects

head(metadata_stn_info)

metadata_stn_info %>% 
  select(Station, Transect, MLD, SAL_MAX, O2_MIN, DCM) %>% 
  pivot_longer(cols = c(MLD, SAL_MAX, O2_MIN, DCM), names_to = "FEATURE", values_to = "Depth") %>% 
  distinct() %>% 
  mutate(STATION = factor(Station, levels = c(1, 2, 3, 4, 5, 9, 8, 7, 11, 12, 14, 15)),
         DEPTH_FEATURES = factor(FEATURE, levels = c("MLD", "DCM", "SAL_MAX", "O2_MIN"), labels = c("MLD", "DCM", "Salinity max", "Oxygen min"))) %>% 
  ggplot(aes(x = as.factor(STATION), y = Depth, group = DEPTH_FEATURES,
             fill = DEPTH_FEATURES, shape = DEPTH_FEATURES, color = DEPTH_FEATURES)) +
  geom_line() +
  geom_point(size = 3, color = "white") +
  scale_color_manual(values = c("#18545e", "#c49235","#b54136", "#603655")) +
  scale_fill_manual(values = c("#18545e", "#c49235","#b54136", "#603655")) +
  scale_shape_manual(values = c(21, 23, 22, 24)) +
  scale_y_reverse() +
  facet_grid(cols = vars(Transect), scales = "free", space = "free") +
  theme_classic() +
    theme(axis.text = element_text(size = 10, colour="black"),
          panel.border = element_rect(fill = NA),
          legend.background = element_blank(),
          panel.grid.major = element_line(color = "grey90", linetype = "dotdash"),
        legend.text = element_text(size = 10, colour="black", face = "bold"),
        legend.title = element_blank(),
        strip.background = element_blank()) +
  labs(x = "Station", y = "Depth (m)")
# ggsave("figures/depth_profiles.svg", width = 8, height = 3, device = "svg", limitsize = FALSE)

5 Curate taxonomic assignment

unique(asv_long_cleaned_mod$Domain) # take EUKS
unique(asv_long_cleaned_mod$Supergroup)
unique(asv_long_cleaned_mod$Division)
# View(asv_long_cleaned_mod %>% filter(Division == "Opisthokonta"))
# View(asv_long_cleaned_mod %>%
       # select(Domain, Supergroup, Division, Subdivision, Class, Order, Family, Genus, Species) %>% distinct())

Isolating Domain to “Eukaryota”, removing bacteria, “Unassigned” and eukaryotes with :nucl. At the Division level, removing “Opisthokonta”

asv_long_wtax <- asv_long_cleaned_mod %>%
  filter(Domain == "Eukaryota") %>%
  filter(Supergroup != "nucl") %>%
  mutate(DIVISION = case_when(
    Subdivision == "X" ~ "Unannotated",
    is.na(Subdivision) ~ "Unannotated",
    TRUE ~ Subdivision
  )) %>%
  # Goal is to combine Subdivision-Class
  mutate(SUPERGROUP_CLASS = case_when(
    ## High level NAs at the supergroup and division level
    (is.na(Division) & !is.na(Supergroup)) ~ paste(Supergroup, "Unannotated", sep = "-"),
    (is.na(Division) & is.na(Supergroup)) ~ "Eukaryote-Unannotated",
    (is.na(Subdivision)) ~ paste(Supergroup, "Unannotated", sep = "-"),
    (is.na(Class) & !is.na(Subdivision)) ~ paste(Subdivision, "Unannotated", sep = "-"),
    #
    (Subdivision == "X") ~ paste(Division, "Unannotated", sep = "-"),
    (Class == "X") ~ paste(Subdivision, "Unannotated", sep = "-"),
    TRUE ~ paste(Subdivision, Class, sep = "-"))) %>%
  mutate(Supergroup_simplified = case_when(
    Supergroup == "TSAR" ~ paste("TSAR", Subdivision, sep = "-"),
    Supergroup == "Obazoa" ~ paste("Obazoa", DIVISION, sep = "-"),
    TRUE ~ Supergroup
  ))

# head(asv_long_wtax)
tax_key <- asv_long_wtax %>%
  select(FeatureID, Taxon, Supergroup, Division, Subdivision, Class, Order, Family, Genus, Species) %>% distinct()

5.0.1 Stats on what was removed

head(asv_long_cleaned_mod)
head(asv_long_wtax)
sum(asv_long_cleaned_mod$SEQUENCE_COUNT)
sum(asv_long_wtax$SEQUENCE_COUNT)
38481175/(41695184 + 38481175)
# 41,695,184
# 38,481,175
length(unique(asv_long_cleaned_mod$FeatureID))
length(unique(asv_long_wtax$FeatureID))

5.1 Average across replicates & apply ASV threshold

head(asv_long_wtax)

asv_long_avg_wtax <- asv_long_wtax %>%
#   # First average ASV sequence count across replicates
  group_by(FeatureID, Taxon, 
           DIVISION, SUPERGROUP_CLASS, Supergroup_simplified, 
           Supergroup, Division, Subdivision, Class, Order, Family, Genus, Species,
           STN_NISKIN) %>%
  summarise(MEAN_REPS_seq = mean(SEQUENCE_COUNT)) %>%
  ungroup()
length(unique(asv_long_avg_wtax$FeatureID))
# sum(asv_long_avg_wtax$MEAN_REPS_seq)
# unique(asv_long_avg_wtax$STN_NISKIN)

5.1.1 Remove ASVs with too low sequences

asvs_total_seq <- asv_long_avg_wtax %>% 
  group_by(FeatureID) %>% 
  summarise(SEQ_SUM = sum(MEAN_REPS_seq)) %>% 
  filter(SEQ_SUM >= 100)
glimpse(asvs_total_seq)
length(unique(asvs_total_seq$FeatureID))
sum(asvs_total_seq$SEQ_SUM)

# Starting with 22 mill
# 38,262 ASVs

asvs_to_keep <- unique(asvs_total_seq$FeatureID)

Cleaned and averaged: 22 million sequences and 38,262 ASVs

Filtered the whole dataset by ASVs that have more than 100 sequences total Leaving us with 21,865,501 total sequences. 10,563 ASVs

5.1.2 Rename data R objects

asv_long_avg_wtax_clean <- asv_long_avg_wtax %>% 
  filter(FeatureID %in% asvs_to_keep)

asv_long_wtax_clean <- asv_long_wtax %>% 
  filter(FeatureID %in% asvs_to_keep)

6 Save

Core data frames to take to additional code. * asv_long_avg_wtax_clean * asv_long_wtax_clean * tax_key * metadata_allsamples * metadata_seqsamples * metadata_stn_info

# | eval: false
# save(asv_long_avg_wtax_clean,
#      asv_long_wtax_clean, tax_key,
#      metadata_allsamples, metadata_seqsamples, metadata_stn_info,
#      file = "input-data/gulf-2023-sequence-analysis_02062026.RData")

7 Session Information

sessionInfo()