ANCOMBC2_two

By Angus Ball

Had enough of me by now? you can also read the devs tutorial here: https://www.bioconductor.org/packages/release/bioc/vignettes/ANCOMBC/inst/doc/ANCOMBC.html

Welcome to the fray! ANCOM-BC-2 is great and has alot of cool features, read about em (and understand them >:[! (These are the same as the ones on the git):

So what’s ANCOM-BC-2’s deal?

-Multiple pairwise comparisons!! This may not sound like a big deal but this is the first program I’ve encountered that promises it’ll do all the comparisons you care about in groups of 3+ rather than just having a reference taxon and you have to fiddle around with it. (it’ll also do binary comparisons)

-Pattern analysis over ordered study groups (effectively pairwise comparisons but with a fuzzer image)

-ANCOM-BC-2 will also account for taxon-specific bias. Remember how choosing your PCR primers was so important? because they will preferentially bias certain taxons? ANCOMBC2 promises to account for this in your data.

-(This one is really cool) ANCOMBC2 will perform sensitivity analysis based on choice of pseudo-count. See pseudo-count can change your results (specifically regarding rare taxa as they are more biased by the psuedo-count), leading to an inflated FDR. ANCOMBC2 will run the analysis with multiple pseudo-counts and see how they effect the results and filter out taxa that are significantly affected by pseudo-count.

-ANCOMBC2 knows about compositionality and how scary it is. See we use relative abundance because microbiome data is compositional, but a change in the absolute abundance of a single taxon can alter the relative abundances of all taxa (This is shockingly not great). Therefore, when we hypothesis test comparing relative abundances is not equivalent to absolute abundances. Relative abundance/ sampling fraction can be affected by the microbial load in the system and library size. Library size is simple too small a library size you don’t hit the important parts of your rarefaction curve and you miss out on species. Microbial load is a bit more complex (and explained well in the ANCOMBC paper). But in short two ecosystems/samples can have the same relative abundance and library sizes of a taxon, but their absolute population size, or microbial load, is different. If you don’t account for this then you’d say the two ecosystems are the same (wrong). ANCOMBC accounts for microbial load by introducing a sample-specific offset term to bias correct for the microbial load and constitutionality of microbiome data.

-Preforms worse when sample size is small (n = 5), and runs perfectly when n >= 10

-Accounts for structural zeros. Structural zeros are defined as taxons present in one sample type and not present in another (i.e. imagine comparing samples from a jungle to a forest, there’ll just be some species that straight up don’t exist in one environment or the other and this is expect). You can’t really do stats on these structural zeros because its just not there in one sample. So we’ll have to examine this differently (and subjectively). Please give Naught all zeros in sequence count data are the same a skim it has some stuff to say about zeros (if you couldn’t guess by the title). Structural zeros are called biological zeros in this paper. Now this being said the way these structural zeros are defined is a bit suspect. See if your library size in some samples is too small it won’t pick up rare taxa at all! so these rare taxa might be erroneously picked up as structural zeros (where they should be called sampling zeros). However, in my humble opinion using the same sequencing tech and primers should mitigate the issues with sampling zeros (and technical zeros) as you’d expect equal various among samples. ANCOMBC skirts some of these issues by first filter taxa to have to exist in a certain amount of samples (10% by default) to get around rare taxa getting picked up just because they are rare.

Cool, so how do we get ANCOMBC running?, well This package is new as hell, and the biocmanaer version isn’t up to date, so we’ll be downloading straight off of the github

#Heres how to install it
devtools::install_github("FrederickHuangLin/ANCOMBC", branch = "bugfix")

it uhh takes a second

You’ll also need these packages (Which you probably don’t have)

#Also download these bad boys
devtools::install_github("vmikk/metagMisc")
devtools::install_github("zdk123/SpiecEasi")
devtools::install_github("GraceYoon/SPRING")

devtools::install_github("stefpeschel/NetCoMi", 
                         dependencies = c("Depends", "Imports", "LinkingTo"),
                         repos = c("https://cloud.r-project.org/",
                                   BiocManager::repositories()))

#if NetCoMi doesn't work, you may need to install git on your computer https://git-scm.com/downloads

In the end this is the libraries you’ll need (You probably have most of em already + there aren’t special download instructions, use biocmanager for phyloseq, and speedyseq)

Load your phylosq object!

physeq_Key <- readRDS("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\Lisas data\\R objects\\physeq_Key.rds")

Now I’m going to run ANCOMBC2 on a binary comparison (i.e. two options), if you have 3+ comparisons please see the other tutorial!

This data has samples from a treated area and control area, we only need year 1 samples for this analysis so lets subset them now

physeq_yr1 <- subset_samples(physeq_Key, Years == "1")

We are also going to run this analysis on a Family level. I recommend running this analysis at a family and genus level (Especially if you are working with fungi, since then you can still use Fungaltraits to look at species comparisons). But you’re welcome to use any taxon level you like!

physeq_glom <- speedyseq::tax_glom(physeq_yr1, #physeq object
                                  taxrank = "Family", #taxa you want to conglomerate to
                                  NArm = FALSE) #dont remove NA taxa

Then we have some data management to do!

#First remove any taxa with a sum of sequence counts equal to zero (You might have this if you have removed specific samples from your analysis)
physeq_glompruned <- prune_taxa(taxa_sums(physeq_glom) >0, physeq_glom)

#Then we are going to rename the taxa, this is why we downloaded netcomi, so especially in fungal taxonomy there are taxa with the same name even at a family level! its all very confusing and annoying. This next function efffectively just renames all the duplicates to taxa_1 taxa_2 ect so we can determine what species actually are differentially abundant 

physeq_Key_renamed <- renameTaxa(physeq_glompruned, 
                                 pat = "<name>", 
                                 substPat = "<name>_<subst_name>(<subst_R>)",
                                 numDupli = "Family") #Change family to taxon of interest

Columns 6 contain NAs only and are ignored.Columns 7 contain NAs only and are ignored.

ANCOM will remove taxa that arent present (or prevalent) in enough samples so theres some statistical trust in the difference (i.e. there could be a massive difference between the relative abundance of a taxa in treated compared to control, but if its only because of one sample does it really mean anything?)

However this seems to only work for the actual statistical test not the pretest structural zero measurement! so we’ll just help ancombc by running the prevelance filter first with yet another package

physeq_prev <- metagMisc::phyloseq_filter_prevalence(physeq_Key_renamed, prev.trh = 0.1) #this makes sure to only include samples with a prevelance over 10% (0.1) (i.e taxa must be present in 10% of samples). 

PS. do you NEED to rename your physeq object after each command? of course not! hell you can run all these things under one command with enough %>% pipes, but since I usually make these code decks first troubleshooting comes up ALOT, so having seperate names lets me easily know where things go wrong.

We’ll also extract the taxa table from the phyloseq object, this becomes important later on

family_taxa <- as(phyloseq::tax_table(physeq_prev), "matrix")
family_taxa <- as.data.frame(family_taxa)

For whatever reason you cant feed ancombc a phyloseq object so uhh make this dumb rando object with mia

year1 = mia::makeTreeSummarizedExperimentFromPhyloseq(physeq_prev)

Most of my names will be year1_something, thats just because of how I did my analysis and what this data says! feel free to change it at your hearts content

Now… run ancombc

output_year1 = ancombc2(data = year1, #data
                  assay_name = "counts", #arbitrary name
                  tax_level = "Family", #level of comparison, like "Family", NULL is base level (ASV)
                  fix_formula = "Treat", #covariate that controls similar samples
                  rand_formula = NULL, #do microbial abundances have reandom effects in data? lets hope not...
                  p_adj_method = "holm", # p value adjust method
                  pseudo_sens = TRUE, #ANCOM-BC2 doesn't use a zero count, but it measures the effect a pseudocount will have on your data, how nice!
                  prv_cut = 0.1, #means taxa need to have non-zero counts in over 10% of samples, BUT doesn't work for structural zeros, therefore we have the change above too
                  lib_cut = 0, #removes samples with less reads than this number
                  s0_perc = 0.05,#this protets you from rare taxa being signifigantly different just because of small standard errors and not differential abundance
                  group = "Treat", #comparison group, must be in covariate list
                  struc_zero = TRUE, #structural zeros mean the sample is fully absent in one sample type compared to the other 
                  #i.e. if you were comparing a rainforest and a desert there would be a bunch of structural zeros clouding the analysis
                  #by saying yes here the program checks for structural zeros and if it finds any it automatically decides they are differentially abundant (bc they are)
                  neg_lb = FALSE,#makes the definition of struct_zero a bit wider, recommended if sample size is over 30 per group 
                  alpha = 0.05, #hopefully you know what this one is ;)
                  n_cl = 1, #how many nodes/clusters to be used, if youre on a windows or mac just say 1
                  verbose = TRUE,#talkative >:]
                  global = FALSE, #perform global test y/n
                  pairwise = FALSE, #hey want to do pairwise tests for multiple group comparisons make this true! 
                  dunnet = FALSE, #dunnets test involves the use of a control group so weird dont use?
                  trend = FALSE, #another random test
                  )

Checking the input data type ...
The input data is of type: TreeSummarizedExperiment
PASS
Checking the sample metadata ...
The specified variables in the formula: Treat
The available variables in the sample metadata: LW., block.site.ID..true.replicates., Treated.or.Control, Plot....repeated.measures.within.each.replicate., Years.post.treatment..controls.are.not.treated.but.coordinate.with.blocks.of.those.treatment.years., Total.Mass.g., Course.Fragment.Mass.g., percent.coarse.fragment, Colour, Texturing, BlockTreat, Years, YearsTreat, Block, Treat, TreatBlock
PASS
Checking other arguments ...
The number of groups of interest is: 2
Warning: The group variable has < 3 categories 
The multi-group comparisons (global/pairwise/dunnet/trend) will be deactivatedThe sample size per group is: control = 22, treated = 23
PASS
Obtaining initial estimates ...
Estimating sample-specific biases ...
Loading required package: foreach

Attaching package: ‘foreach’

The following objects are masked from ‘package:purrr’:

    accumulate, when

Loading required package: rngtools
Sensitivity analysis for pseudo-count addition to 0s: ...
ANCOM-BC2 primary results ...

Extract the structural zero data… (go read the ancombc paper if you don’t know what it is)

#which families are structural zeros?
tab_zero_year1 = output_year1$zero_ind #extract the structural zero data

tab_zero_year1 %>%
  datatable(caption = "The detection of structural zeros") #display it nice

#this only has the family label on it, lets add some more information from the family taxa dataframe we extracted early

tab_zero_year1_bound <- cbind(tab_zero_year1, family_taxa) #this just matches columns with the same information in them

#remove random columns
tab_zero_year1_clean <- select(tab_zero_year1_bound, -"taxon", -"Genus", -"Species")


#now lets make a table of taxons that are only structural zeros

#structural zeros in mineral layer
tab_zero_year1_Sig_control <- dplyr::filter(tab_zero_year1_clean, `structural_zero (Treat = control)` == TRUE )  #im filtering them based on if they were from the control or treated fractions

#Structural zeros in organic layer
tab_zero_year1_Sig_treated <- dplyr::filter(tab_zero_year1_clean, `structural_zero (Treat = treated)` == TRUE ) 

# you can now use these tables and maybe make some kind of graphic to display which of these taxa are structural zeros? 

Now lets extract the meat and potatoes of our data

res_prim_year1 = output_year1$res

Now lets make a nice graph/readable tables from it, because if you look at res_prim_year1 oh boy is there alot of random stuff.

df_sig_y1 = res_prim_year1 %>%
  dplyr::select(taxon, contains("Treat")) %>% #First we want to select samples that are for the treat comparison (this is only important if you have multiple comparisons at play)
  dplyr::filter(diff_Treattreated == "TRUE") %>% #then select taxa that passed the differential abundance test diff_
  dplyr::filter(passed_ss_Treattreated == "TRUE") #then select taxa that passed the sensitivity analysis (again go read the ancombc paper)

#a reminder for a taxa to be included in your results it must pass BOTH the sensitivity analysis and differential abundnace tests


#then we can add some useful information from our family taxa table to this datatable
df_sig_y1 <- merge(family_taxa, df_sig_y1, by.y = "taxon", by.x = "Family")


df_sig_y1 <- df_sig_y1 %>%
  dplyr::arrange(desc(lfc_Treattreated)) %>%
  dplyr::mutate(direct = ifelse(lfc_Treattreated > 0, "Positive LFC", "Negative LFC"))
df_sig_y1

So theres alot going on here… lfc_ represents the lofc fold change from the alternate comparison group lfc_(intercept) can be ignored, its a stats thing that has no bearing on your conclusions

se_ is standard error

W_ and p_ are our test statistics and p value, but only worry about q_ which is our family-wise error rate adjusted p value.

Now lets make a pretty graph!

fig_y1 = df_sig_y1 %>%
  dplyr::mutate(Family = fct_reorder(Family, dplyr::desc(lfc_Treattreated))) %>% #reorders taxa log fold change in a descending pattern 
  ggplot(aes(x = Family, y = lfc_Treattreated, fill = direct)) + 
  geom_bar(stat = "identity", width = 0.7, color = "black", 
           position = position_dodge(width = 0.4)) + #adds the lfc change itself
  geom_errorbar(aes(ymin = lfc_Treattreated - se_Treattreated, ymax = lfc_Treattreated + se_Treattreated),
                width = 0.2, position = position_dodge(0.05), color = "black") + #adds our error bars
  labs(x = NULL, y = "Log fold change", 
       title = "Log fold changes of fungi families after treatment in one year old plots") + #er title...
  scale_fill_discrete(name = NULL) +
  scale_color_discrete(name = NULL) +
  theme_bw() + #ggplot fun theme stuff, feel free to play around
  theme(plot.title = element_text(hjust = 0.5),
        panel.grid.minor.y = element_blank(),
        axis.text.x = element_text(angle = 60, hjust = 1
                                   )) # This makes the names at the bottom of the graph all cool and sideways looking
fig_y1

And you’re done! In the next tutorial I will show the difference between multiple group comparisons