NetCoMi

By Angus Ball

This tutorial is based off of the information provided here: https://github.com/stefpeschel/NetCoMi

install the packages! Note you may have to install git for this to work: https://git-scm.com/downloads

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()))

Load libraries

library(NetCoMi)

Loading required package: SpiecEasi

library(phyloseq)
library(dplyr)

Attaching package: ‘dplyr’

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

    filter, lag

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

    intersect, setdiff, setequal, union

library(ggplot2)

For my dataset i want to construct two networks and compare them one for the inoculcated trials and one for the uninoculated trials. I’m also going to remove all the primary biosolids samples so they don’t screw up my analysis (This is an experimental design choice). NetCoMi can take phyloseq objects as inputs so this is pretty easy.

physeq_Key <- readRDS(file = "C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\physeq_Key_unmerged.rds") #load physeq object
physeq_Key <- subset_samples(physeq_Key, Amd != "Primary sludge and biosolids mixture") #take the inoculated samples
#then ill conglomerate taxa to genus level
physeq_Key <- tax_glom(physeq_Key, taxrank = "Genus")

#then I'll remove all the taxa with zero abundance
otu_table_physeq_Key <- otu_table(physeq_Key) #steals the OTU table from the phyloseq object
otu_table_physeq_Key<-as.data.frame(otu_table_physeq_Key) #converts it to a dataframe

otu_table_physeq_Key$Sum <- rowSums(otu_table_physeq_Key[,])


otu_table_physeq_sparceless <- otu_table_physeq_Key %>% filter(Sum >20)#this removes any taxa that have less than 20 reads across all samples. yes i remove rare taxa, but this hurts our network analysis (which is already pretty inaccurate), so its a necessary evil. Why 20? because each network will be constructed with ~20 samples therefore this would imply that 50% of all samples in the network have this taxa (thats not how it works but thats how it statistically works)

otu_table_physeq_fin <- select(otu_table_physeq_sparceless, -"Sum") #removes the sum column
otu_table_physeq_fin <- as.matrix(otu_table_physeq_fin)
OTU = otu_table(otu_table_physeq_fin, taxa_are_rows = TRUE)

#extract taxa and sample data tables
TAX <- tax_table(physeq_Key)
SAM <- sample_data(physeq_Key)
#merge it
physeq_sparceless_woutPB <- merge_phyloseq(TAX, OTU, SAM)

physeq_inoculated <- subset_samples(physeq_sparceless_woutPB, Inoculum == "Inoculated") #take the inoculated samples


physeq_uninoculated <- subset_samples(physeq_sparceless_woutPB, Inoculum == "Uninoculated")

also note, these groups have to have the same amount of samples to be compared. unfortunately physeq_uninoculated has one less sample, than physeq_inoculated. therefore i will remove one sample from physeq_inoculated. The sample removed is NOT a random choice, but based on my experimental design

physeq_inoculated<- subset_samples(physeq_inoculated, sample_names(physeq_inoculated) != "Tailings-01-Unfertalized-Inoculated")

save these!

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

physeq_uninoculated <- readRDS("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\Robjects\\physeq_uninoculated.rds")
physeq_inoculated <- readRDS("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\Robjects\\physeq_inoculated.rds")

The rest I do in my linux box for SPEED, but the program seems to run pretty fast you might be able to do it on your own computer

netcomi requires each genus to have a unique identifier, which is fair enough, this function finds any genus with a non unique name (i.e all the unlabelled genuses at “g__“) and renames them to their lowest unique identifier. so g__test1 ->“test1” f__test2, g__ ->“test2(f)” (f for family) ect

physeq_inoculated <- renameTaxa(physeq_inoculated, 
                                  pat = "<name>", 
                                  substPat = "<name>_<subst_name>(<subst_R>)",
                                  numDupli = "Genus")
physeq_uninoculated <- renameTaxa(physeq_uninoculated, 
                                  pat = "<name>", 
                                  substPat = "<name>_<subst_name>(<subst_R>)",
                                  numDupli = "Genus")

Now we run the network analysis. I’ve decided to use SPRING it seems like the chillest option out of the bunch, but you’re welcome to read all the papers and make an informed choice too ;)

(An actual less sarcastic answer to the question of “why spring” is follows):

SPRING uses the mclr transformation, its like the clr transformation but without pseudocounts! note it doesn’t do this like the beta diversity packages we’ve been using, its doesn’t care about the wibbly wobbly actual effect adding a pseudocount does to your clr transformation and just knows what hip and happening (the previous comment should make about 0 sense, which is infact more sense than a detailed explanation of the mclr transformation). See it calculates the geometric means of a sample based on all your rows with values in them, and then it orders these by rank and sets anything that was a zero before the transformation to be at the bottom. So its not a “True” transformation in the sense that it treats everything equally, but it does the transformation on the data we need AND is rank preserving on the data we don’t care about. (This is very sexy ngl)

SPRING also incorporates measures for sparsity which alot of network analysis don’t explicitly do.

The Final reason I like SPRING is the fact its a conditional dependence method. This is beneficial because correlation based methods (like SparCC) often fail to differentiate between direct and indirect associations and there are ALOT of indirect associations to be picked up on in microbial network analysis, so mitigating this mitigates false conclusions

net_inoc_spring <- netConstruct(data = physeq_inoculated,#your physeq object, if you only want to make one network, do not include data2 = ..., otherwise this is where we tell the network analysis to make two networks based on these groups
                                data2 =physeq_uninoculated, 
                                taxRank = "Genus", #tells the network analysis to use genus level taxonomic information
                           verbose = 2, #how talkative the program is
                           filtTax = "highestFreq", #this tells the program to look at the taxa with the highest frequency
                           filtTaxPar = list(highestFreq = 85), #the highest 85 taxa
                           measure = "spring", #here is where you say which measure you want spring sparcc ect
                           measurePar = list(nlambda = 50, #these measures are based on what the paper used so thats why im using them
                                             rep.num = 50,
                                             Rmethod = "approx",
                                             ncores=10), #number of computer cores to use
                           normMethod = "none", #handled by spring itseld
                           zeroMethod = "none", #handled by spring
                           sparsMethod = "none", #spring returns a sparse network already so no extra sparsification necessary
                           dissFunc = "signed",
                           seed = 123456)

hey did you get a warning that looks like this? Warning message: In Matrix::nearPD(R, corr = TRUE) : ‘nearPD()’ did not converge in 100 iterations That means your data was too scarce! The easiest way to do this is just decrease the amount of taxa you look at I.e. from 75->65 under filtTaxPar = list(highestFreq = x). This will functionally only ttake the highest abundance taxa and there for solve the sparcity issue

saveRDS(net_inoc_spring, "net_inoc_spring.rds")

Then we can analyze the networks!

props_inoc <- netAnalyze(net_inoc_spring, #your network
                           centrLCC = FALSE, #i dont understand
                           avDissIgnoreInf = TRUE, 
                           sPathNorm = FALSE,
                           clustMethod = "cluster_fast_greedy",
                           hubPar = c("degree", "eigenvector"),
                           hubQuant = 0.95,
                           lnormFit = TRUE,
                           normDeg = FALSE,
                           normBetw = FALSE,
                           normClose = FALSE,
                           normEigen = FALSE,
                           verbose = 2)

These should and can change based on your questions, but honestly I dont understand it very well, I’m using the tutorial specs

then we can plot it out

This will plot out just the networks as they want to exist individually, but this will make comparison hard!

props_inoc <- readRDS("C:\\Users\\angus\\OneDrive - UNBC\\Angus Ball\\Lab work\\Bioinformatics\\ULTRA\\Robjects\\props_inoc.rds")

plot(props_inoc, 
     layout = "spring", #shape, alternative option is circle
     sameLayout = FALSE, #same layout of both graphs?
     nodeColor = "cluster", #how to color nodes
     repulsion = 0.7, #how far apart are nodes
     nodeSize = "mclr", #node size is based off of clr transform size
     labelScale = F, #don't scale names with node size
     shortenLabels = "none",#if you want to shorten the label names, othe roption is "intelligent"
     cexNodes = 2, #size of nodes
     cexLabels = .5,#size of names
     cexHubLabels = .75,#size of hub names
     cexTitle = 2,#size of title
     groupNames = c("Inoculated", "Uninoculated"), #group names
     hubBorderCol  = "gray40")

legend(0.7, 1.1, cex = 2.2, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#009900","red"), 
       bty = "n", horiz = TRUE)

We can instead force the networks to look the same in OTU order for easy comparison

plot(props_inoc, 
     layout = "spring", #shape, alternative option is circle
     sameLayout = TRUE, #same layout of both graphs?
     layoutGroup = "union", #combination of both layouts
     rmSingles = "inboth", #remove single taxa no associated with other taxa in both groups
     nodeColor = "cluster", #how to color nodes
     repulsion = 0.7, #how far apart are nodes
     nodeSize = "mclr", #node size is based off of clr transform size
     labelScale = F, #don't scale names with node size
     shortenLabels = "none",#if you want to shorten the label names, othe roption is "intelligent"
     cexNodes = 2, #size of nodes
     cexLabels = .5,#size of names
     cexHubLabels = .75,#size of hub names
     cexTitle = 2,#size of title
     groupNames = c("Inoculated", "Uninoculated"), #group names
     hubBorderCol  = "gray40")

legend(0.5, 1.1, cex = 2.2, title = "estimated association:",
       legend = c("+","-"), lty = 1, lwd = 3, col = c("#009900","red"), 
       bty = "n", horiz = TRUE)

For some reason ggsave doesn’t seem to want to work with these plots. But hand exporting it with the export button above the graph works so uhh just use that

explanations of the given variables is given at https://rdrr.io/github/stefpeschel/NetCoMi/man/plot.microNetProps.html

we will then use statistical tests to compare the structure of these two networks. NOTE: the p value displayed by this is NOT comparing the two networks. Its suggesting that the absolute difference between the networks is not equal to zero. take this as you will.

Note this process will take 5ever. estimated time 8-10hr

comp_inoc_orig <- netCompare(props_inoc, permTest = TRUE, nPerm = 1000, 
                          storeAssoPerm = TRUE,
                          fileStoreAssoPerm = "assoPerm_comp",
                          storeCountsPerm = FALSE,
                          cores = 1, #got an error for using multiple cores so only using one :(
                          seed = 123456,
                          adjust = "BH",
                          verbose = TRUE)
#we can explore the output in
summary(comp_inoc_orig)

Then we can see if there are differential associations within our network

diff_inoc <- diffnet(net_inoc_spring,
                       diffMethod = "fisherTest", 
                       adjust = "lfdr")
saveRDS(diff_inoc, "diff_inoc.rds")
#I didn't have any but if you do...
plot(diff_inoc, 
     cexNodes = 0.8, 
     cexLegend = 3,
     cexTitle = 4,
     mar = c(2,2,8,5),
     legendGroupnames = c("Inoculated", "Uninoculated"),
     legendPos = c(0.7,1.6))

props_inoc_pears <- netAnalyze(net_inoc_spring, 
                                 clustMethod = "cluster_fast_greedy",
                                 weightDeg = TRUE,
                                 normDeg = FALSE,
                                 gcmHeat = FALSE)
# Identify the differentially associated OTUs
diffmat_sums <- rowSums(diff_inoc$diffAdjustMat)
diff_asso_names <- names(diffmat_sums[diffmat_sums > 0])
plot(props_inoc_spring, 
     nodeFilter = "names",
     nodeFilterPar = diff_asso_names,
     nodeColor = "gray",
     highlightHubs = FALSE,
     sameLayout = TRUE, 
     layoutGroup = "union",
     rmSingles = FALSE, 
     nodeSize = "clr",
     edgeTranspHigh = 20,
     labelScale = FALSE,
     cexNodes = 1.5, 
     cexLabels = 3,
     cexTitle = 3.8,
     groupNames = c("Inoculated", "Uninoculated"),
     hubBorderCol  = "gray40")

legend(-0.15,-0.7, title = "estimated correlation:", legend = c("+","-"), 
       col = c("#009900","red"), inset = 0.05, cex = 4, lty = 1, lwd = 4, 
       bty = "n", horiz = TRUE)

And you’re done! FREE, to live another day