library(knitr)
library(readxl)
#detach(package:plyr)
library(dplyr)
library(ggplot2)
library(reshape2)
library(kableExtra)
library(ggh4x)
library(stringr)
library(seqinr)
library(palettetown)
library(phyloseq)
set.seed(100)
path<-
"/Users/kylielanglois/SCCWRP/KL data - General/sequencing/data/mASCI-FHAB_071526"

d<-read.csv(file.path(path, 
    "BIO_mASCI_rcbL_071526/fastqc_multiqc_dada2_filteronly_results_V2 2.csv"))

f<-read.csv(file.path(path, 
"BGC_FHAB_16Scyano_071326/BGC_FHAB_16Scyano_071326_fastqc_multiqc_dada2_filteronly_results.csv"))

v<-read.delim(file.path(path, 
"BIO_mASCI_rcbL_071526/BIO_mASCI_rcbL_071526_trimmed_no_adpt_demux_Vmerged_derep_denovo_table.tsv"), 
 header = F)

no.amp<-read_excel("/Users/kylielanglois/SCCWRP/eDNA Sample Management - eDNA/mASCI/Sequencing plate maps/SCCWRP - mASCI-FHAB_May2026-CGASH251017-5-PCR.xlsx", sheet = 1)
no.amp.rcbl<-subset(no.amp, no.amp$`rcbL band?`=="n")
no.amp.samp.rcbl<-no.amp.rcbl$`Sample Name`
no.amp.cyano<-subset(no.amp, no.amp$`16Scyano band?`=="n")
no.amp.samp.cyano<-no.amp.cyano$`Sample Name`

rcbL

compare vsearch to dada2

#plot sample types and PCR amplification results---------------
t.m<-melt(d.pick, id.vars = "sample.names", 
          measure.vars = c("reads.out", "no.chimera"))
t.m$value<-as.numeric(t.m$value)
t.m$sample.names<-sub(".*-rcbl_", "", t.m$sample.names)
t.m$SampleType<-ifelse(grepl("EB_", t.m$sample.names), "EB", 
                       ifelse(grepl("_FB", t.m$sample.names), "FB", 
                              ifelse(grepl("DUP", t.m$sample.names), "DUP", "sample")))

pg<-ggplot(data = t.m)+
  geom_bar(aes(x=sample.names, y=value, 
               fill=variable, color=SampleType), 
           stat = "identity", position = position_dodge2())+
  annotate(geom = "text", x=no.amp.samp.rcbl, y=0, label="*", 
           vjust=1, size=6)+
  scale_fill_manual(values=c( "grey", "green"), 
                    labels=c("filtered", "non.chimeric"))+
  scale_color_manual(values=c("darkblue",  "darkred", "purple", "darkblue"))+
  theme_bw()+
  theme(axis.text.x = element_text(angle=45, hjust=1, size=6))
pg  

d.v.compare<-merge(d.pick, v.pick, by="sample.names")
d.v.compare$perc.retain.v<-
  round((as.numeric(d.v.compare$v.out)/as.numeric(d.v.compare$reads.in)), 4)*100 
d.v.m<-melt(d.v.compare, id.vars = "sample.names", 
            measure.vars = c("perc.retain", "perc.retain.v"))
d.v.m$blank<-ifelse(grepl("FB|EB", d.v.m$sample.names), "blank", "sample")

tg<-ggplot(data=d.v.m)+
  geom_bar(aes(x=sample.names, y=value, fill=variable), 
           stat = "identity", position = position_dodge2())+
  scale_fill_manual(values=c("#f1a340", "#542788"), 
                    labels = c("dada2", "vsearch"), 
                    name="")+
  facet_grid(~blank, scales = "free_x", space = "free_x")+
  labs(x="sample", y="sequences retained (%)",
       title = "mASCI rcbL 7-15-26 sequence processing")+
  theme_bw()+
  theme(axis.text.x = element_blank(), 
        legend.position = "bottom")
tg

edit dada2 processed file

dat<-read.csv(file.path(
  path, "BIO_mASCI_rcbL_071526/BIO_mASCI_rcbL_071526_trimmed_noadpt_V2_7_table.csv"))

tax<-read.delim(file.path(
  path, 
  "BIO_mASCI_rcbL_071526/BIO_mASCI_rcbL_071526_trimmed_noadpt_V2_7_diat.rcbL.v12_bl_tax_99.tsv"))
tax<-tax[, !grepl("Consens|num", colnames(tax))]

map<-read.csv(
  "/Users/kylielanglois/OneDrive - SCCWRP/algae/mASCI/eDNAmaster_inventory_simple_2026-07-17.csv")
dat.p<-dat
rownames(dat.p)<-dat.p$ASV #make rows into taxa
dat.p<-dat.p[, -1]
OTU = otu_table(dat.p, taxa_are_rows = T)

crosswalk.rcbl<-data.frame(SEQ_ID=colnames(dat.p)) #get sequencer name
crosswalk.rcbl$SCCWRP_ID<-sub(".*\\.rcbl_", "", crosswalk.rcbl$SEQ_ID)
#get SCCWRP_ID from sequencing name
map1<-map[map$SCCWRP_ID %in% crosswalk.rcbl$SCCWRP_ID, ] 
#only some samples from complete metadata of algae
map1<-merge(map1, crosswalk.rcbl, by="SCCWRP_ID") #add to metadata
map1$no.amp.rcbl<-ifelse(map1$SCCWRP_ID %in% no.amp.rcbl$`Sample Name`, 
                         "no_amp", "")
#add PCR amplification to metadata
row.names(map1)<-map1$SEQ_ID
samp.no.fb<-map1$SEQ_ID[map1$Replicate.!="FB"]
MAP = sample_data(map1)

physeq=phyloseq(OTU, MAP) #rownmaes of map must match colnames of otu
physeq1<-prune_samples(samp.no.fb, physeq)

n.o<-ordinate(physeq1, method = "NMDS", distance = "bray", trymax=1000) 
FALSE Square root transformation
FALSE Wisconsin double standardization
FALSE Run 0 stress 0.1737745 
FALSE Run 1 stress 0.1863245 
FALSE Run 2 stress 0.1781276 
FALSE Run 3 stress 0.1976341 
FALSE Run 4 stress 0.1991994 
FALSE Run 5 stress 0.1832061 
FALSE Run 6 stress 0.1789341 
FALSE Run 7 stress 0.1870938 
FALSE Run 8 stress 0.2070798 
FALSE Run 9 stress 0.1855384 
FALSE Run 10 stress 0.2133659 
FALSE Run 11 stress 0.1895204 
FALSE Run 12 stress 0.1841769 
FALSE Run 13 stress 0.1895609 
FALSE Run 14 stress 0.173128 
FALSE ... New best solution
FALSE ... Procrustes: rmse 0.02578621  max resid 0.1842313 
FALSE Run 15 stress 0.1870165 
FALSE Run 16 stress 0.1853513 
FALSE Run 17 stress 0.1753686 
FALSE Run 18 stress 0.1965167 
FALSE Run 19 stress 0.2213667 
FALSE Run 20 stress 0.2147681 
FALSE Run 21 stress 0.1856002 
FALSE Run 22 stress 0.187203 
FALSE Run 23 stress 0.1840341 
FALSE Run 24 stress 0.2012099 
FALSE Run 25 stress 0.1980938 
FALSE Run 26 stress 0.202871 
FALSE Run 27 stress 0.1831125 
FALSE Run 28 stress 0.1948417 
FALSE Run 29 stress 0.2140616 
FALSE Run 30 stress 0.1953493 
FALSE Run 31 stress 0.1990194 
FALSE Run 32 stress 0.1751878 
FALSE Run 33 stress 0.1972706 
FALSE Run 34 stress 0.1928177 
FALSE Run 35 stress 0.1775912 
FALSE Run 36 stress 0.2130295 
FALSE Run 37 stress 0.1905586 
FALSE Run 38 stress 0.1956952 
FALSE Run 39 stress 0.197765 
FALSE Run 40 stress 0.186018 
FALSE Run 41 stress 0.2142366 
FALSE Run 42 stress 0.173429 
FALSE ... Procrustes: rmse 0.02127945  max resid 0.1827075 
FALSE Run 43 stress 0.1955851 
FALSE Run 44 stress 0.1737506 
FALSE Run 45 stress 0.1772925 
FALSE Run 46 stress 0.1751539 
FALSE Run 47 stress 0.193521 
FALSE Run 48 stress 0.1868334 
FALSE Run 49 stress 0.2020367 
FALSE Run 50 stress 0.1860774 
FALSE Run 51 stress 0.1777607 
FALSE Run 52 stress 0.2004483 
FALSE Run 53 stress 0.1926412 
FALSE Run 54 stress 0.1852696 
FALSE Run 55 stress 0.1896 
FALSE Run 56 stress 0.1869707 
FALSE Run 57 stress 0.1819038 
FALSE Run 58 stress 0.1731276 
FALSE ... New best solution
FALSE ... Procrustes: rmse 0.00100322  max resid 0.006409925 
FALSE ... Similar to previous best
FALSE *** Solution reached
#try "NMDS" and "PCoA"

## [1] "# of ASV: 3148"
## [1] "total # of seqs: 2409743"

### 16Scyano compare AL vs KL analyze

dat.c<-read.csv(file.path(
  path, "BGC_FHAB_16Scyano_071326/BGC_FHAB_16Scyano_071326_V2_2_table.csv"))

tax.c<-read.delim(file.path(
  path, 
  "BGC_FHAB_16Scyano_071326/BGC_FHAB_16Scyano_071326_Monchamp_bl_tax.tsv"))
tax.c<-tax.c[, !grepl("Consens|num", colnames(tax.c))]

dat.a<-read.delim(file.path(
  path, "BGC_FHAB_16Scyano_071326/BGC_FHAB_16Scyano_071326_table_AL.tsv"), 
  header = F)
colnames(dat.a)<-dat.a[2, ]
dat.a<-dat.a[3:nrow(dat.a), ]

tax.a<-read.delim(file.path(
  path, 
  "BGC_FHAB_16Scyano_071326/BGC_FHAB_16Scyano_071326_rep-set-seqs_AL_Monchamp_bl_tax.tsv"))
tax.a<-tax.a[, !grepl("Consens|num", colnames(tax.a))]


map<-read.csv(
  "/Users/kylielanglois/OneDrive - SCCWRP/algae/mASCI/eDNAmaster_inventory_simple_2026-07-17.csv")
## [1] "# of ASV: 36694"
## [1] "total # of seqs: 1922601"

## [1] "# of ASV: 34542"
## [1] "total # of seqs: 1857662"