#### Load necessary packages ####
library(tidyverse)
library(phyloseq)
library(ape)
library(vegan)
#### Upload files exported from QIIME2 pipeline ####
metafp <- "parkinsons_metadata.txt"
meta <- read_delim(metafp, delim="\t")
otufp <- "feature-table_new.txt"
otu <- read_delim(file = otufp, delim="\t", skip=1)
taxfp <- "taxonomy.tsv"
tax <- read_delim(taxfp, delim="\t")
phylotree <- "rooted_tree.nwk"
phylotree <- read.tree(phylotree)
#### Format Feature Table and Metadata, and Taxonomy data ahead of generating a Phyloseq object ####
## Modify feature table:
otu_mat <- as.matrix(otu[,-1])
rownames(otu_mat) <- otu$`#OTU ID`
OTU <- otu_table(otu_mat, taxa_are_rows = T)
samp_df <- as.data.frame(meta[,-1])
rownames(samp_df)<- meta$`#SampleID`
SAMP <- sample_data(samp_df)
# Reformat taxa table; already split by domain, phylum, class etc,
# but includes confd. values
taxa <- tax %>%
separate(col=`Feature ID Taxon Confidence`, sep=" "
, into = c("Feature ID","Domain"))
taxa <- taxa %>%
separate(col=`...5`, sep=" "
, into = c("...5","ConsensusFamily")) %>%
separate(col=`...6`, sep=" "
, into = c("...6","ConsensusGenus")) %>%
separate(col=`...7`, sep=" "
, into = c("...7","ConsensusSpecies"))
taxa <- taxa %>%
select(-ConsensusFamily, -ConsensusGenus, -ConsensusSpecies)
# Rename column names from the taxa table
colnames(taxa)[3] = "Phylum"
colnames(taxa)[4] = "Class"
colnames(taxa)[5] = "Order"
colnames(taxa)[6] = "Family"
colnames(taxa)[7] = "Genus"
colnames(taxa)[8] ="Species"
taxa <- filter(taxa,
Domain == "d__Bacteria" &
Class!="c__Chloroplast" & Family !="f__Mitochondria")
# Reformat taxa table
taxa_matrix <- as.matrix(taxa[,-1])
rownames(taxa_matrix) <- taxa$`Feature ID`
TAX <- tax_table(taxa_matrix)
#### Make Phyloseq object ####
PARKINSONS <- phyloseq(OTU, SAMP, TAX, phylotree)
# Filter to include only those with PD; we are interested in possible confounders
# among the PD group only
DISEASE_ONLY <- subset_samples(PARKINSONS, Disease == "PD")
# Rarefy samples (Depth based on parameters used by Cirstea et al. (2020)
# for the same dataset)
DISEASE_RARE <- rarefy_even_depth(DISEASE_ONLY, rngseed = 1, sample.size = 3797)
sample_data <- as.data.frame(sample_data(DISEASE_RARE))
