1. Load data and construct network

Set up phyloseq

tree_file = '~/Dropbox/RoL/pitcher-16s/qiime2-stuff/rooted-tree/tree.nwk'
df = '~/Dropbox/RoL/pitcher-16s/site-data-bioclim-pca.txt'

ps = merge_phyloseq(import_biom('~/Dropbox/RoL/pitcher-16s/table-with-taxonomy-converted.biom', parseFunction = parse_taxonomy_greengenes), 
                    sample_data(read.csv(df, sep = '\t', row.names = 1)), 
                    read_tree(tree_file))
tax_table(ps) = tax_table(ps)[,1:7]
colnames(tax_table(ps)) = c('Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species')

# filter out taxa unassigned at Kingdom level
ps = ps %>% subset_taxa( Kingdom != "Unassigned")

# drop florida because it's a different subspecies
ps = ps %>% subset_samples(site_id != 'FLK')

# drop sample with ph = 1
# ps = ps %>% subset_samples(ph>1)

# agglomerate into OTUs
ps_OTU = ps %>% tip_glom(h = 0.05, hcfun=hclust)

# add OTUID column
tax_table(ps_OTU) = cbind(tax_table(ps_OTU), OTUID=paste0("OTU_",1:length(taxa_names(ps_OTU))) )

ps6 = filter_taxa(ps_OTU, function(x) sum(x >= 1) >= 6, TRUE) %>% rarefy_even_depth(rngseed = 777)
flist = ps6 %>% transform_sample_counts(function(x) x / sum(x)) %>% filter_taxa(function(x) sum(x>= 0.01) > 0)
ps6Top = prune_taxa(flist, ps6)

Make networks

# se.ps6Top <- spiec.easi(ps6Top, method='mb', lambda.min.ratio=1e-1,
#                            nlambda=30, pulsar.params=list(rep.num=99, ncores=4))
# saveRDS(se.ps6Top, "speiec-easi-OTU-network.rds")
se.ps6Top = readRDS("speiec-easi-OTU-network.rds")
sc.ps6Net =make_sparcc_net(ps6Top, 0.25)
# plot the two networks
se.ps6Net <- adj2igraph(getRefit(se.ps6Top), vertex.attr=list(name=as.vector(taxa_names(ps6Top))))

plot_network(se.ps6Net, ps6Top, type='taxa')
plot_network(sc.ps6Net, ps6Top, type='taxa')

Intersection of networks

sc.se = intersection(sc.ps6Net, se.ps6Net)
# taxa=ps6Top %>% subset_taxa(taxa_names(ps6Top) %in% names(V(sc.se))) %>% tax_table(.) 
# vertex_attr(sc.se, "name") <- as.vector(taxa[,'Genus'])

# remove some verticies to make it easier to see
verticesToRemove <- V(sc.se)[igraph::degree(sc.se) < 2]
# These edges are removed from the graph
sc.se <- delete.vertices(sc.se, verticesToRemove)

plot_network(sc.se, ps6Top)

#vertex_attr(sc.ps6Net, "name") <- as.vector(tax_table(ps6Top)[,'Genus'])

Network parameters

network_results = function(net) {
  E(net)$weight = NA
  #  modularity
  fc = cluster_fast_greedy(net,weights =NULL)
  modularity = modularity(net,membership(fc))
  # graph properties
  results = data.frame(clustering.coefficient=transitivity(net), # clustering.coefficient
                modularity=modularity,# modularity
                mean.degree= mean(igraph::degree(net)),# mean degree
                size=length(E(net)),# size - number of edges
                order=length(V(net)),# order - number of vertices
                mean.distance=mean_distance(net),#average path length
                norm.degree=mean(igraph::degree(net)/length(V(net))),#normalized degree
                betweenness.centrality=centralization.betweenness(net)$centralization #betweenness centrality
                )
  return(results)
}
network_results(sc.se)

Is the network scale-free? Yes. A significant p-value indicates that the network does NOT follow a power law.

fit_power_law(igraph::degree(sc.se))$KS.p
## [1] 0.9990253

2. Associate modules with variables

Option A: Fast greedy plus merging with WGCNA

Build network

# https://stackoverflow.com/questions/51271730/how-to-remove-small-communities-using-igraph-in-r
SizeToKeep = 5
fc = cluster_fast_greedy(sc.se)

N.Mods = length(which(table(fc$membership) >= SizeToKeep))
# change color hex code to name so I can see what module is identified
my_colors = RColorBrewer::brewer.pal(N.Mods,"Paired")
# set small module nodes to "grey"
color_names = data.frame(color.name = sapply(my_colors, color.id))$color.name
color_names = c(color_names, rep("grey", length(which(table(fc$membership) < SizeToKeep))))

comps = membership(fc)
V(sc.se)$color = color_names[comps]
plot(sc.se, vertex.label=NA,vertex.size=5)

# prep data
datExpr = data.frame(t(otu_table(ps6Top %>% 
                                   transform_sample_counts(function(x) x / sum(x)) %>%
                                   subset_taxa(taxa_names(ps6Top) %in% names(V(sc.se)))
                                   )))

TOM = TOMsimilarity(as_adjacency_matrix(sc.se, sparse = F))
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
dissTOM = 1-TOM
# Call the hierarchical clustering function
geneTree = hclust(as.dist(dissTOM), method = "average")

# Get Module Eigengenes per cluster
MEs0 <- moduleEigengenes(datExpr, colors = V(sc.se)$color)$eigengenes
# Reorder modules so similar modules are next to each other
MEs0 <- orderMEs(MEs0)
module_order = names(MEs0) %>% gsub("ME","", .)

# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs0)
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")

Plot module similarity

# Plot the result
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
# choose threshold: corresponds to correlation of 1-MEDissThres
MEDissThres = 0.4
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")

# Call an automatic merging function
merge = mergeCloseModules(datExpr, V(sc.se)$color, cutHeight = MEDissThres, verbose = 3)
##  mergeCloseModules: Merging modules whose distance is less than 0.4
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 9 module eigengenes in given set.
##    Calculating new MEs...
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 9 module eigengenes in given set.
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs

plotDendroAndColors(geneTree, cbind( V(sc.se)$color, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

Network with merged modules

V(sc.se)$color = mergedColors
plot(sc.se, vertex.label=NA,vertex.size=5)

Module-site relationships

MEs0 <- mergedMEs
# Reorder modules so similar modules are next to each other
MEs0 <- orderMEs(MEs0)
module_order = names(MEs0) %>% gsub("ME","", .)

MEs0$treatment = row.names(MEs0)

# arrange by a variable here to get a sense of geography
site_order = rownames(arrange(data.frame(sample_data(ps6)), bioclim_pca1))

# tidy & plot data
mME = MEs0 %>%
  pivot_longer(-treatment) %>%
  mutate(
    name = gsub("ME", "", name),
    name = factor(name, levels = module_order),
    treatment = factor(treatment, levels=site_order)
  )
mME %>% ggplot(., aes(x=treatment, y=name, fill=value)) +
  geom_tile() +
  theme_bw() +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(-1,1)) +
  theme(axis.text.x = element_text(angle=90)) +
  labs(title = "Module-site relationships sorted by bioclim_pca1", y = "Modules", fill="corr")

Relatedness within modules (SES MPD)

# add a "module" column to the tax table
ps.network = ps6Top %>% subset_taxa(taxa_names(ps6Top) %in% names(V(sc.se)))
tax_table(ps.network) = cbind(tax_table(ps.network), Module=mergedColors)
module_order = names(mergedMEs) %>% gsub("ME","", .)
# subset the tax table to the Module variable
mods = data.frame(tax_table(ps.network)[,"Module"])
# manually add module ID to taxa
mods = mods %>% mutate(mod1 = ifelse(Module==module_order[1],1,0),
                mod2 = ifelse(Module==module_order[2],1,0),
                mod3 = ifelse(Module==module_order[3],1,0),
                mod4 = ifelse(Module==module_order[4],1,0),
                mod5 = ifelse(Module==module_order[5],1,0),
                mod6 = ifelse(Module==module_order[6],1,0),
                mod7 = ifelse(Module==module_order[7],1,0),
                mod8 = ifelse(Module==module_order[8],1,0),
                mod9 = ifelse(Module==module_order[9],1,0)) %>%
  select(-Module)
names(mods) = module_order
# set up data 
rooted = phy_tree(ps.network)
# otu = data.frame(otu_table(ps.network))
otu = mods
match <- match.phylo.comm(phy = rooted, comm = t(otu))

# SES MPD of modules
mod_SESMPD = ses.mpd(match$comm, cophenetic(match$phy), null.model = "taxa.labels")
kable(mod_SESMPD %>% select(c(ntaxa,mpd.obs.z,mpd.obs.p)))
ntaxa mpd.obs.z mpd.obs.p
lightpink2 10 0.0297792 0.501
thistle3 9 -0.3550796 0.349
firebrick2 13 -0.5960268 0.276
darkorange1 8 -2.5098436 0.009
forestgreen 5 -1.5470449 0.062
darkseagreen3 11 1.3614358 0.919
sandybrown 19 0.0888211 0.518
dodgerblue3 7 -2.6901916 0.007
lightskyblue2 18 -1.1173234 0.138

Correlation plot

library(corrplot); library(psych)
## corrplot 0.92 loaded
dat = data.frame(sample_data(ps6Top)) %>%
  cbind(mergedMEs) %>%
  na.omit(.)

dat2 = data.frame(sample_data(ps6Top)) %>%
  select(-c(subspec,site_id,Description,site_id.1))

ct = corr.test(dat2, mergedMEs)

corrplot(ct$r, p.mat = ct$p.adj, sig.level = 0.1 , insig = "blank")

Regression of modules with environmental data

Latitude

ggplot(dat, aes(latitude, MEfirebrick2)) +geom_point()

mod1 = lm(MEfirebrick2~latitude, data=dat)
anova(mod1)
summary(mod1)
## 
## Call:
## lm(formula = MEfirebrick2 ~ latitude, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14464 -0.05201 -0.01827  0.02477  0.28768 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.371126   0.068298   5.434  4.4e-07 ***
## latitude    -0.008355   0.001526  -5.475  3.7e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.087 on 93 degrees of freedom
## Multiple R-squared:  0.2437, Adjusted R-squared:  0.2356 
## F-statistic: 29.97 on 1 and 93 DF,  p-value: 3.703e-07

Pitcher size and prey

ggplot(dat, aes(inverts, MEdarkorange1)) +geom_point()

mod1 = lm(MEdarkorange1~inverts, data=dat)
anova(mod1)
summary(mod1)
## 
## Call:
## lm(formula = MEdarkorange1 ~ inverts, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.02161 -0.01999 -0.01882 -0.00748  0.21243 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0090783  0.0065672  -1.382    0.170
## inverts      0.0002657  0.0023052   0.115    0.908
## 
## Residual standard error: 0.04716 on 93 degrees of freedom
## Multiple R-squared:  0.0001428,  Adjusted R-squared:  -0.01061 
## F-statistic: 0.01328 on 1 and 93 DF,  p-value: 0.9085

OTUs in the modules

latitudeOTUs = ps.network %>% 
  transform_sample_counts(function(x) x / sum(x)) %>%
  subset_taxa(Module=="firebrick2") 
kable(tax_table(latitudeOTUs))
Kingdom Phylum Class Order Family Genus Species OTUID Module
90086558ccc079f6a635830df3a4ebb9 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter NA OTU_64 firebrick2
b1f0f8eff3d0ea0d7a19fc6178e9ebeb Bacteria Bacteroidota Bacteroidia Flavobacteriales Weeksellaceae NA NA OTU_214 firebrick2
df6b18712acf2c44bc154778f6b3ef09 Bacteria Acidobacteriota Acidobacteriae Acidobacteriales Acidobacteriaceae_(Subgroup_1) Granulicella NA OTU_650 firebrick2
572e3dfbe64c22ad15717c464f1534f0 Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae NA NA OTU_841 firebrick2
31b69392eb47ad5bf4cc470deeb175da Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae NA NA OTU_989 firebrick2
b71ab0ec9fae32ebaf2d8ec2ccc312b9 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Novosphingobium NA OTU_1001 firebrick2
b868361da89d5b4a8941afce5c860246 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Chromobacteriaceae Aquitalea NA OTU_1284 firebrick2
a93cbf059fd5f2636658ba59355e73d8 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Comamonadaceae NA NA OTU_1354 firebrick2
2010445510c2a3c9db056f3e365441df Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Comamonadaceae NA NA OTU_1360 firebrick2
0a17be34f84b48eecd0b82aec06a0929 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Comamonadaceae NA NA OTU_1398 firebrick2
882330144d2c06e8c47fc6a4f8988141 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1434 firebrick2
72a610c3ddd4199e03399a0ae4e4c7a1 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1480 firebrick2
504b1ad547b40a63f01095f8fbb593a4 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1481 firebrick2 
#write.table(tax_table(latitudeOTUs), "latitudeOTUs2.csv", sep = ',')
invertOTUs = ps.network %>% 
  transform_sample_counts(function(x) x / sum(x)) %>%
  subset_taxa(Module=="darkorange1")
sort(taxa_sums(invertOTUs))
## de1f0a9055c8bbecdda152cdfa8d8c83 d570f86ee0e1e76d190a80c9b1b5333f 
##                       0.05022154                       0.05956822 
## 19dd45f7f7adb9434de602dfe88fcceb 37ae5606defe4e69cb6de58d6f433034 
##                       0.13323196                       0.15154168 
## 132f7843a1ce20e8017e9d47b84eee2a 3d708fbb9bfd8f6e8cabfade38047a42 
##                       0.30303093                       0.30472115 
## 7833a2353efad879be2a75379abfba46 e7734250ef53cf043e4216b6e7590383 
##                       0.58715039                       0.83225779
kable(tax_table(invertOTUs))
Kingdom Phylum Class Order Family Genus Species OTUID Module
d570f86ee0e1e76d190a80c9b1b5333f Bacteria Firmicutes Clostridia Oscillospirales NA NA NA OTU_287 darkorange1
7833a2353efad879be2a75379abfba46 Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae NA NA OTU_320 darkorange1
e7734250ef53cf043e4216b6e7590383 Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae Clostridium_sensu_stricto_12 NA OTU_326 darkorange1
19dd45f7f7adb9434de602dfe88fcceb Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae Clostridium_sensu_stricto_10 NA OTU_334 darkorange1
37ae5606defe4e69cb6de58d6f433034 Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae NA NA OTU_336 darkorange1
132f7843a1ce20e8017e9d47b84eee2a Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae Lachnoclostridium NA OTU_358 darkorange1
3d708fbb9bfd8f6e8cabfade38047a42 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Xanthobacteraceae NA NA OTU_1054 darkorange1
de1f0a9055c8bbecdda152cdfa8d8c83 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Aquaspirillaceae Microvirgula NA OTU_1339 darkorange1 
#write.table(tax_table(invertOTUs), "invertpreyOTUs.csv", sep = ',')

Stacked bar plot of module membership

bp1 = psmelt(latitudeOTUs %>% filter_taxa(function(x) sum(x) >= 0.1, prune=T) ) %>% 
  arrange(latitude) %>% 
  mutate(across(Family, ~coalesce(.,paste(Order)))) %>% 
  mutate(across(Genus, ~coalesce(.,paste(Family)))) 

# arrange by a variable
site_order = rownames(arrange(data.frame(sample_data(latitudeOTUs)), latitude))
bp1 = bp1 %>% mutate(Sample = factor(Sample, levels=site_order))

bp1$site_id2 = factor(bp1$site_id, levels = unique(bp1$site_id))

plotcolors = c(RColorBrewer::brewer.pal(12,"Paired"), RColorBrewer::brewer.pal(8,"Spectral"))
ggplot(bp1) + geom_col(aes(Sample, Abundance, fill=Genus)) + 
  theme(axis.text.x = element_blank()) +
  scale_fill_manual(values = plotcolors) +
  facet_wrap(~site_id2, scales = "free")

Option B: Topological overlap method

Build network

# Turn adjacency into topological overlap
TOM = TOMsimilarity(as_adjacency_matrix(sc.se, sparse = F))
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
dissTOM = 1-TOM
# Call the hierarchical clustering function
geneTree = hclust(as.dist(dissTOM), method = "average")
# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 10
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, 
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize);
##  ..cutHeight not given, setting it to 0.984  ===>  99% of the (truncated) height range in dendro.
##  ..done.
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
# Plot the resulting clustering tree (dendrogram)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")

datExpr = data.frame(t(otu_table(ps6Top %>% 
                                   transform_sample_counts(function(x) x / sum(x)) %>%
                                   subset_taxa(taxa_names(ps6Top) %in% names(V(sc.se)))
                                   )))

# Get Module Eigengenes per cluster
MEs1 <- moduleEigengenes(datExpr, colors = dynamicColors)$eigengenes
# Reorder modules so similar modules are next to each other
MEs1 <- orderMEs(MEs1)
module_order = names(MEs1) %>% gsub("ME","", .)

# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs1)
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")

Plot module similarity

# Plot the result
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
# choose threshold: corresponds to correlation of 1-MEDissThres
MEDissThres = 0.4
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")

# Call an automatic merging function
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
##  mergeCloseModules: Merging modules whose distance is less than 0.4
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 7 module eigengenes in given set.
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 6 module eigengenes in given set.
##    Calculating new MEs...
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 6 module eigengenes in given set.
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs;

plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), 
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

Network with merged modules

V(sc.se)$color = mergedColors
plot(sc.se, vertex.label=NA,vertex.size=5)

Module-site relationships

MEs1 <- mergedMEs
# Reorder modules so similar modules are next to each other
MEs1 <- orderMEs(MEs1)
module_order = names(MEs1) %>% gsub("ME","", .)

MEs1$treatment = row.names(MEs1)
#MEs1$treatment = get_variable(ps6Top, "site_id")

# tidy & plot data
mME = MEs1 %>%
  pivot_longer(-treatment) %>%
  mutate(
    name = gsub("ME", "", name),
    name = factor(name, levels = module_order),
    treatment = factor(treatment, levels=site_order)
  )
mME %>% ggplot(., aes(x=treatment, y=name, fill=value)) +
  geom_tile() +
  theme_bw() +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(-1,1)) +
  theme(axis.text.x = element_text(angle=90)) +
  labs(title = "Module-site relationships sorted by bioclim_pca1", y = "Modules", fill="corr")

Relatedness within modules (SES MPD)

# add a "module" column to the tax table
ps.network = ps6Top %>% subset_taxa(taxa_names(ps6Top) %in% names(V(sc.se)))
tax_table(ps.network) = cbind(tax_table(ps.network), Module=mergedColors)

module_order = names(mergedMEs) %>% gsub("ME","", .)
# subset the tax table to the Module variable
mods = data.frame(tax_table(ps.network)[,"Module"])
# manually add module ID to taxa
mods = mods %>% mutate(mod1 = ifelse(Module==module_order[1],1,0),
                mod2 = ifelse(Module==module_order[2],1,0),
                mod3 = ifelse(Module==module_order[3],1,0),
                mod4 = ifelse(Module==module_order[4],1,0),
                mod5 = ifelse(Module==module_order[5],1,0),
                mod6 = ifelse(Module==module_order[6],1,0)) %>%
  select(-Module)
names(mods) = module_order
# set up data 
rooted = phy_tree(ps.network)
# otu = data.frame(otu_table(ps.network))
otu = mods
match <- match.phylo.comm(phy = rooted, comm = t(otu))

# SES MPD of modules
mod_SESMPD = ses.mpd(match$comm, cophenetic(match$phy), null.model = "taxa.labels")
kable(mod_SESMPD %>% select(c(ntaxa,mpd.obs.z,mpd.obs.p)))
ntaxa mpd.obs.z mpd.obs.p
brown 14 -2.4381088 0.011
turquoise 21 -1.1337982 0.122
blue 31 0.7085226 0.766
green 13 -1.4838197 0.079
black 10 0.0263959 0.512
red 11 0.2948844 0.611

Correlation plot

library(corrplot); library(psych)
dat = data.frame(sample_data(ps6Top)) %>%
  cbind(MEs1) %>%
  na.omit(.)

dat2 = data.frame(sample_data(ps6Top)) %>%
  select(-c(subspec,site_id,Description,site_id.1))

ct = corr.test(dat2, mergedMEs)

corrplot(ct$r, p.mat = ct$p.adj, sig.level = 0.1 , insig = "blank")

Regression of modules with environmental data

Latitude

ggplot(dat, aes(latitude, MEturquoise)) +geom_point()

mod1 = lm(MEturquoise~latitude, data=dat)
anova(mod1)
summary(mod1)
## 
## Call:
## lm(formula = MEturquoise ~ latitude, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.13670 -0.05399 -0.01955  0.02715  0.39253 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.341727   0.070826   4.825 5.47e-06 ***
## latitude    -0.007702   0.001583  -4.867 4.62e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09022 on 93 degrees of freedom
## Multiple R-squared:  0.203,  Adjusted R-squared:  0.1944 
## F-statistic: 23.69 on 1 and 93 DF,  p-value: 4.62e-06

Pitcher size and prey

ggplot(dat, aes(inverts, MEgreen)) +geom_point()

mod1 = lm(MEgreen~inverts, data=dat)
anova(mod1)
summary(mod1)
## 
## Call:
## lm(formula = MEgreen ~ inverts, data = dat)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.108191 -0.018466 -0.016542 -0.009679  0.214443 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0082204  0.0067655  -1.215    0.227
## inverts      0.0003765  0.0023748   0.159    0.874
## 
## Residual standard error: 0.04858 on 93 degrees of freedom
## Multiple R-squared:  0.0002702,  Adjusted R-squared:  -0.01048 
## F-statistic: 0.02514 on 1 and 93 DF,  p-value: 0.8744

OTUs in the modules

latitudeOTUs = ps.network %>% 
  transform_sample_counts(function(x) x / sum(x)) %>%
  subset_taxa(Module=="turquoise")
kable(tax_table(latitudeOTUs))
Kingdom Phylum Class Order Family Genus Species OTUID Module
92fb4682c111dd6e6f48af9df136e4cd Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter NA OTU_62 turquoise
38e25bbe6d992aa9f3a16960d46a6ee7 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter NA OTU_63 turquoise
90086558ccc079f6a635830df3a4ebb9 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter NA OTU_64 turquoise
b1f0f8eff3d0ea0d7a19fc6178e9ebeb Bacteria Bacteroidota Bacteroidia Flavobacteriales Weeksellaceae NA NA OTU_214 turquoise
168166cab31e0595c52ce7352f0db892 Bacteria Proteobacteria Alphaproteobacteria Acetobacterales Acetobacteraceae Acidisoma NA OTU_609 turquoise
df6b18712acf2c44bc154778f6b3ef09 Bacteria Acidobacteriota Acidobacteriae Acidobacteriales Acidobacteriaceae_(Subgroup_1) Granulicella NA OTU_650 turquoise
b57cb6018b602be504eca9417a0abe2c Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae Brevundimonas NA OTU_822 turquoise
572e3dfbe64c22ad15717c464f1534f0 Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae NA NA OTU_841 turquoise
31b69392eb47ad5bf4cc470deeb175da Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae NA NA OTU_989 turquoise
b71ab0ec9fae32ebaf2d8ec2ccc312b9 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Novosphingobium NA OTU_1001 turquoise
4b268cdb712389eb8bf03da6d48fe99f Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Rhodanobacteraceae NA NA OTU_1132 turquoise
aba772f0b9dfb11a3f3efac51e12148f Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Alkanindiges NA OTU_1139 turquoise
b868361da89d5b4a8941afce5c860246 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Chromobacteriaceae Aquitalea NA OTU_1284 turquoise
a93cbf059fd5f2636658ba59355e73d8 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Comamonadaceae NA NA OTU_1354 turquoise
2010445510c2a3c9db056f3e365441df Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Comamonadaceae NA NA OTU_1360 turquoise
0a17be34f84b48eecd0b82aec06a0929 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Comamonadaceae NA NA OTU_1398 turquoise
66d5abe6a7e075e65aec51076e7b5f05 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Burkholderiaceae Burkholderia-Caballeronia-Paraburkholderia NA OTU_1422 turquoise
882330144d2c06e8c47fc6a4f8988141 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1434 turquoise
72a610c3ddd4199e03399a0ae4e4c7a1 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1480 turquoise
504b1ad547b40a63f01095f8fbb593a4 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1481 turquoise
340375bcf95c24b77f70ba3af0dd6b6f Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA OTU_1497 turquoise 
#write.table(tax_table(yellowOTUs), "yellow-OTUs.csv", sep = ',')
invertOTUs = ps.network %>% 
  transform_sample_counts(function(x) x / sum(x)) %>%
  subset_taxa(Module=="green")
sort(taxa_sums(invertOTUs))
## de1f0a9055c8bbecdda152cdfa8d8c83 d570f86ee0e1e76d190a80c9b1b5333f 
##                       0.05022154                       0.05956822 
## 19dd45f7f7adb9434de602dfe88fcceb 37ae5606defe4e69cb6de58d6f433034 
##                       0.13323196                       0.15154168 
## 132f7843a1ce20e8017e9d47b84eee2a 3d708fbb9bfd8f6e8cabfade38047a42 
##                       0.30303093                       0.30472115 
## f09b6938fe1e93bb6c54874bc8d8c137 7a1b88dd4eff593cd8c2c13ab14da61e 
##                       0.35406299                       0.39515648 
## df3e11a2dea797be38778cc561a8a7f2 570d521e8740d2d782febcdd83fc19e8 
##                       0.42764880                       0.47467745 
## 7833a2353efad879be2a75379abfba46 e7734250ef53cf043e4216b6e7590383 
##                       0.58715039                       0.83225779 
## 2d7f2ebdb1abe1e1167c94d007f440d0 
##                       0.96257544
kable(tax_table(invertOTUs))
Kingdom Phylum Class Order Family Genus Species OTUID Module
d570f86ee0e1e76d190a80c9b1b5333f Bacteria Firmicutes Clostridia Oscillospirales NA NA NA OTU_287 green
7833a2353efad879be2a75379abfba46 Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae NA NA OTU_320 green
e7734250ef53cf043e4216b6e7590383 Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae Clostridium_sensu_stricto_12 NA OTU_326 green
19dd45f7f7adb9434de602dfe88fcceb Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae Clostridium_sensu_stricto_10 NA OTU_334 green
37ae5606defe4e69cb6de58d6f433034 Bacteria Firmicutes Clostridia Clostridiales Clostridiaceae NA NA OTU_336 green
132f7843a1ce20e8017e9d47b84eee2a Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae Lachnoclostridium NA OTU_358 green
df3e11a2dea797be38778cc561a8a7f2 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium NA OTU_913 green
7a1b88dd4eff593cd8c2c13ab14da61e Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA OTU_936 green
2d7f2ebdb1abe1e1167c94d007f440d0 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA OTU_963 green
3d708fbb9bfd8f6e8cabfade38047a42 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Xanthobacteraceae NA NA OTU_1054 green
570d521e8740d2d782febcdd83fc19e8 Bacteria Proteobacteria Gammaproteobacteria Enterobacterales NA NA NA OTU_1096 green
f09b6938fe1e93bb6c54874bc8d8c137 Bacteria Proteobacteria Gammaproteobacteria Enterobacterales Erwiniaceae NA NA OTU_1097 green
de1f0a9055c8bbecdda152cdfa8d8c83 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Aquaspirillaceae Microvirgula NA OTU_1339 green

Stacked bar plot of module membership

bp1 = psmelt(latitudeOTUs %>% filter_taxa(function(x) sum(x) >= 0.1, prune=T) ) %>% 
  arrange(latitude) %>% 
  mutate(across(Family, ~coalesce(.,paste(Order)))) %>% 
  mutate(across(Genus, ~coalesce(.,paste(Family)))) 

# arrange by a variable
site_order = rownames(arrange(data.frame(sample_data(latitudeOTUs)), latitude))
bp1 = bp1 %>% mutate(Sample = factor(Sample, levels=site_order))

bp1$site_id2 = factor(bp1$site_id, levels = unique(bp1$site_id))

plotcolors = c(RColorBrewer::brewer.pal(12,"Paired"), RColorBrewer::brewer.pal(8,"Spectral"))
ggplot(bp1) + geom_col(aes(Sample, Abundance, fill=Genus)) + 
  theme(axis.text.x = element_blank()) +
  scale_fill_manual(values = plotcolors) +
  facet_wrap(~site_id2, scales = "free")