Importing data from DADA2 via Phyloseq one sample

By Angus Ball

Hi, by now you’ve processed your data through DADA2 and now its time for some data analysis!!! haha sike, we have to format it first, get rekt

These are the packages you’ll need

library(phyloseq)
library(ggplot2)
library(janitor)
library(dplyr)
library(tidyverse)
library(ggpubr)
library(ape)
library(splitstackshape)
library(Biostrings)

Fun fact: backups! hey i’ve noticed that you didn’t immediately create a backup and a working csv file, Shame on you! Its important to always create backup files you NEVER modify in excel when you get a new file you plan on modifying. It’d be a shame is you were to delete something you shouldn’t have and then restart from the beginning.

Firstly your data can exist in a variety of ways. If you have a txt file follow this (if you have a csv file skip it):

Now its time to import our data from DATA2, unfortunately we need to work with a .csv file and our data come off in text files Fun fact: .csv files are just basic excel files

To turn your text files into a .csv file open excel and then try to open your .txt file. Excel will go “hey…” and try to convert your txt file into an excel file. This is perfect! Our file is delimited (the autopicked option), via tab (the autopicked option after you click next once), so all you have to do is meander through the pop up and excel will do all the heavy lifting for you! Simply save your file as a .csv file where ever you’d like!

Now go off into to world and create a csv file of you 16S_R[name]_dada2_nochim_tax.txt file! Okay sweet, on the left should be a row of SV_0..Sv_10 ect, in the middle is you sample name with all the count data per species and on the right is tax.vector.

You’re 16S_R[name]_dada2_nochim_tax.csv file is now ready to be imported into R

Lets import this table

MetaG_1_total <- read_csv("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\ULTRA\\They call me... data\\Angus-16S\\16S-dada2_nochim_tax.csv")
Seq_1 <- read_csv("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\ULTRA\\They call me... data\\Angus-16S\\16S-sv_seqs.csv")

If I gave you an R object try

#MetaG_1 <- readRDS(filepath\\MetaG_1)
##Seq_1 <- readRDS(filepath\\Seq_1)

This is what the data looks like (note the tax.vector column at the very end)

head(MetaG_1_total)

head(Seq_1)

Congrats! We now have the data in R

Then we have to remove the SVs so there isnt multiple of the same value in the final table

#remove the sV value
MetaG_1_total <- select(MetaG_1_total, -(1))

First things first, we have to format the tables in a way that phyloseq is expecting them and unfortunately, this is a bit of a pain. Phyloseq expects a couple tables. First a count table (this is the MetaG table), A taxonomy table (This is tax.vector), and a sample data table (This ones based off of the group of data that exist in your sample)

Lets first start by creating the taxonomy table from the MetaG tables!

Tax<-select(MetaG_1_total, ("tax.vector")) 
#then we can remove tax.vector from the OTU table
MetaG_1_total <- select(MetaG_1_total, -"tax.vector")

Oh bother! Phyloseq expects the taxonomy to be separated by taxons i.e. one column for kingdom, one for family ect. But everything is under one column. This means we have to deliminate again, and looking at the how to column is created from DADA2, we need to deliminate by colon

cSplit deliminates within a column and then the rename functions make all the columns names phyloseq expects!

Tax <- cSplit(Tax, "tax.vector", ":")

Warning: 'as.is' should be specified by the caller; using TRUEWarning: 'as.is' should be specified by the caller; using TRUEWarning: 'as.is' should be specified by the caller; using TRUEWarning: 'as.is' should be specified by the caller; using TRUEWarning: 'as.is' should be specified by the caller; using TRUEWarning: 'as.is' should be specified by the caller; using TRUEWarning: 'as.is' should be specified by the caller; using TRUE

Tax <- dplyr::rename(Tax,
 Kingdom =  tax.vector_1,
 Phylum = tax.vector_2,
 Class = tax.vector_3,
 Order = tax.vector_4,
 Family = tax.vector_5,
 Genus = tax.vector_6,
 Species = tax.vector_7,
  
  )

#our datasets kinda wibbly wobbly because we've been adding things together, so sometimes if there isnt a designated taxon it writes NA and sometimes it write S__ or F__, lets fix that so it always just writes S__ ect
#first set to nas to zero
Tax[is.na(Tax)]<-0
#then find and replace with a similar command as above
Tax <- Tax %>%
  mutate(Kingdom = ifelse(Kingdom == 0, "k__",Kingdom))%>%
  mutate(Phylum = ifelse(Phylum == 0, "p__NA",Phylum))%>%
  mutate(Class = ifelse(Class == 0, "c__",Class))%>%
  mutate(Order = ifelse(Order == 0, "o__",Order))%>%
  mutate(Family = ifelse(Family == 0, "f__",Family))%>%
  mutate(Genus = ifelse(Genus == 0, "g__",Genus))%>%
  mutate(Species = ifelse(Species == 0, "s__",Species))

We have to take our tables (that exist as a dataframe) and convert them into a maxtrix, why? because the phyloseq class uses matrixes and not dataframes (in the example above we couldn’t feed our class two int values, it wants and int and a char!)

#now real stuff
#convert these bad boys to matrixes
Tax <- as.matrix(Tax)
MetaG <- as.matrix(MetaG_1_total)

You can check the class of the objects with a simple command (oh yes the tables have been objects all along, we’re creating objects within classes within classes! neat isnt it?)

#check
class(Tax)

[1] "matrix" "array" 

class(MetaG)

[1] "matrix" "array" 

Phyloseq expects the far left column to be called OTU, this is how to do it, I’m naming the OTUs 1 and 2 because of how phyloseq merges row names later on

#call row values OTU
rownames(MetaG) <- paste0("OTU", 1:nrow(MetaG))
rownames(Tax) <- paste0("OTU", 1:nrow(Tax))

Okay hopefully you know a bit about objects/classes, if not Funless fact: please go read the first couple chapters of a coding in java textbook Funner fact: data structure is neat, say you want to store a number, “1”, you can do so in java by telling it a datatype (int for integer), and then a value, for example int x=1; this tells the computer that x is equal to one. Great but int only works for numbers (and not even all of them ;), so what is you wanted to store a letter? well you use char (for character), so char y=“b”; this is saying that y is equal to the letter “b”. Obviously this quick example looses alot of the subtly and complexity of the system (re: go read the first couple chapters of a java textbook), but alas.

So now lets say you wanted a thing, lets call it a class, that contained a number and a character, you could create a class that used x and it used y to create an object. This is a poor explanation that has likely made several computer scientists sad but whatever, phyloseq is a class that creates phyloseq objects, each object having at a minimum an abundance table (metaG_) and a taxonomy table (Tax_). Excellent you are now equally confused as earlier but I feel like I’ve done my due dillegence, onward.

Now phyloseq wants us to turn these matrices into a different class created by phyloseq (not isnt a phyloseq object yet)

OTU = otu_table(MetaG, taxa_are_rows = TRUE)
TAX = tax_table(Tax)

Now we’ll create a phyloseq object with each OTU table and TAXonomy table

#combine that data
physeq = phyloseq(OTU, TAX)

You can see it worked by

physeq

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 47334 taxa and 79 samples ]
tax_table()   Taxonomy Table:    [ 47334 taxa by 7 taxonomic ranks ]

Excellent, almost done!

Now if you inspect the sample names they have a variety of meanings

physeq@otu_table

Things like “municipal compost”, “fertalized” ect, this is metadata that are in the sample names, but the phyloseq object doesn’t know about yet. Fun fact: if you want to know what this actually is or what the names mean ask to read my manuscript but we must import this metadata.

Key <- read_csv("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\ULTRA\\They call me... data\\Angus-16S\\ultra to categories key.csv")

Rows: 52 Columns: 5── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): Name, Amd, Conc, Fertalized, Inoculum
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

if you explore this document, you’ll realize it has all the samples in order they appear on column names within the MetaG table, but they are in the rows instead. The meta data in other columns

for the same reasons we converted to OTU earlier (phyloseq being fickle), we have to make the sample names be the left most row name

Key <- data.frame(Key[,-1], row.names=Key$Name)
sampledata = sample_data(Key)

Now we will add the key data

physeq_Key = merge_phyloseq(physeq, sampledata)

Hey remember how we imported the sequence information earlier but haven’t done anything with it? well lets add that information to the phyloseq object just like the key

We will pull the sequences from the file using biostrings, then the taxon names from the phyloseq object, then combine this information and add it to the phyloseq object

sequences <- Biostrings::DNAStringSet(Seq_1$Seq)
names(sequences) <- taxa_names(physeq_Key)
physeq_Key <- merge_phyloseq(physeq_Key, sequences)

Lets check if it worked

physeq_Key

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 47334 taxa and 51 samples ]
sample_data() Sample Data:       [ 51 samples by 4 sample variables ]
tax_table()   Taxonomy Table:    [ 47334 taxa by 7 taxonomic ranks ]
refseq()      DNAStringSet:      [ 47334 reference sequences ]

And hey refseq is there whoot! this isn’t terribly important for most of your analysis but it’ll be necessary for the odd thing

Congrats you are done (lie)!!! You now have a phyloseq object and can move on to transformations and data analysis

lets save it

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

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