Moving PICRUSt2 data into a phyloseq object

library(phyloseq)
library(tidyverse)

When you finish running the picrust2 pipeline you’ll get an output file called path_abun_unstrat_descrip.tsv. take this file and unzip it to get path_abun_unstrat_descrip.csv. This is the file you’ll want to use for the stats.

It looks like this!

EC_table <- read.csv("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\picrust\\pathways_out\\path_abun_unstrat_descrip.tsv\\path_abun_unstrat_descrip.csv")
EC_table

We’ll say the taxa table is the pathway and description columns, and everything else is your OTU table

EC_taxa <- select(EC_table, c("pathway", "description")) #move columns to taxa table


EC_table <- select(EC_table, -c("pathway", "description")) #remove the sample columns from your OTU table

HEY remember when we removed all the spaces from the names to run picrust? well lets put those back! Picrust also seemed to eat all the dashes and replace them with periods so we’ll change that back too

spacemore <- function(x) {colnames(x) <- gsub("_", " ", colnames(x));x} #making some quick gsup functions because that the way random internet answer I found did it
dashmore <- function(x) {colnames(x) <- gsub("[.]", "-", colnames(x));x}


EC_table <- spacemore(EC_table)
EC_table <- dashmore(EC_table)

#then this is normal phyloseq stuff you know! 
EC_taxa <- as.matrix(EC_taxa)
EC_table <- as.matrix(EC_table)

#check
class(EC_taxa)

[1] "matrix" "array" 

class(EC_table) 

[1] "matrix" "array" 

rownames(EC_table) <- paste0("OTU", 1:nrow(EC_table))
rownames(EC_taxa) <- paste0("OTU", 1:nrow(EC_taxa))

OTU_EC = otu_table(EC_table, taxa_are_rows = TRUE)
TAX_EC = tax_table(EC_taxa)
physeq_EC = phyloseq(OTU_EC, TAX_EC)

Then we can add the sample data and save the file

Key <- read_csv("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\ULTRA\\They call me... data\\Angus-16S\\ultra to categories key.csv")
Key <- data.frame(Key[,-1], row.names=Key$Name)
sampledata = sample_data(Key)

physeq_EC = merge_phyloseq(physeq_EC, sampledata)

#saveRDS(physeq_EC, file = "C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\physeq_EC.rds")

#physeq_EC <- readRDS(file = "C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\physeq_EC.rds")

And done! feel free to look at physeq_EC with the relative abundance graphs or beta diversity, but be warned. Some programs get ANGRY that this is data isnt count data (because it has decimals in it) (Specifically divnet - alpha diversity). Luckily there is a simple solution! just times everything by a large number!

This is as simple as…

EC_table <-10000000*EC_table
#run this command before you convert EC_table to a matrix