knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(ggplot2)
library(readr)
library(ggpubr)
library(vegan)
library(scales)
library(lme4)
library(lmerTest)
library(pbkrtest)
library(emmeans)
library(rstatix)
library(forcats)
library(scales)
metadata <- read_tsv("../03_metadata/metadata.txt")
# Alpha diversity metrics
shannon <- read_tsv("/Users/victoria/Documents/ANEQ505_test/r_decomp_tutorial/04_code/alpha_div/shannon.tsv")
evenness <- read_tsv("alpha_div/evenness.tsv")
observed_features <- read_tsv("alpha_div/observed_features.tsv")
faiths_pd <- read_tsv("alpha_div/faith_pd.tsv")
# Shannon
shannon <- shannon %>%
rename(sample_name = 1)  #rename 1st column to to match metadata
# Evenness
evenness <- evenness %>%
rename(sample_name = 1)  #rename 1st column to to match metadata
# Observed Features
observed_features <- observed_features %>%
rename(sample_name = 1)  #rename 1st column to to match metadata
# Faiths PD
faiths_pd <- faiths_pd %>%
rename(sample_name = 1)  #rename 1st column to to match metadata
# Add to metadata
metadata <- metadata %>%
left_join(shannon, by = "sample_name") %>%
left_join(evenness, by = "sample_name") %>%
left_join(observed_features, by = "sample_name") %>%
left_join(faiths_pd, by = "sample_name")
metadata <- metadata %>%
drop_na(shannon_entropy, pielou_evenness, observed_features, faith_pd)
bray_curtis <- read_tsv("beta_div/bray_curtis.txt")
jaccard <- read_tsv("beta_div/jaccard.txt")
uw_unifrac <- read_tsv("beta_div/unweighted_unifrac.txt")
w_unifrac <- read_tsv("beta_div/weighted_unifrac.txt")
View(bray_curtis)
View(observed_features)
View(uw_unifrac)
View(w_unifrac)
View(uw_unifrac)
bray_curtis_var <- bray_curtis %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
View(bray_curtis_var)
bray_curtis <- bray_curtis %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
View(bray_curtis)
# Jaccard
# get variance
jaccard_var <- jaccard %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Clean ordinace
jaccard <- jaccard %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_jaccard <- jaccard %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
# Unweighted Unifrac
# Get variance
uw_var <- uw_unifrac %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Clean ordinance
uw_unifrac <- uw_unifrac %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_uwunifrac <- uw_unifrac %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
# Weighted Unifrac
# get variance
wunifrac_var <- w_unifrac %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Tidy ordinace
w_unifrac <- w_unifrac %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_wunifrac <- w_unifrac %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
knitr::opts_chunk$set(echo = TRUE)
jaccard <- read_tsv("beta_div/jaccard.txt")
library(tidyverse)
library(ggplot2)
library(readr)
library(ggpubr)
library(vegan)
library(scales)
library(lme4)
library(lmerTest)
library(pbkrtest)
library(emmeans)
library(rstatix)
library(forcats)
library(scales)
jaccard <- read_tsv("beta_div/jaccard.txt")
jaccard_var <- jaccard %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(column_last(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
jaccard_var <- jaccard %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Read in beta div metrics
bray_curtis <- read_tsv("beta_div/bray_curtis.txt")
jaccard <- read_tsv("beta_div/jaccard.txt")
uw_unifrac <- read_tsv("beta_div/unweighted_unifrac.txt")
w_unifrac <- read_tsv("beta_div/weighted_unifrac.txt")
# Clean metrics
# bray curtis
# get variance
bray_curtis_var <- bray_curtis %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(lost_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Read in beta div metrics
bray_curtis <- read_tsv("beta_div/bray_curtis.txt")
jaccard <- read_tsv("beta_div/jaccard.txt")
uw_unifrac <- read_tsv("beta_div/unweighted_unifrac.txt")
w_unifrac <- read_tsv("beta_div/weighted_unifrac.txt")
# Clean metrics
# bray curtis
# get variance
bray_curtis_var <- bray_curtis %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# clean ordination
bray_curtis <- bray_curtis %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_bc <- bray_curtis %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
metadata <- read_tsv("../03_metadata/metadata.txt")
# Read in beta div metrics
bray_curtis <- read_tsv("beta_div/bray_curtis.txt")
jaccard <- read_tsv("beta_div/jaccard.txt")
uw_unifrac <- read_tsv("beta_div/unweighted_unifrac.txt")
w_unifrac <- read_tsv("beta_div/weighted_unifrac.txt")
# Clean metrics
# bray curtis
# get variance
bray_curtis_var <- bray_curtis %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# clean ordination
bray_curtis <- bray_curtis %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_bc <- bray_curtis %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
# Jaccard
# get variance
jaccard_var <- jaccard %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Clean ordinace
jaccard <- jaccard %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_jaccard <- jaccard %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
# Unweighted Unifrac
# Get variance
uw_var <- uw_unifrac %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Clean ordinance
uw_unifrac <- uw_unifrac %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_uwunifrac <- uw_unifrac %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
# Weighted Unifrac
# get variance
wunifrac_var <- w_unifrac %>%
slice(3) %>%   # row 3 contains the numeric proportion explained values
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%  # split tab-separated values
mutate(across(starts_with("PC"), ~ as.numeric(.x) * 100)) %>%  # convert to percent
select(PC1, PC2, PC3)
# Tidy ordinace
w_unifrac <- w_unifrac %>%
slice(6:n()) %>%
filter(!Eigvals %in% c("Biplot", "Site constraints")) %>%
separate(last_col(), into = paste0("PC", 1:10), sep = "\t") %>%
rename(sample_name = Eigvals) #change eigenvals to match the sampleid, sample_name etc in your metadata
metadata_wunifrac <- w_unifrac %>%
select(sample_name, PC1, PC2, PC3) %>%
left_join(metadata) %>%
mutate(across(c(PC1, PC2, PC3), as.numeric))
View(metadata_jaccard)
knitr::opts_chunk$set(echo = TRUE)
packages_to_install <- c("tidyverse", "ggplot2","readr","ggpubr","vegan","scales","lme4", "lmerTest","pbkrtest","emmeans", "rstatix", "forcats", "scales")
# Each of the packages we will use are CRAN packages
packages_to_install <- c("tidyverse", "ggplot2","readr","ggpubr","vegan","scales","lme4", "lmerTest","pbkrtest","emmeans", "rstatix", "forcats", "scales")
