Chapter 11 Dominance-microbiome analysis

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

11.0.1 Posterior estimates

The next step is to get the posterior estimates of beta, which is the parameter that links dominance with the microbiome. We employ a support threshold of ‘0.05’ to assign significance, meaning that parameters with posterior overlaps of <10% are considered significantly different.

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree %>%
        keep.tip(., tip=m$spNames) 

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    rename(intercept=2) %>%
    select(genome,dominance,groupvariable) %>%
    column_to_rownames(var="genome")

# Aggregate basal GIFT into elements
function_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap
postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.2, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

#Add functions heatmap
postestimates_tree <- gheatmap(postestimates_tree, function_table, offset=0.6, width=3.5, colnames=FALSE) +
            vexpand(.08) +
            coord_cartesian(clip = "off") +
            scale_fill_gradient(low = "#f4f4f4", high = "steelblue", na.value="white") +
            labs(fill="Function fullness")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add completeness barplots
postestimates_tree <- postestimates_tree +
            geom_fruit(data=genome_metadata,
            geom=geom_bar,
            grid.params=list(axis="x", text.size=2, nbreak = 1),
            axis.params=list(vline=TRUE),
            mapping = aes(x=length, y=genome, fill=completeness),
                 offset = 3.85,
                 orientation="y",
                 stat="identity") +
            scale_fill_gradient(low = "#cf8888", high = "#a2cc87") +
            labs(fill="Genome completeness")

postestimates_tree +
        vexpand(.25, 1) # expand top 

11.0.2 Positively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="dominance", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    select(genome,phylum,class,order,species,value) %>%
    tt()

genome	phylum	class	order	species	value
A_DRC05_bin.59	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910577345	1.000
C_HT_bin.102	p__Bacillota	c__Bacilli	o__Lactobacillales	s__Limosilactobacillus reuteri	0.998
A_HTC12_bin.48	p__Bacillota	c__Bacilli	o__Lactobacillales	s__Lactobacillus johnsonii	0.995
C_HT_bin.9	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__QXXE01 sp910589115	0.994
A_DTC03_bin.77	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp910586975	0.992
B_C04M3_bin.80	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.991
B_C05M2_bin.12	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp910578095	0.991
A_CDC10_bin.30	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910589485	0.990
B_C03M5_bin.33	p__Bacillota_B	c__Dehalobacteriia	o__UBA4068	s__	0.990
A_CDC13_bin.65	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp009774535	0.987
A_T3C12_bin.34	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-485 sp009775355	0.985
C_DR_bin.56	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-485 sp910578075	0.985
B_C04M3_bin.108	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter sp009936035	0.984
B_C06M3_bin.61	p__Bacillota	c__Bacilli	o__Erysipelotrichales	s__C-19 sp910576855	0.981
B_C12F1_bin.72	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910578515	0.980
C_HR_bin.153	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.980
A_T1C10_bin.3	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__UBA7173 sp002491305	0.979
B_C11F5_bin.68	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__UBA7173 sp001689485	0.979
B_C12F1_bin.61	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp009775535	0.977
A_DRC04_bin.49	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Enterenecus sp910587255	0.975
C_DR_bin.4	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Avispirillum sp011957885	0.974
A_HTC03_bin.32	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910587645	0.973
B_C04M5_bin.101	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Paramuribaculum sp910579675	0.973
B_C11F3_bin.29	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__CAG-95 sp009911035	0.973
A_CDC12_bin.65	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__Coproplasma sp003979335	0.972
B_C13F3_bin.77	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp910587915	0.972
A_CDC13_bin.13	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp002358555	0.971
B_C05M3_bin.38	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp910574915	0.971
B_C13F3_bin.53	p__Bacillota	c__Bacilli	o__Erysipelotrichales	s__Erysipelatoclostridium cocleatum	0.970
A_DTC13_bin.55	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__RACS-045 sp910579785	0.969
A_T2C12_bin.11	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__NSJ-51 sp910585415	0.969
B_C06M5_bin.62	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__CAG-475 sp910577815	0.968
B_C13F3_bin.5	p__Bacillota_A	c__Clostridia	o__TANB77	s__	0.967
B_C13F3_bin.97	p__Bacillota_A	c__Clostridia	o__TANB77	s__	0.967
C_CR_bin.18	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__VSOB01 sp910589075	0.966
C_DT_bin.159	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__	0.965
A_T2C05_bin.3	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Muribaculum gordoncarteri	0.964
B_C13F2_bin.43	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Enterenecus sp910587255	0.963
B_C04M3_bin.128	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Anaerotruncus sp000403395	0.962
B_C13F3_bin.35	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__MGBC136627 sp910585935	0.961
A_OPC13_bin.50	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__UBA1405 sp910589145	0.960
B_C05M1_bin.75	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Clostridium_Q sp009911305	0.959
B_C13F3_bin.70	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.959
A_CDC05_bin.89	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__	0.958
A_OPC10_bin.79	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910574345	0.957
C_OP_bin.44	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp910586205	0.957
C_HT_bin.137	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.956
C_CR_bin.157	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__	0.954
C_T3_bin.64	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__NSJ-51 sp910585415	0.954
C_DR_bin.31	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Anaerotignum sp910576545	0.953
A_OPC10_bin.40	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.952
B_C12F2_bin.47	p__Bacillota	c__Bacilli	o__Lactobacillales	s__Enterococcus_D casseliflavus	0.952
C_HR_bin.13	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Evtepia sp910584805	0.952
C_T2_bin.55	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.952
B_C04M2_bin.33	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Enterenecus sp910585265	0.951
B_C04M3_bin.31	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Angelakisella sp013316495	0.951
A_OPC03_bin.64	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Acetatifactor sp910589655	0.950
A_T1C11_bin.9	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-485 sp002362485	0.949
B_C04M3_bin.18	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp910586445	0.949
B_C03M3_bin.127	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter sp009936035	0.948
B_C04M4_bin.45	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.948
A_T3C03_bin.17	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__QANA01 sp910588425	0.947
B_C04M2_bin.76	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__CAG-56 sp004793585	0.945
A_HRC11_bin.35	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.944
B_C04M4_bin.5	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Choladocola sp009774145	0.944
C_DT_bin.37	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter sp009936035	0.944
C_DR_bin.21	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__MD308 sp910578455	0.943
A_CDC03_bin.78	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Lawsonibacter sp910577805	0.942
B_C13F5_bin.28	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910576785	0.942
A_HRC10_bin.28	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__Gallimonas sp910578065	0.941
A_T3C11_bin.23	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__RGIG4284 sp910576205	0.941
A_HTC12_bin.21	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910575585	0.940
A_T2C05_bin.28	p__Bacillota	c__Bacilli	o__Erysipelotrichales	s__MGBC163490 sp910588075	0.940
B_C10F2_bin.5	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.940
C_HR_bin.146	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Choladocola sp910575445	0.940
B_C05M5_bin.57	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910587785	0.939
B_C10F5_bin.45	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Dysosmobacter sp000403435	0.939
B_C12F3_bin.66	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__	0.939
B_C11F5_bin.87	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.938
A_CRC05_bin.32	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp910579075	0.937
B_C04M1_bin.34	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910579965	0.937
B_C04M3_bin.90	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__14-2 sp910580015	0.937
C_T2_bin.30	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.934
C_T3_bin.41	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Acetatifactor sp910584235	0.934
C_CD_bin.124	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910579965	0.933
C_OP_bin.73	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910585435	0.933
C_HR_bin.50	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910588325	0.932
A_CDC03_bin.64	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.930
B_C03M3_bin.121	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Lachnoclostridium_B sp910574185	0.930
B_C10F1_bin.59	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Choladocola sp910583895	0.929
B_C10F5_bin.36	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3282 sp003611805	0.928
B_C11F5_bin.21	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__COE1 sp000403335	0.928
A_CDC10_bin.13	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC113161 sp910587565	0.927
A_DRC13_bin.41	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter muris	0.927
C_CD_bin.62	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__1XD42-69 sp009911505	0.925
C_OP_bin.15	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Roseburia sp910584215	0.925
A_CDC05_bin.91	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.923
C_T2_bin.92	p__Bacillota_A	c__Clostridia	o__TANB77	s__	0.923
A_CDC05_bin.30	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.922
A_CRC04_bin.71	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Sporofaciens sp910574885	0.922
B_C13F3_bin.99	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Merdisoma sp910574255	0.922
B_C04M1_bin.8	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Merdisoma sp910576325	0.921
B_C03M2_bin.44	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella muris	0.920
B_C06M5_bin.16	p__Actinomycetota	c__Coriobacteriia	o__Coriobacteriales	s__Adlercreutzia muris	0.920
C_T3_bin.99	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910577915	0.920
B_C05M1_bin.6	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp910584845	0.919
A_DRC06_bin.37	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp009775265	0.918
B_C12F1_bin.38	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Dysosmobacter sp910588005	0.918
B_C10F1_bin.10	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp910577765	0.917
B_C13F3_bin.93	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Anaerotruncus sp003612625	0.917
A_HRC05_bin.71	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC100320 sp910577995	0.915
A_T3C13_bin.3	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.914
B_C03M1_bin.60	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC100798 sp910584975	0.911
C_OP_bin.100	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.910
C_DT_bin.155	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.909
A_T1C03_bin.24	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UMGS1370 sp910577645	0.907
B_C12F1_bin.23	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910579965	0.907
A_T2C13_bin.14	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__UBA7173 sp001689685	0.906
A_CRC05_bin.73	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC100320 sp910588305	0.905
B_C03M1_bin.29	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp009911065	0.905
B_C04M5_bin.14	p__Actinomycetota	c__Coriobacteriia	o__Coriobacteriales	s__Adlercreutzia caecimuris	0.905
C_HR_bin.28	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Sporofaciens sp910585725	0.903
B_C04M5_bin.93	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__Gallimonas sp910588465	0.902

11.0.3 Negatively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="dominance", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    select(genome,phylum,class,order,species,value) %>%
    tt()

genome	phylum	class	order	species	value
A_T3C05_bin.41	p__Pseudomonadota	c__Gammaproteobacteria	o__Burkholderiales	s__Parasutterella excrementihominis	0.095

11.0.4 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="dominance") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="dominance") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% estimate$genome) %>%
    arrange(match(genome, estimate$genome)) %>%
     select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    select(colors) %>%
    pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

11.1 Predict responses to dominance

Based on the fitted model, we can predict the structure of the microbial community for each level of dominance for a certain value of sequencing depth.

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = seq(0,1,0.1)
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "dominance", 
                      non.focalVariables = list(logseqdepth=list(1)), 
                      ngrid=gradientlength) %>%
             predict(m, Gradient = ., expected = TRUE)

# weights:  3 (2 variable)
initial  value 162.196440 
final  value 133.205769 
converged

predY <- pred %>%
        do.call(rbind,.) %>%
        as.data.frame() %>%
        mutate(dominance=rep(gradient,1000)) %>%
        pivot_longer(-c(dominance), names_to = "genome", values_to = "value")

11.1.1 Plot responses to dominance

The estimated responses of genomes exhibiting significant positive and negative association with dominance.

# Select desired support threshold
support=0.9
negsupport=1-support

#Get phylum colors from the EHI standard
phylum_colors <- genome_metadata %>%
    left_join(read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv"), by=join_by(phylum == phylum)) %>%
    arrange(match(genome, genome_tree$tip.label)) %>%
    select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    #slice(2:5) %>%
    select(colors) %>%
    pull()

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="dominance") %>%
    select(genome,trend) %>%
    left_join(predY, by=join_by(genome==genome)) %>%
    #filter(trend != "Neutral") %>%
    #filter(genome %in% predY_asc)  %>% #only display genomes with contrasting dynamics across treatments
    group_by(genome, trend, dominance) %>%
    summarize(value = mean(value, na.rm = TRUE)) %>%
    left_join(genome_metadata, by=join_by(genome == genome)) %>%
    ggplot(aes(x=dominance, y=value, group=genome, color=phylum, linetype=trend)) + 
        geom_line() +
        scale_linetype_manual(values=c("solid","dashed","solid")) +
        scale_color_manual(values=phylum_colors) +
        facet_grid(fct_rev(trend) ~ phylum) +
        labs(y="Genome abundance (log)",x="Dominance") +
        theme(legend.position = "none") +
        theme_minimal() +
         theme(legend.position = "none",
               axis.text.x = element_text(angle = 45, hjust = 0.8,),
               axis.line.x = element_line(size = 0.5, linetype = "solid", colour = "black"),
)

11.2 Plot traits linked to trends

11.2.1 Test functional trait differences

Using Wilcoxon test with Benjamini-Hochberg procedure for FDR, identify functional traits that differ between genomes that are either positively or negatively association with dominance.

postestimate_test <- getPostEstimate(hM=m, parName="Beta")$support %>%
      as.data.frame() %>%
      mutate(variable=m$covNames) %>%
      pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
      mutate(trend = case_when(
            value >= support ~ "Positive",
            value <= negsupport ~ "Negative",
            TRUE ~ "Neutral")) %>%
      filter(variable=="dominance", trend!="Neutral") %>%
      select(genome,trend) %>%
      left_join(function_table %>% rownames_to_column(var="genome"), by=join_by(genome==genome)) %>%
      pivot_longer(-c(genome,trend), names_to = "trait", values_to = "value") %>%
      group_by(trait) %>%
      summarise(p_value = wilcox.test(value ~ trend)$p.value) %>%
      mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
      filter(p_adjust < 0.05)

11.2.2 Community functional predictions

11.2.2.1 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- predY %>%
  group_by(dominance, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    select(-row_id) %>%
    column_to_rownames(var = "dominance") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="dominance")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "dominance") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
B0103	4.104165e-04	-3.291650e-05	7.890253e-04	0.833	0.167
D0305	1.666043e-03	-5.160041e-04	4.749368e-03	0.818	0.182
B0708	2.496614e-03	-1.254454e-03	6.849490e-03	0.809	0.191
D0511	3.346787e-03	-1.331733e-03	9.059830e-03	0.807	0.193
B0706	3.490214e-03	-1.338794e-03	9.164635e-03	0.793	0.207
D0602	2.314528e-03	-1.953434e-03	6.825266e-03	0.777	0.223
B0102	7.873893e-04	-1.022802e-03	3.281036e-03	0.757	0.243
B0705	9.483230e-04	-9.210338e-04	2.906552e-03	0.751	0.249
D0512	2.232714e-03	-2.285578e-03	7.141727e-03	0.741	0.259
B0214	2.183718e-03	-2.327279e-03	7.105594e-03	0.738	0.262
B0803	5.248397e-04	-3.941343e-04	1.356366e-03	0.735	0.265
D0513	1.780906e-03	-1.771384e-03	5.921903e-03	0.733	0.267
D0203	1.778707e-03	-1.126479e-03	5.615073e-03	0.729	0.271
D0910	2.284603e-03	-2.433390e-03	7.635939e-03	0.726	0.274
B0805	9.467965e-04	-6.959032e-04	3.405693e-03	0.702	0.298
B0106	8.383317e-04	-1.385429e-03	2.879423e-03	0.695	0.305
B0802	1.100765e-03	-1.174843e-03	3.581643e-03	0.691	0.309
D0613	2.110746e-03	-4.926596e-03	8.185827e-03	0.691	0.309
D0706	6.825220e-04	-1.505567e-03	3.069765e-03	0.691	0.309
B0707	6.780900e-04	-9.208398e-04	2.444123e-03	0.678	0.322
B0104	9.731324e-04	-2.550485e-03	4.280839e-03	0.666	0.334
D0501	7.046608e-04	-2.116070e-03	3.379259e-03	0.651	0.349
B0711	2.470479e-04	-2.860633e-03	3.009486e-03	0.642	0.358
D0209	8.235682e-04	-2.423924e-03	4.158671e-03	0.640	0.360
B0703	4.234373e-04	-9.192862e-04	1.949271e-03	0.638	0.362
B0702	9.399523e-04	-2.801397e-03	4.548685e-03	0.637	0.363
B0704	1.053597e-03	-3.933395e-03	5.536359e-03	0.633	0.367
S0105	4.898645e-04	-5.154987e-03	4.877545e-03	0.631	0.369
D0309	8.851932e-04	-1.838238e-03	4.166754e-03	0.627	0.373
D0210	8.947192e-04	-2.485486e-03	4.495952e-03	0.626	0.374
B0101	2.938283e-04	-1.436322e-03	2.071892e-03	0.619	0.381
D0509	1.587614e-04	-1.685350e-03	1.914757e-03	0.616	0.384
D0308	1.089061e-03	-4.002174e-03	6.379185e-03	0.609	0.391
D0505	1.644630e-04	-4.384325e-04	8.456321e-04	0.606	0.394
B0216	2.189220e-04	-2.577254e-03	2.838512e-03	0.603	0.397
B0601	6.568203e-05	-4.143817e-03	3.400514e-03	0.600	0.400
B0712	1.047982e-04	-6.265841e-04	7.828061e-04	0.600	0.400
D0702	6.593381e-04	-2.406870e-03	4.153419e-03	0.592	0.408
D0205	6.862388e-04	-2.074949e-03	3.811739e-03	0.590	0.410
B0701	5.829055e-04	-4.109473e-03	4.706115e-03	0.585	0.415
B0105	2.354402e-04	-2.015111e-03	2.103642e-03	0.582	0.418
D0708	1.264386e-04	-9.750279e-04	1.247018e-03	0.577	0.423
D0302	4.654686e-04	-1.831433e-03	2.975143e-03	0.566	0.434
D0212	3.251211e-04	-2.898410e-03	3.765002e-03	0.557	0.443
B0207	2.698546e-04	-2.666829e-03	3.620693e-03	0.553	0.447
D0801	3.365227e-06	-4.240681e-05	5.356688e-05	0.542	0.458
D0802	3.365227e-06	-4.240681e-05	5.356688e-05	0.542	0.458
D0306	4.203091e-04	-3.016939e-03	4.039724e-03	0.530	0.470
D0208	6.868435e-05	-5.145522e-04	6.837883e-04	0.527	0.473
D0207	1.451571e-03	-2.996464e-03	7.084396e-03	0.522	0.478
B0604	1.273429e-03	-3.898953e-03	7.676924e-03	0.511	0.489
B0401	8.192590e-04	-1.684257e-03	4.410619e-03	0.507	0.493
D0705	4.009476e-04	-1.299421e-03	2.675762e-03	0.503	0.497
D0310	6.911104e-05	-4.605457e-03	5.099046e-03	0.499	0.501
D0516	1.170211e-05	-1.431965e-04	1.830712e-04	0.475	0.525
D0206	5.240233e-04	-4.955146e-03	6.503960e-03	0.472	0.528
B0403	2.591203e-04	-1.257287e-03	2.162790e-03	0.462	0.538
D0213	3.623353e-04	-2.930587e-03	4.142586e-03	0.450	0.550
B0402	2.095399e-04	-2.148903e-03	3.402101e-03	0.443	0.557

# Negatively associated
element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D0607	-8.291600e-04	-2.269564e-03	-4.498889e-05	0.023	0.977
D0503	-1.421521e-03	-4.523134e-03	-6.849123e-06	0.026	0.974
D0601	-2.841810e-03	-9.009419e-03	-9.672588e-06	0.051	0.949
B0218	-2.763354e-03	-7.009181e-03	-1.883848e-06	0.084	0.916
D0104	-9.962497e-04	-2.089788e-03	-1.657664e-07	0.089	0.911
D0907	-4.481490e-03	-9.939284e-03	1.482275e-05	0.101	0.899
B0221	-3.925856e-03	-9.035304e-03	1.408680e-04	0.110	0.890
B0307	-2.345379e-03	-5.271998e-03	1.514725e-04	0.112	0.888
D0902	-2.847800e-03	-9.003818e-03	2.066314e-05	0.117	0.883
D0103	-4.144738e-04	-9.605540e-04	2.374533e-06	0.121	0.879
B0902	-3.958902e-07	-9.677130e-07	7.274520e-10	0.129	0.871
B0219	-4.930725e-05	-1.223234e-04	1.034796e-06	0.154	0.846
D0903	-4.975555e-05	-1.227435e-04	1.031276e-06	0.154	0.846
D0805	-4.844271e-07	-1.267004e-06	8.087070e-09	0.155	0.845
D0806	-1.598609e-07	-4.181113e-07	2.668733e-09	0.155	0.845
D0904	-4.844271e-07	-1.267004e-06	8.087070e-09	0.155	0.845
D0603	-1.674390e-05	-4.082512e-05	2.538289e-07	0.156	0.844
D0611	-4.927112e-05	-1.218486e-04	1.072277e-06	0.156	0.844
D0811	-1.666549e-04	-5.410216e-04	5.817357e-07	0.156	0.844
B0801	-3.005124e-05	-7.223503e-05	3.493701e-07	0.157	0.843
B0212	-3.099222e-03	-8.050103e-03	1.211155e-03	0.171	0.829
D0517	-2.494875e-05	-6.237695e-05	9.129985e-07	0.179	0.821
D0101	-1.719170e-03	-4.309367e-03	4.446062e-04	0.180	0.820
B0903	-1.238816e-04	-2.904086e-04	3.612286e-05	0.191	0.809
B0217	-2.355016e-03	-6.195672e-03	5.156159e-04	0.192	0.808
D0508	-1.453102e-04	-3.489720e-04	4.932707e-05	0.200	0.800
B0303	-5.689171e-04	-1.693970e-03	2.476785e-04	0.202	0.798
B0309	-2.406975e-05	-5.922252e-05	4.563999e-06	0.207	0.793
B0209	-2.300098e-03	-6.685690e-03	1.699779e-03	0.214	0.786
B0210	-3.035285e-03	-8.721838e-03	2.095731e-03	0.223	0.777
B0208	-2.249476e-03	-6.605608e-03	1.757571e-03	0.231	0.769
B0213	-2.467121e-03	-7.453025e-03	1.806678e-03	0.231	0.769
B0206	-1.728677e-03	-5.374166e-03	1.637505e-03	0.236	0.764
B1028	-4.016906e-04	-1.330101e-03	1.395929e-04	0.239	0.761
D0807	-2.676669e-04	-7.656125e-04	1.502580e-04	0.241	0.759
D0301	-4.778863e-05	-1.553416e-04	5.215710e-05	0.243	0.757
D0817	-2.891783e-04	-7.957855e-04	1.956856e-04	0.247	0.753
B0205	-1.688480e-03	-5.329931e-03	1.387674e-03	0.248	0.752
B0215	-2.461911e-03	-7.783948e-03	2.135384e-03	0.253	0.747
D0816	-6.115945e-04	-1.728626e-03	5.384457e-04	0.253	0.747
D0606	-1.856943e-04	-5.267146e-04	1.291961e-04	0.254	0.746
D0908	-1.741660e-03	-5.548600e-03	2.346721e-03	0.261	0.739
D0704	-2.432267e-03	-7.507592e-03	2.002404e-03	0.262	0.738
D0510	-5.647592e-05	-1.933528e-04	6.450299e-05	0.264	0.736
D0819	-1.086629e-08	-3.131216e-08	4.582416e-09	0.264	0.736
B0302	-5.661645e-05	-1.488243e-04	1.825840e-05	0.268	0.732
B0602	-7.347080e-04	-2.809640e-03	1.014026e-03	0.277	0.723
D0502	-2.001295e-03	-7.106487e-03	2.716307e-03	0.279	0.721
B1012	-2.555086e-04	-9.382364e-04	3.963025e-04	0.281	0.719
B0204	-8.212644e-04	-2.986972e-03	1.055787e-03	0.288	0.712
D0518	-1.869769e-04	-5.865696e-04	1.782459e-04	0.292	0.708
D0912	-5.608830e-04	-1.599417e-03	2.192870e-04	0.296	0.704
B0211	-8.700513e-04	-3.370054e-03	1.183226e-03	0.297	0.703
B0901	-2.927109e-05	-1.030717e-04	4.511132e-05	0.297	0.703
B0220	-3.421852e-04	-1.148521e-03	4.808089e-04	0.298	0.702
B0709	-8.196714e-05	-3.169548e-04	1.688408e-04	0.311	0.689
S0104	-6.336098e-04	-3.145604e-03	2.086668e-03	0.315	0.685
D0307	-7.036201e-04	-2.693563e-03	1.193173e-03	0.316	0.684
S0201	-1.252561e-03	-5.825400e-03	3.750881e-03	0.317	0.683
D0610	-2.612261e-04	-9.770819e-04	4.358967e-04	0.322	0.678
D0701	-4.259483e-05	-1.519154e-04	5.245908e-05	0.335	0.665
D0507	-1.433011e-04	-6.875030e-04	3.570201e-04	0.353	0.647
D0201	-1.644238e-03	-5.941041e-03	9.391648e-04	0.354	0.646
S0103	-4.589567e-05	-1.298387e-03	1.365016e-03	0.360	0.640
S0301	-1.031175e-03	-5.708423e-03	4.142674e-03	0.367	0.633
D0604	-1.062937e-03	-7.024071e-03	4.896637e-03	0.368	0.632
B0804	-6.188490e-04	-3.370393e-03	1.509535e-03	0.373	0.627
D0211	-2.689049e-04	-1.646414e-03	1.112083e-03	0.377	0.623
D0202	-6.174943e-04	-2.561183e-03	1.020619e-03	0.380	0.620
B1004	-5.911958e-05	-4.524567e-04	3.073879e-04	0.386	0.614
D0608	-2.152204e-05	-1.205862e-04	7.870969e-05	0.388	0.612
S0202	-9.007486e-05	-4.947459e-04	3.222497e-04	0.402	0.598
D0906	-1.036392e-06	-4.975518e-06	2.627358e-06	0.406	0.594
D0504	-7.944653e-05	-4.601832e-04	2.629917e-04	0.409	0.591
B0605	-8.736315e-04	-5.535142e-03	3.899984e-03	0.426	0.574
D0506	-1.005735e-04	-1.237638e-03	9.332254e-04	0.427	0.573
D0102	-2.421593e-04	-1.847734e-03	1.172253e-03	0.446	0.554
D0609	-6.601686e-05	-1.798551e-03	1.907520e-03	0.453	0.547
B0603	-1.413073e-04	-2.275988e-03	1.748526e-03	0.464	0.536
D0905	-8.787333e-05	-4.054964e-03	4.302750e-03	0.466	0.534
D0911	-8.869413e-05	-4.089514e-03	4.404582e-03	0.469	0.531
D0304	-1.179077e-04	-2.406748e-03	2.089449e-03	0.473	0.527
S0101	-7.349656e-05	-1.347788e-03	9.602682e-04	0.489	0.511
D0901	-5.324917e-04	-7.268231e-03	5.322711e-03	0.521	0.479
D0204	-2.675759e-04	-2.863755e-03	2.203051e-03	0.524	0.476
B0710	-4.912803e-04	-2.893019e-03	1.274017e-03	0.571	0.429

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  select(-negative_support) %>%
  rename(support=positive_support) %>%
  slice(1:5)

# Negatively associated
negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  select(-positive_support) %>%
  rename(support=negative_support) %>%
  mutate(support=-support) %>%
  slice(1:10) %>%
  arrange(-mean)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),negative$trait))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

11.2.2.2 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- predY %>%
  group_by(dominance, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    select(-row_id) %>%
    column_to_rownames(var = "dominance") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="dominance")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "dominance") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
B08	0.0006145788	-0.0002240063	0.001682819	0.816	0.184
B07	0.0011878162	-0.0006215992	0.003301166	0.760	0.240
D05	0.0004984042	-0.0004281026	0.001607224	0.745	0.255
B01	0.0005897564	-0.0010033494	0.002038586	0.673	0.327
S01	0.0001177181	-0.0024796017	0.002041800	0.599	0.401
D02	0.0003166735	-0.0015917805	0.002259158	0.563	0.437
D03	0.0004249427	-0.0012783718	0.002415569	0.536	0.464
B04	0.0003926786	-0.0016173340	0.003010112	0.438	0.562

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
B03	-5.085998e-04	-0.0011041092	-7.200935e-09	0.096	0.904
D01	-1.018245e-03	-0.0025075048	1.925184e-05	0.106	0.894
B10	-5.323854e-05	-0.0001321572	1.014008e-05	0.149	0.851
B02	-1.545790e-03	-0.0043634642	7.360975e-04	0.198	0.802
D08	-7.589549e-05	-0.0001818120	1.929862e-05	0.211	0.789
D09	-7.370514e-04	-0.0022990535	5.608475e-04	0.222	0.778
B09	-3.453682e-05	-0.0000961876	2.775886e-05	0.245	0.755
D07	-6.584616e-04	-0.0025220716	1.046351e-03	0.302	0.698
S02	-5.550692e-04	-0.0025357452	1.590999e-03	0.320	0.680
S03	-1.031175e-03	-0.0057084232	4.142674e-03	0.367	0.633
B06	-3.095813e-04	-0.0021675065	1.250320e-03	0.405	0.595
D06	-3.055253e-04	-0.0015702819	8.006510e-04	0.466	0.534

#lm option (incoherent with Bayesian framework)
#function_predictions <- community_functions %>%
#   bind_rows() %>%
#   pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
#  group_by(trait) %>%
#  mutate(dominance=as.numeric(dominance)) %>%
#  summarize(model=list(lm(value ~ dominance))) %>%
#  ungroup() %>%
#  select(trait,model) %>%
#  mutate(estimate = map_dbl(model, ~broom::tidy(.) %>% filter(term == "dominance") %>% pull(estimate))) %>%
#  mutate(p_value = map_dbl(model, ~broom::tidy(.) %>% filter(term == "dominance") %>% pull(p.value))) %>%
#  mutate(p_value_adj = p.adjust(p_value, method = "bonferroni")) %>%
#  select(-model) %>%
#  arrange(-estimate) %>%
#  filter(p_value_adj < 0.05) 

#brms option (very lengthy) > even more if group is added as random effect (not in this example)
#function_predictions <- community_functions %>%
#   bind_rows() %>%
#   pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
#  mutate(dominance=as.numeric(dominance)) %>%
#  group_by(trait) %>%
#  do(model=list(brm(value ~ dominance, data=.))) %>%
#  ungroup() %>%
#  select(trait,model) %>%
#  mutate(estimate = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(estimate))) %>%
#  mutate(error = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(std.error))) %>%
#  mutate(error = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(conf.low))) %>%
#  mutate(error = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(conf.high))) %>%
#  arrange(-estimate)

community_functions %>%
    bind_rows() %>%
    pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
  filter(trait %in% function_predictions$trait) %>%
  mutate(trait=factor(trait, levels=function_predictions$trait)) %>%
    mutate(dominance=as.numeric(dominance)) %>%
    ggplot(aes(x=dominance, y=value)) +
          geom_smooth(method = lm, formula = y ~ x, se = TRUE) +
          #geom_smooth(method = lm, formula = y ~ splines::bs(x, 3), se = TRUE) +
          facet_wrap(~trait, ncol=5, scales="free") +
          theme_minimal() +
          labs(x="Dominance",y="Metabolic Capacity Index")