set up

packages

library(knitr)
library(readxl)
#detach(package:plyr)
library(dplyr)
library(ggplot2)
library(reshape2)
library(kableExtra)
library(ggh4x)
library(stringr)
library(Biostrings)
library(palettetown)
library(phyloseq)
set.seed(100)

files

path<-
"/Users/kylielanglois/SCCWRP/KL data - General/sequencing/data/MICRO_zymo_072726"

#zymo 7/27/26 
dat<-read.delim(file.path(path, 
    "MICRO_zymo_072726_V2_4_table.csv"))
dat$ASV<-trimws(dat$ASV)
dat.seq<-read.csv(file.path(path, 
    "MICRO_zymo_072726_merger_table_V2_4.csv"))

map.zy<-read.csv(file.path(path, "MICRO_zymo_072726_map.csv"))
map.zy1<-map.zy[, match(c("SEQUENCING_NAME", "SAMPLE_TYPE"), colnames(map.zy))]

#BLAST taxonomy
tax<-read.delim(file.path(path, "MICRO_zymo_072726_silva-138-99_bl_tax.tsv"))
tax$Feature.ID<-trimws(tax$Feature.ID)
dat.tax<-merge(dat, tax,  by.x = "ASV", by.y = "Feature.ID", all.x = T)

#NB taxonomy
tax.nb<-read.delim(file.path(path, "MICRO_zymo_072726_V2_4_silva_NB_taxonomy.tsv"))
tax.nb$Feature.ID<-trimws(tax.nb$Feature.ID)
dat.tax.nb<-merge(dat, tax.nb,  by.x = "ASV", by.y = "Feature.ID", all.x = T)
#previously combined table
#combo<-read.csv("/Users/kylielanglois/OneDrive - SCCWRP/SD I-O data/all SD IO combo 2024/SDIO_fastq_combo_table_norep_nobl_nozero.csv")

#previously combined fasta
#combo.seq<-readDNAStringSet(
#  "/Users/kylielanglois/OneDrive - SCCWRP/SD I-O data/all SD IO combo 2024/SDIO_fastq_combo_rep_set.fasta")
#combo.seq.repset<-data.frame(ASV=names(combo.seq), 
#                             seqs=paste(combo.seq))
#combo.tab.seq<-merge(combo.seq.repset, combo, by="ASV", all.y=T)
#combo.tab.seq<-combo.tab.seq[, -1]

map.24<-read.csv("/Users/kylielanglois/OneDrive - SCCWRP/SD I-O data/all SD IO combo 2024/SD_IO_metadata_combo_012624_norep_nobl_nozero_latlong.csv")
map.24.concise<-map.24[, match(c("SEQUENCING_NAME", "SAMPLE_TYPE"), 
                               colnames(map.24))]

#merged in qiime
comb26<-read.delim(file.path(path, 
                             "MICRO_zymo_072726_fastq_combo24_merged.txt"), 
                   header = F)
colnames(comb26)<-comb26[2, ]
colnames(comb26)[1]<-"Feature.ID"
comb26<-comb26[3:nrow(comb26), ]

basic for Zymo 7/27/26

## [1] "total sequences: 412922"
## [1] "62.06% sequences BLAST unassigned"
## [1] "0.1% sequences NB unassigned"

compare to SD-IO combined (2024)

comb26.phy<-comb26
comb26.phy<-comb26.phy[, -1]
comb26.phy<-as.data.frame(sapply(comb26.phy, as.numeric))
rownames(comb26.phy)<-comb26$Feature.ID

#phyloseq------
colnames(comb26.phy)<-trimws(colnames(comb26.phy)) #double check this ALWAYS
OTU=otu_table(comb26.phy, taxa_are_rows = T) 

map1<-rbind(map.24.concise, map.zy1)
rownames(map1)<-map1$SEQUENCING_NAME #rownames(map) must match colnames(otu)
MAP=sample_data(map1)

physeq=phyloseq(OTU, MAP)

#remove samples types--------
samples_keep<-map1[!grepl("POLLUTOGRAPH|ENCAMP|EXPERIMENT", map1$SAMPLE_TYPE), ]
physeq1<-prune_samples(samples_keep$SEQUENCING_NAME, physeq)

#phyloseq ordination---------
n.o<-ordinate(physeq1, method = "PCoA", distance = "bray", trymax=1000)

coords<-as.data.frame(n.o$vectors)
coords$SEQUENCING_NAME<-rownames(coords)
coords2<-coords
coords2<-merge(coords2, map1, by="SEQUENCING_NAME")