Seed_analysis_sequence_decontam.Rmd

---
title: "Argonne2023_seed_analysis"
author: "Abby Sulesky-Grieb"
date: "2023-07-11"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Steps for analysis:
- Create phyloseq objects for both projects
- Decontam by extraction group
- Rarefy


# Load packages
```{r packages, results = FALSE, message=FALSE, warning=FALSE}
require(microViz)
require(corncob)
require(ggraph)
require(DT)
library(phyloseq)
library(ggplot2)
library(tidyverse)
library(decontam)
library(scales)
library(vegan)
library(microbiome)
library(dplyr)

```


# Phyloseq Objects
### Create Full Phyloseq object
```{r phyloseq}
otu <- read.csv(file="/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/otu_table_ArgonneSeeds_June2023.csv", sep=",", row.names=1)

tax <- read.csv(file="/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/taxonomy_ArgonneSeeds_June2023.csv", sep=",", row.names = 1)
# have to include this matrix line for phyloseq to build tax table properly
tax <- as.matrix(tax)

metadata <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/Multigen_seeds_metadata_NewArgonne.csv", sep=",", row.names=1) 

OTU = otu_table(otu, taxa_are_rows = TRUE)
TAX = tax_table(tax)
meta = sample_data(metadata)

# check outputs
#tax_table(TAX)
#otu_table(OTU)
#sample_data(meta)

full_phyloseq <- merge_phyloseq(OTU, meta, TAX)
full_phyloseq


# Subset by project
rainout_phyloseq <- full_phyloseq %>% ps_filter(Project == "rainout", .keep_all_taxa = TRUE)
multigen_phyloseq <- full_phyloseq %>% ps_filter(Project == "multigen", .keep_all_taxa = TRUE)
rainout_phyloseq
multigen_phyloseq

# Save phyloseq objects
saveRDS(rainout_phyloseq, file="/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/rainout_phyloseq.rds", compress = TRUE)

saveRDS(multigen_phyloseq, file="/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/multigen_phyloseq.rds", compress = TRUE)
```

# Decontam the multigen experiment
### Start by removing host reads
```{r host_remove}
multigen_phyloseq <- readRDS(file="/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/Seed_transmission_Common_Bean/R_Analysis_Files/multigen_phyloseq.rds")

multigen_metadata <- sample_data(multigen_phyloseq)

sort(sample_sums(multigen_phyloseq))
sum(sample_sums(multigen_phyloseq))
# 126,813,614

# start with removing mitochondria and chloroplast
multigen_phyloseq_clean <- multigen_phyloseq %>%
  subset_taxa(
      Family   != "f__Chloroplast" &
      Family  != "f__Mitochondria" &
      Phylum != ""
  )
multigen_phyloseq_clean
sum(sample_sums(multigen_phyloseq_clean))
# 12,305,095
sort(sample_sums(multigen_phyloseq_clean))
((126813614-12305095)/126813614)*100

#  % taxa that were removed:
percent_taxa_removed <- ((9515-8276)/9515)*100
percent_taxa_removed

## inspect library sizes
# before host removal
full_multigen_df <- as.data.frame(multigen_metadata) # Put sample_data into a ggplot-friendly data.frame
full_multigen_df$LibrarySize <- sample_sums(multigen_phyloseq)
full_multigen_df <- full_multigen_df[order(full_multigen_df$LibrarySize),]
full_multigen_df$Index <- seq(nrow(full_multigen_df))
ggplot(data=full_multigen_df, aes(x=Index, y=LibrarySize, color=sample_type)) + geom_point() + ggtitle("Reads before host removal")

# after host removal
df <- as.data.frame(sample_data(multigen_phyloseq_clean)) 
df$LibrarySize <- sample_sums(multigen_phyloseq_clean)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=sample_type)) + geom_point() + ggtitle("Reads after host removal")

write.csv(otu_table(multigen_phyloseq_clean), "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/ASVtable_hostremoved_predecontam.csv")

# 2000-5000 reads per sample after host removal


# decontam with seqcenter negatives  & sample_data(decontam_root_phyloseq_clean)$Compartment == "Rhizosphere"
sample_data(multigen_phyloseq_clean)$is.neg <- sample_data(multigen_phyloseq_clean)$sample_type == "seqcenter_negative"
contamdf.prev0.1_seqcenter <- isContaminant(multigen_phyloseq_clean, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_seqcenter$contaminant)

# FALSE  TRUE 0.1 threshold
# 7557     719 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(multigen_phyloseq_clean, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_seqcenter$contaminant)

seqcenter_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="Seqcenter Decontam")

seqcenter_decontam_plot
  
ggsave("/Users/Abby/OneDrive - Michigan State University/Epigenetics/sequencing/Rhizo_0.5_decontam_plot.png", plot=rhizodecontam_0.5_plot,width = 7,height = 5,units = "in")

# not saving the seqcenter decontam into the full dataset right now



```

### Extraction group g0_1
```{r g0_1}
# Decontam by extraction group - g0_1
phyloseq_g0_1 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g0_1", .keep_all_taxa = TRUE)

sample_data(phyloseq_g0_1)$is.neg <- sample_data(phyloseq_g0_1)$sample_type == "negative"
contamdf.prev0.1_g0_1 <- isContaminant(phyloseq_g0_1, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g0_1$contaminant)

# FALSE  TRUE 0.1 threshold
# 7973     303 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g0_1, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g0_1$contaminant)

g0_1_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g0_1 Decontam")

g0_1_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_1.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_1.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g0_1 <- prune_taxa(keep.taxa, phyloseq_g0_1)
decontam_phyloseq_g0_1

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g0_1), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g0_1 <- prune_samples(keep.samples, decontam_phyloseq_g0_1)
decontam_phyloseq_g0_1

write.csv(otu_table(decontam_phyloseq_g0_1), "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/decontam_ASVg0_1.csv")


# decontam with positive control

sample_data(decontam_phyloseq_g0_1)$is.neg <- sample_data(decontam_phyloseq_g0_1)$sample_type == "positive"
contamdf.prev0.1_g0_1 <- isContaminant(decontam_phyloseq_g0_1, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g0_1$contaminant)

# FALSE  TRUE 0.1 threshold
# 7962     11 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g0_1, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g0_1$contaminant)

g0_1pos_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g0_1 Decontam")

g0_1pos_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_1pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_1pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g0_1 <- prune_taxa(keep.taxa, decontam_phyloseq_g0_1)
decontam_phyloseq_g0_1

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g0_1), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g0_1 <- prune_samples(keep.samples, decontam_phyloseq_g0_1)
decontam_phyloseq_g0_1

write.csv(otu_table(decontam_phyloseq_g0_1), "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/decontam_ASVg0_1pos.csv")

```

# Extraction group g0_2
```{r g0_2}
# Decontam by extraction group - g0_2
phyloseq_g0_2 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g0_2", .keep_all_taxa = TRUE)

sample_data(phyloseq_g0_2)$is.neg <- sample_data(phyloseq_g0_2)$sample_type == "negative"
contamdf.prev0.1_g0_2 <- isContaminant(phyloseq_g0_2, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g0_2$contaminant)

# FALSE  TRUE 0.1 threshold
# 8200     76 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g0_2, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g0_2$contaminant)

g0_2_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g0_2 Decontam")

g0_2_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_2.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_2.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g0_2 <- prune_taxa(keep.taxa, phyloseq_g0_2)
decontam_phyloseq_g0_2

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g0_2), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g0_2 <- prune_samples(keep.samples, decontam_phyloseq_g0_2)
decontam_phyloseq_g0_2


# decontam with positive control

sample_data(decontam_phyloseq_g0_2)$is.neg <- sample_data(decontam_phyloseq_g0_2)$sample_type == "positive"
contamdf.prev0.1_g0_2 <- isContaminant(decontam_phyloseq_g0_2, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g0_2$contaminant)

# FALSE  TRUE 0.1 threshold
# 8185     15 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g0_2, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g0_2$contaminant)

g0_2pos_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g0_2 Decontam")

g0_2pos_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_2pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g0_2pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g0_2 <- prune_taxa(keep.taxa, decontam_phyloseq_g0_2)
decontam_phyloseq_g0_2

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g0_2), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g0_2 <- prune_samples(keep.samples, decontam_phyloseq_g0_2)
decontam_phyloseq_g0_2

write.csv(otu_table(decontam_phyloseq_g0_1), "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/decontam_ASVg0_1pos.csv")
```
# Extraction group g1_1
```{r g1_1}
# Decontam by extraction group - g1_1
phyloseq_g1_1 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g1_1", .keep_all_taxa = TRUE)

sample_data(phyloseq_g1_1)$is.neg <- sample_data(phyloseq_g1_1)$sample_type == "negative"
contamdf.prev0.1_g1_1 <- isContaminant(phyloseq_g1_1, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g1_1$contaminant)

# FALSE  TRUE 0.1 threshold
# 8253     23 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g1_1, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g1_1$contaminant)

g1_1_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g1_1 Decontam")

g1_1_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_1.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_1.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g1_1 <- prune_taxa(keep.taxa, phyloseq_g1_1)
decontam_phyloseq_g1_1

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g1_1), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g1_1 <- prune_samples(keep.samples, decontam_phyloseq_g1_1)
decontam_phyloseq_g1_1


# decontam with positive control

sample_data(decontam_phyloseq_g1_1)$is.neg <- sample_data(decontam_phyloseq_g1_1)$sample_type == "positive"
contamdf.prev0.1_g1_1 <- isContaminant(decontam_phyloseq_g1_1, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g1_1$contaminant)

# FALSE  TRUE 0.1 threshold
# 8243     10 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g1_1, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g1_1$contaminant)

g1_1pos_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g1_1 pos Decontam")

g1_1pos_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_1pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_1pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g1_1 <- prune_taxa(keep.taxa, decontam_phyloseq_g1_1)
decontam_phyloseq_g1_1

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g1_1), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g1_1 <- prune_samples(keep.samples, decontam_phyloseq_g1_1)
decontam_phyloseq_g1_1


```
# Extraction group g1_2
```{r g1_2}
# Decontam by extraction group - g1_2
phyloseq_g1_2 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g1_2", .keep_all_taxa = TRUE)

sample_data(phyloseq_g1_2)$is.neg <- sample_data(phyloseq_g1_2)$sample_type == "negative"
contamdf.prev0.1_g1_2 <- isContaminant(phyloseq_g1_2, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g1_2$contaminant)

# FALSE  TRUE 0.1 threshold
# 8269     7 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g1_2, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g1_2$contaminant)

g1_2_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g1_2 Decontam")

g1_2_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_2.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_2.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g1_2 <- prune_taxa(keep.taxa, phyloseq_g1_2)
decontam_phyloseq_g1_2

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g1_2), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g1_2 <- prune_samples(keep.samples, decontam_phyloseq_g1_2)
decontam_phyloseq_g1_2


# decontam with positive control

sample_data(decontam_phyloseq_g1_2)$is.neg <- sample_data(decontam_phyloseq_g1_2)$sample_type == "positive"
contamdf.prev0.1_g1_2 <- isContaminant(decontam_phyloseq_g1_2, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g1_2$contaminant)

# FALSE  TRUE 0.1 threshold
# 8257     12 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g1_2, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g1_2$contaminant)


write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_2pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_2pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g1_2 <- prune_taxa(keep.taxa, decontam_phyloseq_g1_2)
decontam_phyloseq_g1_2

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g1_2), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g1_2 <- prune_samples(keep.samples, decontam_phyloseq_g1_2)
decontam_phyloseq_g1_2


```

# Extraction group g1_3
```{r g1_3}
# Decontam by extraction group - g1_3
phyloseq_g1_3 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g1_3", .keep_all_taxa = TRUE)

sample_data(phyloseq_g1_3)$is.neg <- sample_data(phyloseq_g1_3)$sample_type == "negative"
contamdf.prev0.1_g1_3 <- isContaminant(phyloseq_g1_3, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g1_3$contaminant)

# FALSE  TRUE 0.1 threshold
# 8271     5 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g1_3, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g1_3$contaminant)

g1_3_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g1_3 Decontam")

g1_3_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_3.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_3.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g1_3 <- prune_taxa(keep.taxa, phyloseq_g1_3)
decontam_phyloseq_g1_3

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g1_3), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g1_3 <- prune_samples(keep.samples, decontam_phyloseq_g1_3)
decontam_phyloseq_g1_3


# decontam with positive control

sample_data(decontam_phyloseq_g1_3)$is.neg <- sample_data(decontam_phyloseq_g1_3)$sample_type == "positive"
contamdf.prev0.1_g1_3 <- isContaminant(decontam_phyloseq_g1_3, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g1_3$contaminant)

# FALSE  TRUE 0.1 threshold
# 8263     8 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g1_3, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g1_3$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_3pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g1_3pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g1_3 <- prune_taxa(keep.taxa, decontam_phyloseq_g1_3)
decontam_phyloseq_g1_3

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g1_3), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g1_3 <- prune_samples(keep.samples, decontam_phyloseq_g1_3)
decontam_phyloseq_g1_3


```

# Extraction group g2_0714
```{r g2_0714}
# Decontam by extraction group - g2_0714
phyloseq_g2_0714 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0714", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0714)$is.neg <- sample_data(phyloseq_g2_0714)$sample_type == "negative"
contamdf.prev0.1_g2_0714 <- isContaminant(phyloseq_g2_0714, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0714$contaminant)

# FALSE  TRUE 0.1 threshold
# 8275     1 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0714, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0714$contaminant)

g2_0714_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0714 Decontam")

g2_0714_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0714.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0714.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0714 <- prune_taxa(keep.taxa, phyloseq_g2_0714)
decontam_phyloseq_g2_0714

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0714), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0714 <- prune_samples(keep.samples, decontam_phyloseq_g2_0714)
decontam_phyloseq_g2_0714


# decontam with positive control

sample_data(decontam_phyloseq_g2_0714)$is.neg <- sample_data(decontam_phyloseq_g2_0714)$sample_type == "positive"
contamdf.prev0.1_g2_0714 <- isContaminant(decontam_phyloseq_g2_0714, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0714$contaminant)

# FALSE  TRUE 0.1 threshold
# 8258     17 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g2_0714, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0714$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0714pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0714pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0714 <- prune_taxa(keep.taxa, decontam_phyloseq_g2_0714)
decontam_phyloseq_g2_0714

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0714), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0714 <- prune_samples(keep.samples, decontam_phyloseq_g2_0714)
decontam_phyloseq_g2_0714
```

# Extraction group g2_0114
```{r g2_0114}
# Decontam by extraction group - g2_0114
phyloseq_g2_0114 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0114", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0114)$is.neg <- sample_data(phyloseq_g2_0114)$sample_type == "negative"
contamdf.prev0.1_g2_0114 <- isContaminant(phyloseq_g2_0114, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0114$contaminant)

# FALSE  TRUE 0.1 threshold
# 8226     50 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0114, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0114$contaminant)

g2_0114_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0114 Decontam")

g2_0114_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0114.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0114.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0114 <- prune_taxa(keep.taxa, phyloseq_g2_0114)
decontam_phyloseq_g2_0114

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0114), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0114 <- prune_samples(keep.samples, decontam_phyloseq_g2_0114)
decontam_phyloseq_g2_0114


# decontam with positive control

sample_data(decontam_phyloseq_g2_0114)$is.neg <- sample_data(decontam_phyloseq_g2_0114)$sample_type == "positive"
contamdf.prev0.1_g2_0114 <- isContaminant(decontam_phyloseq_g2_0114, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0114$contaminant)

# FALSE  TRUE 0.1 threshold
# 8214     12 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g2_0114, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0114$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0114pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0114pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0114 <- prune_taxa(keep.taxa, decontam_phyloseq_g2_0114)
decontam_phyloseq_g2_0114

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0114), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0114 <- prune_samples(keep.samples, decontam_phyloseq_g2_0114)
decontam_phyloseq_g2_0114
```

# Extraction group g2_0121
```{r g2_0121}
# Decontam by extraction group - g2_0121
phyloseq_g2_0121 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0121", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0121)$is.neg <- sample_data(phyloseq_g2_0121)$sample_type == "negative"
contamdf.prev0.1_g2_0121 <- isContaminant(phyloseq_g2_0121, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0121$contaminant)

# FALSE  TRUE 0.1 threshold
# 8246     30 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0121, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0121$contaminant)

g2_0121_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0121 Decontam")

g2_0121_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0121.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0121.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0121 <- prune_taxa(keep.taxa, phyloseq_g2_0121)
decontam_phyloseq_g2_0121

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0121), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0121 <- prune_samples(keep.samples, decontam_phyloseq_g2_0121)
decontam_phyloseq_g2_0121


# decontam with positive control

sample_data(decontam_phyloseq_g2_0121)$is.neg <- sample_data(decontam_phyloseq_g2_0121)$sample_type == "positive"
contamdf.prev0.1_g2_0121 <- isContaminant(decontam_phyloseq_g2_0121, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0121$contaminant)

# FALSE  TRUE 0.1 threshold
# 8239     7 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g2_0121, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0121$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0121pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0121pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0121 <- prune_taxa(keep.taxa, decontam_phyloseq_g2_0121)
decontam_phyloseq_g2_0121

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0121), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0121 <- prune_samples(keep.samples, decontam_phyloseq_g2_0121)
decontam_phyloseq_g2_0121
```

# Extraction group g2_0125
```{r g2_0125}
# Decontam by extraction group - g2_0125
phyloseq_g2_0125 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0125", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0125)$is.neg <- sample_data(phyloseq_g2_0125)$sample_type == "negative"
contamdf.prev0.1_g2_0125 <- isContaminant(phyloseq_g2_0125, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0125$contaminant)

# FALSE  TRUE 0.1 threshold
# 8273     3 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0125, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0125$contaminant)

g2_0125_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0125 Decontam")

g2_0125_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0125.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0125.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0125 <- prune_taxa(keep.taxa, phyloseq_g2_0125)
decontam_phyloseq_g2_0125

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0125), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0125 <- prune_samples(keep.samples, decontam_phyloseq_g2_0125)
decontam_phyloseq_g2_0125


# decontam with positive control

sample_data(decontam_phyloseq_g2_0125)$is.neg <- sample_data(decontam_phyloseq_g2_0125)$sample_type == "positive"
contamdf.prev0.1_g2_0125 <- isContaminant(decontam_phyloseq_g2_0125, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0125$contaminant)

# FALSE  TRUE 0.1 threshold
# 8251     22 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g2_0125, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0125$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0125pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0125pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0125 <- prune_taxa(keep.taxa, decontam_phyloseq_g2_0125)
decontam_phyloseq_g2_0125

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0125), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0125 <- prune_samples(keep.samples, decontam_phyloseq_g2_0125)
decontam_phyloseq_g2_0125
```

# Extraction group g2_0709
```{r g2_0709}
# Decontam by extraction group - g2_0709
phyloseq_g2_0709 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0709", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0709)$is.neg <- sample_data(phyloseq_g2_0709)$sample_type == "negative"
contamdf.prev0.1_g2_0709 <- isContaminant(phyloseq_g2_0709, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0709$contaminant)

# FALSE  TRUE 0.1 threshold
# 8203     73 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0709, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0709$contaminant)

g2_0709_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0709 Decontam")

g2_0709_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0709.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0709.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0709 <- prune_taxa(keep.taxa, phyloseq_g2_0709)
decontam_phyloseq_g2_0709

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0709), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0709 <- prune_samples(keep.samples, decontam_phyloseq_g2_0709)
decontam_phyloseq_g2_0709


# decontam with positive control

sample_data(decontam_phyloseq_g2_0709)$is.neg <- sample_data(decontam_phyloseq_g2_0709)$sample_type == "positive"
contamdf.prev0.1_g2_0709 <- isContaminant(decontam_phyloseq_g2_0709, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0709$contaminant)

# FALSE  TRUE 0.1 threshold
# 8190     13 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g2_0709, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0709$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0709pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0709pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0709 <- prune_taxa(keep.taxa, decontam_phyloseq_g2_0709)
decontam_phyloseq_g2_0709

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0709), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0709 <- prune_samples(keep.samples, decontam_phyloseq_g2_0709)
decontam_phyloseq_g2_0709
```

# Extraction group g2_0720
```{r g2_0720}
# Decontam by extraction group - g2_0720
phyloseq_g2_0720 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0720", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0720)$is.neg <- sample_data(phyloseq_g2_0720)$sample_type == "negative"
contamdf.prev0.1_g2_0720 <- isContaminant(phyloseq_g2_0720, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0720$contaminant)

# FALSE  TRUE 0.1 threshold
# 8208     68 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0720, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0720$contaminant)

g2_0720_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0720 Decontam")

g2_0720_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0720.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0720.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0720 <- prune_taxa(keep.taxa, phyloseq_g2_0720)
decontam_phyloseq_g2_0720

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0720), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0720 <- prune_samples(keep.samples, decontam_phyloseq_g2_0720)
decontam_phyloseq_g2_0720


# decontam with positive control

sample_data(decontam_phyloseq_g2_0720)$is.neg <- sample_data(decontam_phyloseq_g2_0720)$sample_type == "positive"
contamdf.prev0.1_g2_0720 <- isContaminant(decontam_phyloseq_g2_0720, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0720$contaminant)

# FALSE  TRUE 0.1 threshold
# 8193     15 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(decontam_phyloseq_g2_0720, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0720$contaminant)

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0720pos.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0720pos.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0720 <- prune_taxa(keep.taxa, decontam_phyloseq_g2_0720)
decontam_phyloseq_g2_0720

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0720), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0720 <- prune_samples(keep.samples, decontam_phyloseq_g2_0720)
decontam_phyloseq_g2_0720
```

# Extraction group g2_0804
```{r g2_0804}
# Decontam by extraction group - g2_0804
phyloseq_g2_0804 <- multigen_phyloseq_clean %>% 
  ps_filter(extr_group == "g2_0804", .keep_all_taxa = TRUE)

sample_data(phyloseq_g2_0804)$is.neg <- sample_data(phyloseq_g2_0804)$sample_type == "negative"
contamdf.prev0.1_g2_0804 <- isContaminant(phyloseq_g2_0804, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0804$contaminant)

# FALSE  TRUE 0.1 threshold
# 8239     37 

# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(phyloseq_g2_0804, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)

# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
                    contaminant=contamdf.prev0.1_g2_0804$contaminant)

g2_0804_decontam_plot <- ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + labs(x="Prevalence (Negative Controls)", y="Prevalence (True Samples)", title="g2_0804 Decontam")

g2_0804_decontam_plot

write.csv(df.pa, "/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0804.csv")

# open contaminant-table in excel and add column name OTUID, then read back in below

##removing contaminants from phyloseq object
df.pa <- read.csv("/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/contaminant-table-g2_0804.csv")
subset.df <- subset(df.pa, contaminant== "FALSE")
keep.taxa <- as.vector(subset.df$OTUID)

subset.df.remove <- subset(df.pa, contaminant== "TRUE")
remove.taxa.list <- as.vector(subset.df.remove)
remove.taxa.df <- data.frame(remove.taxa.list)

decontam_phyloseq_g2_0804 <- prune_taxa(keep.taxa, phyloseq_g2_0804)
decontam_phyloseq_g2_0804

##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0804), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0804 <- prune_samples(keep.samples, decontam_phyloseq_g2_0804)
decontam_phyloseq_g2_0804


# decontam with positive control

sample_data(decontam_phyloseq_g2_0804)$is.neg <- sample_data(decontam_phyloseq_g2_0804)$sample_type == "positive"
contamdf.prev0.1_g2_0804 <- isContaminant(decontam_phyloseq_g2_0804, method="prevalence", neg="is.neg")
table(contamdf.prev0.1_g2_0804$contaminant)

# FALSE  TRUE 0.1 threshold
# 8239     0 


##filter metadata to true samples only then merge with phyloseq decontam (with removed OTU contaminants)
subset.metadata<- subset(sample_data(decontam_phyloseq_g2_0804), is.neg ==FALSE)

keep.samples <- as.vector(subset.metadata$SampleID)
keep.samples

decontam_phyloseq_g2_0804 <- prune_samples(keep.samples, decontam_phyloseq_g2_0804)
decontam_phyloseq_g2_0804
```

# Merge phyloseq objects back together
```{r re-merge}
seed_phyloseq_decontam <- merge_phyloseq(decontam_phyloseq_g0_1, decontam_phyloseq_g0_2, decontam_phyloseq_g1_1, decontam_phyloseq_g1_2, decontam_phyloseq_g1_3, decontam_phyloseq_g2_0114, decontam_phyloseq_g2_0121, decontam_phyloseq_g2_0125, decontam_phyloseq_g2_0709, decontam_phyloseq_g2_0714, decontam_phyloseq_g2_0720, decontam_phyloseq_g2_0804)

seed_phyloseq_decontam

saveRDS(seed_phyloseq_decontam, file="/Users/Abby/OneDrive - Michigan State University/multigen_analysis/seeds/New_argonne_seqs/seed_phyloseq_decontam.rds", compress = TRUE)
```