Aldex2 differential abundance

By Angus Ball

library’s

library(phyloseq)
library(ALDEx2)
library(microbiomeMarker)

Differential abundance time! First load the data

filepathing <- "C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\Kenzies Data\\physeq_Key.rds" #add your file location here
physeq_Key <- readRDS(file = filepathing)

Loading required package: phyloseq
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Go read this massive scary as hell document, or at least skim the interesting parts???? https://www.bioconductor.org/packages/devel/bioc/vignettes/ALDEx2/inst/doc/ALDEx2_vignette.html

Note Aldex the package doesn’t have an option for adding psuedocounts during the clr transformation. Therefore we are going to apply the psuedocount then run aldex without the transformation.

But HEY Angus you said the robust clr transformation was better! Why are we using pseudocounts here? The answer is because all the papers that I’ve read that compare differential abundance methodologies all use a pseudocount. Therefore, the philosophy is using something that could be worse than the best option, but has been proven better than the most alternatives.

Now, there are many ways to add a pseudocount, but hopefully you’ve come to appreciate my hamfisted approaches so…

physeq_psuedo <- physeq_Key #copy of physeq_Key
otu_table_pseudo <- otu_table(physeq_psuedo)
otu_table_pseudo<-as.data.frame(otu_table_pseudo)
otu_table_pseudo<-otu_table_pseudo +1 # literally justs adds one to every value in otu_table

# Then we need to take our sample data and convert it to a vector
sd = data.frame(sample_data(physeq_Key)) #takes the sample data out
vect <- as.vector(sd[,"Layer"]) #switch layer with whatever column you're interested in

Fun fact, the aldex package proper doesn’t use phyloseq objects, just an honest to god OTU table, and the vector. BTW in our OTU table all the species are called by OTU, i.e OTU1 OTU2 ect. you can change this here and now on the otu table or be arsed with it later (what Ill show)

Time to run aldex

physeq.aldex<-aldex(otu_table_pseudo, # OTU table
                    vect, #sample data in a vector
                    mc.samples = 128, #how many monte carlo simulations to run, 128 is recommended
                    denom = "all", #which OTU's to use, i.e. all of them, but there may be cases where you don't want to use all of em.
                    verbose = TRUE, #tells you progress
                    test = "t", #t test because I'm looking at a binary comparision here, use "kw" if you have multiple comparisons 
                    effect = TRUE) #determines if you want to calculate abundances and effect size

See https://www.bioconductor.org/packages/devel/bioc/vignettes/ALDEx2/inst/doc/ALDEx2_vignette.html#722_aldexttest(x)

for explanations on variable names

#extract taxa information
taxa <- tax_table(physeq_Key)
taxa <- as.data.frame(taxa)

aldex.treat.bound <- cbind(physeq.aldex, taxa)#bind taxa informaton to the aldex object

aldex_sig_treat <- aldex.treat.bound[ aldex.treat.bound$"wi.eBH" < 0.05, ] #extract significant values
aldex_sig_treat #there wasnt any signifigant values :(  better luck next time!

Now we did this at a species level which typically means there are just so many comparisons that nothing will ever be significant. use taxa_glom to conglomerate all the taxa to different levels (i.e. genus and class ect) and see if something pops up!

#e.g.
physeq_genus <- tax_glom(physeq_Key, taxrank="Genus")

PS. theres a method to do linear combinations of variable comparisons but its beyond the scope of my time right now