PalaueDNA_Nov2021_2023July4.Rmd

---
title: "PalaueDNA_Nov2021"
author: "Collin J. Closek"
date: "2023-02-27"
output: html_document
---

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


## Load required packages 
```{r }
#install.packages("tidyverse")
#install.packages("vegan")
#install.packages("devtools")
#devtools::install_github("rlbarter/superheat")
#install.packages("here")
install.packages("phyloseq")
#install.packages("owmr")
library(tidyverse)
library(tidyr)
library(ggplot2)
library(dplyr)
library(vegan)
library(data.table)
library(superheat)
library(here)

#if(!requireNamespace("BiocManager")){
#  install.packages("BiocManager")
#}
#BiocManager::install("phyloseq")
library(phyloseq)
#packageVersion("phyloseq")
library(owmr)

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")

#BiocManager::install("multtest")
#BiocManager::install("phyloseq")

#if (!requireNamespace("devtools", quietly = TRUE))
#      install.packages('devtools')
library(devtools)

#devtools::install_github("gauravsk/ranacapa")
library(ranacapa)
#ranacapa::runRanacapaApp()
BiocManager::install("MicrobiotaProcess")
library(MicrobiotaProcess)
library(patchwork)
```

##Load data
```{r}
#esv_data_all <- (read.delim("./data/JVB1710-MiFishU-read-data.csv", sep=",", check.names=FALSE, header = TRUE))

esv_data_all <- read.csv(here("data/JVB1710-MiFishU-read-data.csv"),header = TRUE)

head (esv_data_all)
#glimpse (esv_data_all)
```

##Remove preliminary run (i.e., remove duplicated samples "S046...1" that have newer ".2" version) & cleanup columns
```{r}
data <- select(esv_data_all,-c(S046420.1, S046419.1, S046417.1, S046416.1, S046415.1, S046414.1, X))
data <- data %>%
  rename(
    Match_percent=X..match,
    Species_number=X..species
  ) 

#glimpse(data)
```


##Cobine header row cells into one 
Takes the following row headers: "TestId", "ESVId", "sequence", "Kingdom",  "Phylum", "Class", "Order", "Family", "Genus", "Species", "% match", "# species"
Merges them into one "MergedID" with the separator ","
```{r}
#esv_data_all_mergedrows <- unite(esv_data_all, col='MergedID', c('TestId', 'ESVId', 'sequence', 'Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species', '% match', '# species'), sep=',')
```


##Create long data version
```{r}
data_long <-
  data %>%
  pivot_longer(
    "S049508.1":"S049391.0",
    names_to = "sample",
    values_to = "obs"
  )

#glimpse(data_long)
```



##Bring in samples key & merge with long data
```{r}
samples <- read.csv(here("data/JVB1710-samples_rv3.csv"),header=TRUE)
glimpse(samples)

left_joined <- merge(data_long, samples, by = "sample")
#glimpse(left_joined)
```

##Create full sample ID - merge original Sample name w/ the JV sample version name
```{r}
left_joined_unite <- left_joined %>%
  unite(col="ID_v", SampleID, sample, sep= "-", 
        remove = TRUE)

#glimpse(left_joined_unite)
```

##Bring data frame back into Wide Format
```{r}
esv_read_data <- 
  left_joined_unite %>% 
  pivot_wider(
    names_from = "ID_v", 
    values_from = "obs"
  )

write.csv(esv_read_data, file = "/Users/closek/Documents/Nov2021_PalaueDNA/R/products/esv_read_data_rv.csv")
```

##Identify and verify that higher taxonomic rankings check out 
```{r}
unique(esv_read_data$Phylum)
#Should only bring back Chordata
#brings back: "Chordata"    ""
#    "" is n=179 unannotated ESVs
unique(esv_read_data$Class)
unique(esv_read_data$Class)
#"Actinopteri"    ""               "Mammalia" (n=3spp - false killer whale, spinner dolphin, and pantropical spotted dolphin)      "Chondrichthyes" (n=1spp - Blue shark), "Aves" (n=1spp ~ game fowl/pheasant, chicken)           "Amphibia"  (n=1spp - giant slamander)
#Select for only Classes: Actinopteri and Mammalia, which will remove Classes: Aves and Amphibia, as well as the unannotated ESVs 
#Within Mammalia remove non-marine mammals in  Family: Bovidae, Suidae, Canidae, Felindae, Hominidae and Order: Primates
```


##Create Nearshore & Offshore dataframes
```{r}
offshore_esv <- esv_read_data %>% select(2, 8, 9, 10, contains('COS'), contains('AFZ'), contains('DFZ'), contains('NMS'))

nearshore_esv <- esv_read_data %>% select(2, 8, 9, 10, contains('COS'), contains('SHR'), contains('LGN'), contains('FRR'))

write.csv(offshore_esv, file = "/Users/closek/Documents/Nov2021_PalaueDNA/R/products/offshore_esv.csv")
write.csv(nearshore_esv, file = "/Users/closek/Documents/Nov2021_PalaueDNA/R/products/nearshore_esv.csv")
```

##Create Nearshore & Offshore species lists
```{r}
offshore_esv_sp <- esv_read_data %>% select(2, 10, contains('COS'), contains('AFZ'), contains('DFZ'), contains('NMS'))

nearshore_esv_sp <- esv_read_data %>% select(2, 10, contains('COS'), contains('SHR'), contains('LGN'), contains('FRR'))

write.csv(offshore_esv_sp, file = "/Users/closek/Documents/Nov2021_PalaueDNA/R/products/offshore_esv_sp.csv")
write.csv(nearshore_esv_sp, file = "/Users/closek/Documents/Nov2021_PalaueDNA/R/products/nearshore_esv_sp.csv")
```


##Unique nearshore species
```{r}
ns_esv <- nearshore_esv[rowSums(nearshore_esv==0, na.rm=TRUE)<ncol(nearshore_esv), ]
#nearshore_esv[rowSums(nearshore_esv[])>0,]
n_distinct(nearshore_esv$Species)
unique(nearshore_esv$Species)


nearshore_esv_sp_wo0 <- read.csv(here("products/nearshore_esv_sp_wo0.csv"),header = TRUE)
unique(nearshore_esv_sp_wo0$Species)
n_distinct(nearshore_esv_sp_wo0$Species)
```


##Unique offshore species
```{r}
n_distinct(offshore_esv$Species)
unique(offshore_esv$Species)

offshore_esv_sp_wo0 <- read.csv(here("products/offshore_esv_sp_wo0.csv"),header = TRUE)
unique(offshore_esv_sp_wo0$Species)
n_distinct(offshore_esv_sp_wo0$Species)
```


##Unique species Nearshore & Offshore
```{r}
n_distinct(esv_read_data$Species)
```


##Unique species
```{r}
#esv_list <- read.csv(here("data/JVB1710-MiFishU-esv-data.csv"),header = TRUE)
#glimpse(esv_list)
#str(esv_list)

#n_distinct(esv_list$Species)
#unique(esv_list$Species)

n_distinct(esv_read_data$Species)
unique(esv_read_data$Species)
```

##Unique Genera
```{r}
#n_distinct(esv_list$Genus)
#unique(esv_list$Genus)

n_distinct(esv_read_data$Genus)
unique(esv_read_data$Genus)
```

##Unique Families
```{r}
#n_distinct(esv_list$Family)
#unique(esv_list$Family)
#unique

n_distinct(esv_read_data$Family)
unique(esv_read_data$Family)

```


##Transform esv_read_data into Phyloseq Object
```{r}
names(esv_read_data)

#create an OTU/ESV count table
esv_read_count <- select(esv_read_data, -c(1))
esv_read_count <- select(esv_read_count, -c(2:18))
esv_read_count <- data.frame(esv_read_count, row.names = 1)
#names(esv_read_count)
#remove_prefix(esv_read_count, c("X")) %>% names()
class(esv_read_count)
esv_read_count_m <- as.matrix(esv_read_count)
esv_read_count_m
class(esv_read_count_m)

rowSums(esv_read_count)
colSums(esv_read_count)

ESV <- otu_table(esv_read_count_m, taxa_are_rows = TRUE)

write.csv(ESV, file = "/Users/closek/Documents/Nov2021_PalaueDNA/R/products/ESV.csv")

#create a taxonomy table
esv_read_tax <- select(esv_read_data, -c(1,3))
esv_read_tax <- select(esv_read_tax, -c(9:194))
esv_read_tax <- data.frame(esv_read_tax, row.names = 1)
class(esv_read_tax)
esv_read_tax_m <- as.matrix(esv_read_tax)
class(esv_read_tax_m)
names(esv_read_tax_m)
TAX <- tax_table(esv_read_tax_m)



#create a sample metadata 
samp_meta <- read.csv(here("data/SampleMetaData.csv"),header = TRUE)
samp_meta_df <- data.frame(samp_meta, row.names = 1)
samp_metadata <- sample_data(samp_meta_df)
class(samp_metadata)

#create phyloseq object
physeq <- phyloseq(ESV, TAX)
physeq

physeq1 <- merge_phyloseq(physeq, samp_metadata)
#physeq1 <- phyloseq(ESV, TAX, samp_metadata)
physeq1
```

##Components of phyloseq object
```{r}
sample_names(physeq1)
rank_names(physeq1)
sample_variables(physeq1)
```

#Subset samples to keep only Environmental samples
```{r}
physeq1 <- subset_samples(physeq1, Control_or_Envtl =="Envtl")
physeq1
```


```{r}

plot_bar(physeq1, fill = "Family")

```

##Subset to keep only ESVs with marine fish and mammal annotations
```{r}
#Select for only Classes: Actinopteri and Mammalia, which will remove Classes: Aves and Amphibia, as well as the unannotated ESVs 
#Within Mammalia remove non-marine mammals in  Family: Bovidae, Suidae, Canidae, Felindae, Hominidae and Order: Primates
physeq1
physeq1_mar <- subset_taxa(physeq1, Class %in% c("Actinopteri", "Mammalia"))
physeq1_mar
physeq1_mar <- subset_taxa(physeq1, !(Order %in% c("Primates")))
physeq1_mar
physeq1_mar <- subset_taxa(physeq1, !(Family %in% c("Bovidae", "Canidae", "Felindae", "Hominidae", "Suidae")))
physeq1_mar
```



```{r}
sample_sums(physeq1_mar)
```

```{r}
sample_sums(physeq1_mar)<=10000
```


```{r}
physeq1_mar_pruned <- prune_samples(sample_sums(physeq1_mar)>=10000, physeq1_mar)
physeq1_mar_pruned
```


##Normalize number of reads
```{r}
total = median(sample_sums(physeq1_mar_pruned))
standf = function(x, t=total) round(t *(x / sum(x)))
physeq2 = transform_sample_counts(physeq1_mar_pruned, standf)
physeq2
```


```{r}

#rcurve(physeq2, add_sample_data = TRUE)
#p <- ggrare(physeq2, step = 500, color = "Zone", label = "Sample", se = FALSE)

sample_data(physeq2)
#TopNESVs <- names(sort(taxa_sums(physeq2), TRUE)[1:30])
#physeq2_30 <- prune_taxa(TopNESVs, physeq2$Species)
#print(physeq2_30)
```



```{r}
plot_bar(physeq2, fill = "Family")
```

```{r}
plot_bar(physeq2, fill = "Order")
```


```{r}
plot_heatmap(physeq1, taxa.label="Phylum")
```


```{r}
plot_heatmap(physeq2, taxa.label="Phylum")
```


```{r}
plot_bar(physeq2, x= "Zone", fill = "Family")
```

```{r}
plot_bar(physeq2, x= "Zone", fill = "Family") + facet_wrap(~Cardinal_direction)
```


```{r}
#plot_net(physeq2, maxdist=0.4, color="Cardinal_direction", shape="Near_or_Offshore")
```

```{r}
plot_richness(physeq2, measures = c("Chao1", "Shannon"), x="Near_or_Offshore", color="Near_or_Offshore")

```


```{r}
plot_richness(physeq2, measures=c("Chao1", "Shannon"), x="Near_or_Offshore", color="Zone")
```

```{r}
physeq2.ord <- ordinate(physeq2, "PCoA", "bray")
```


```{r}
po <- plot_ordination(physeq2, physeq2.ord, type = "samples", color = "Zone", shape = "Near_or_Offshore", title = "Samples")
po

jpeg(file="saving_plot1.jpeg")
po <- plot_ordination(physeq2, physeq2.ord, type = "samples", color = "Zone", shape = "Near_or_Offshore", title = "Samples")
dev.off()
```


```{r}
po +facet_wrap(~Near_or_Offshore)
```


```{r}
po +facet_wrap(~Cardinal_direction)
```


```{r}
po +facet_wrap(~Zone)
```


```{r}
sf = filter_taxa(physeq2, function(x) sd(x)/mean(x) > 3.0, TRUE)
sfAct = subset_taxa(sf, Class=="Actinopteri")
sfMam = subset_taxa(sf, Class=="Mammalia")
sfCon = sunbset_taxa(sf, Class=="Chondrichthyes")
```

##Fish
```{r}
plot_bar(sfAct, fill = "Family")


```

##Mammals
```{r}
plot_bar(sfMam, fill = "Species")

jpeg(file="saving_plot2.jpeg")
plot_bar(sfAct, fill = "Family")
dev.off()
```


```{r}
OTUnames10 = names(sort(taxa_sums(sf), TRUE)[1:10])
sf10  = prune_taxa(OTUnames10,  sf)

sf10near = names(subset(sample_data(sf), Near_or_Offshore=="Nearshore"))
sf10off = names(subset(sample_data(sf), Near_or_Offshore=="Offshore"))

top10 <- names(sort(taxa_sums(sf), decreasing = TRUE)[1:10])

tax_table(sf)[top10,]
print(sf10near)
print(sf10off)
```


```{r}
print(sf10near)
```


Rarefaction
Across all locations, what are you identifying?
Samples from one location
Barplot of locations with counts of taxa (species, or rather genera or families)
Top spp in locations
Add tuna spp names to ESVs in TAX?
Add common names/Palauan names


##from https://joey711.github.io/phyloseq/plot_ordination-examples.html
```{r}
#GP = physeq1
```

```{r}
#wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A=0.5*nsamples(GP))
#GP1 = prune_taxa(wh0, GP)
```

```{r}
#GP1 = transform_sample_counts(GP1, function(x) 1E6 * x/sum(x))
```

```{r}
#phylum.sum = tapply(taxa_sums(GP1), tax_table(GP1)[, "Class"], sum, na.rm=TRUE)
#top5phyla = names(sort(phylum.sum, TRUE))[1:5]
#GP1 = prune_taxa((tax_table(GP1)[, "Class"] %in% top5phyla), GP1)
```

```{r}
#GP.ord <- ordinate(GP1, "NMDS", "bray")
#p1 = plot_ordination(GP1, GP.ord, type="taxa", color="Class", title="taxa")
#print(p1)
```

###
```{r}
#metaMDS(GP1)
```

Option to link phyloseq to ranacapa
```{r}

```