MICRO_zymo_072726_seq_edit
K_Langlois
2026-08-06
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)
library(ggrepel)
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"))
#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]
combo.tax<-read.delim("/Users/kylielanglois/OneDrive - SCCWRP/SD I-O data/all SD IO combo 2024/SDIO_fastq_combo_silva_tax.tsv")
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", "LOC_EVENT_ID"),
colnames(map.24))]
colnames(map.24.concise)<-gsub("LOC_EVENT_ID", "SCCWRP_ID", colnames(map.24.concise))
#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)
#first change blanks to NA
combo.tax<-apply(combo.tax, 2, function(x) gsub("^$|^ $", NA, x))
#change back to data frame
combo.tax<-as.data.frame(combo.tax)
#get last non-blank cell (exclude last taxa column)
combo.tax <- data.frame(lapply(combo.tax, as.character), stringsAsFactors=FALSE)
combo.tax$last_tax <- ifelse(!is.na(combo.tax$Taxon6),
combo.tax$Taxon6,
ifelse(!is.na(combo.tax$Taxon5),
combo.tax$Taxon5,
ifelse(!is.na(combo.tax$Taxon4),
combo.tax$Taxon4,
ifelse(!is.na(combo.tax$Taxon3),
combo.tax$Taxon3,
ifelse(!is.na(combo.tax$Taxon2),
combo.tax$Taxon2,
combo.tax$Taxon1)))))
#first change blanks to NA
tax.nb<-apply(tax.nb, 2, function(x) gsub("^$|^ $", NA, x))
#change back to data frame
tax.nb<-as.data.frame(tax.nb)
#get last non-blank cell (exclude last taxa column)
tax.nb <- data.frame(lapply(tax.nb, as.character), stringsAsFactors=FALSE)
tax.nb$last_tax <- ifelse(!is.na(tax.nb$Taxon6),
tax.nb$Taxon6,
ifelse(!is.na(tax.nb$Taxon5),
tax.nb$Taxon5,
ifelse(!is.na(tax.nb$Taxon4),
tax.nb$Taxon4,
ifelse(!is.na(tax.nb$Taxon3),
tax.nb$Taxon3,
ifelse(!is.na(tax.nb$Taxon2),
tax.nb$Taxon2,
tax.nb$Taxon1)))))
dat.tax.1<-merge(dat, tax.nb, by.x="ASV", by.y="Feature.ID")
library(plyr)
combo.tab.tax<-merge(combo.tax, combo, by.x="Feature.ID", by.y="ASV", all.y=T)
combo.tab.tax.consol.g<-ddply(combo.tab.tax, "Taxon6", numcolwise(sum)) #combine at genus level
combo.tab.tax.consol.s<-ddply(combo.tab.tax, "Taxon7", numcolwise(sum)) #combine at species level
combo.tab.tax.consol.l<-ddply(combo.tab.tax, "last_tax", numcolwise(sum)) #combine at last_tax level
dat.tax.nb.consol.g<-ddply(dat.tax.nb, "Taxon6", numcolwise(sum)) #combine at genus level
dat.tax.nb.consol.s<-ddply(dat.tax.nb, "Taxon7", numcolwise(sum)) #combine at species level
dat.tax.nb.consol.l<-ddply(dat.tax.1, "last_tax", numcolwise(sum)) #combine at last_tax level
#combine at genus level
comb26.g<-merge(combo.tab.tax.consol.g, dat.tax.nb.consol.g, by="Taxon6", all=T)
comb26.g<-comb26.g[, !grepl("Confidence", colnames(comb26.g))]
comb26.g$Taxon6<-ifelse(comb26.g$Taxon6=="", "noID", comb26.g$Taxon6)
#combine at species level
comb26.s<-merge(combo.tab.tax.consol.s, dat.tax.nb.consol.s, by="Taxon7", all=T)
comb26.s<-comb26.s[, !grepl("Confidence", colnames(comb26.s))]
comb26.s$Taxon7<-ifelse(comb26.s$Taxon7=="", "noID", comb26.s$Taxon7)
#combine at last_tax level
comb26.l<-merge(combo.tab.tax.consol.l, dat.tax.nb.consol.l, by="last_tax", all=T)
comb26.l<-comb26.l[, !grepl("Confidence", colnames(comb26.l))]comb26.phy<-comb26.g #CHOOSE WHICH TABLE TO USE <<<<<<<<<<<<<
comb26.phy<-comb26.phy[, -1]
comb26.phy<-as.data.frame(sapply(comb26.phy, as.numeric))
rownames(comb26.phy)<-comb26.l$Taxon6 #MUST MATCH WHICH TABLE YOU CHOSE <<<<<<<<<<<<<
#phyloseq------
colnames(comb26.phy)<-trimws(colnames(comb26.phy)) #double check this ALWAYS
comb26.phy[is.na(comb26.phy)]<-0 #MUST DO
OTU=otu_table(comb26.phy, taxa_are_rows = T)
map1<-rbind(map.24.concise, map.zy)
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 = "NMDS", distance = "bray", trymax=1000, na.rm=T)FALSE Square root transformation
FALSE Wisconsin double standardization
FALSE Run 0 stress 0.1151946
FALSE Run 1 stress 0.1151944
FALSE ... New best solution
FALSE ... Procrustes: rmse 0.0003151682 max resid 0.00303019
FALSE ... Similar to previous best
FALSE Run 2 stress 0.1197269
FALSE Run 3 stress 0.1431903
FALSE Run 4 stress 0.1197267
FALSE Run 5 stress 0.1285245
FALSE Run 6 stress 0.1151944
FALSE ... Procrustes: rmse 5.108233e-05 max resid 0.0004576405
FALSE ... Similar to previous best
FALSE Run 7 stress 0.1436037
FALSE Run 8 stress 0.1151944
FALSE ... New best solution
FALSE ... Procrustes: rmse 2.262675e-05 max resid 0.0002036514
FALSE ... Similar to previous best
FALSE Run 9 stress 0.1151944
FALSE ... New best solution
FALSE ... Procrustes: rmse 0.000298305 max resid 0.002926746
FALSE ... Similar to previous best
FALSE Run 10 stress 0.1386458
FALSE Run 11 stress 0.413419
FALSE Run 12 stress 0.1197267
FALSE Run 13 stress 0.1151941
FALSE ... New best solution
FALSE ... Procrustes: rmse 0.0001815373 max resid 0.001810974
FALSE ... Similar to previous best
FALSE Run 14 stress 0.1254591
FALSE Run 15 stress 0.1343593
FALSE Run 16 stress 0.1372859
FALSE Run 17 stress 0.1151942
FALSE ... Procrustes: rmse 8.332262e-05 max resid 0.000805532
FALSE ... Similar to previous best
FALSE Run 18 stress 0.143374
FALSE Run 19 stress 0.1283221
FALSE Run 20 stress 0.1474628
FALSE *** Solution reachedcoords<-as.data.frame(n.o$points) #"points" for NMD "vectors" for PCoA
coords$SEQUENCING_NAME<-rownames(coords)
coords2<-coords
coords2<-merge(coords2, map1, by="SEQUENCING_NAME")