Chapter 9 Functional differences
load("data/data.Rdata")
sample_metadata$fox_behaviour <- factor(sample_metadata$fox_behaviour, levels=c("tame", "aggr", "unsel"))
sample_metadata$gut_location <- factor(sample_metadata$gut_location, levels=c("F_colon", "M_colon", "D_ileum", "K_ileum"))
behaviour_colors <- c(
tame = "#1f77b4",
aggr = "#d62728",
unsel = "#2ca02c")
gut_location_colors<- c(
F_colon = "#d6604d",
M_colon = "#542788",
D_ileum = "#fdb863",
K_ileum = "#bb99d8")
genome_counts_filt <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),]
rownames(genome_counts_filt) <- NULL#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts[, !grepl("^S", colnames(genome_gifts))],GIFT_db) # Remove structure traits
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
#Get community-weighed average GIFTs per sample
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
#genome_counts_row <- rownames_to_column(genome_counts_row, "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)9.0.1 Metabolic capacity index (MCI)
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = "sample") %>%
group_by(fox_behaviour) %>%
summarise(MCI = mean(value), sd = sd(value)) %>%
tt()| fox_behaviour | MCI | sd |
|---|---|---|
| tame | 0.3810984 | 0.02235157 |
| aggr | 0.3850806 | 0.02104977 |
| unsel | 0.3739909 | 0.02607507 |
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = "sample") %>%
group_by(gut_location) %>%
summarise(MCI = mean(value), sd = sd(value)) %>%
tt()| gut_location | MCI | sd |
|---|---|---|
| F_colon | 0.3642809 | 0.012792940 |
| M_colon | 0.3864882 | 0.014667724 |
| D_ileum | 0.3701253 | 0.025372898 |
| K_ileum | 0.4053539 | 0.001127473 |
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = "sample")
shapiro.test(MCI$value) %>%
tidy()# A tibble: 1 × 3
statistic p.value method
<dbl> <dbl> <chr>
1 0.908 0.00150 Shapiro-Wilk normality test# A tibble: 1 × 4
statistic p.value parameter method
<dbl> <dbl> <int> <chr>
1 1.35 0.510 2 Kruskal-Wallis rank sum test# A tibble: 1 × 4
statistic p.value parameter method
<dbl> <dbl> <int> <chr>
1 23.6 0.0000296 3 Kruskal-Wallis rank sum testWarning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute exact p-value with ties
Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute exact p-value with ties
Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute exact p-value with ties
Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute exact p-value with ties
Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute exact p-value with ties
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: MCI$value and MCI$gut_location
F_colon M_colon D_ileum
M_colon 0.00462 - -
D_ileum 0.44914 0.18855 -
K_ileum 0.00047 0.00478 0.00232
P value adjustment method: fdr The Shapiro test indicates the data is not normally distributed (p = 0.0015), so we use Kruskal-Wallis (non-parametric ANOVA) to test differences between the groups.
Fox behaviour has no signficant effect (p-value = 0.51).
Gut location has significant effect (p-value = 2.956672e-05), so we will test the pairwise comparisons.
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = GIFTs_functions_community_dist ~ gut_location, data = sample_metadata_filtered, permutations = 999, by = "terms")
Df SumOfSqs R2 F Pr(>F)
gut_location 3 1.6619 0.50017 14.009 0.001 ***
Residual 42 1.6607 0.49983
Total 45 3.3226 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MCI %>%
ggplot(aes(x=gut_location, y=value, color = gut_location, fill=gut_location)) +
geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
geom_jitter(width = 0.1, show.legend = TRUE) +
scale_color_manual(values=gut_location_colors)+
scale_fill_manual(values=gut_location_colors)+
theme_minimal() +
theme(
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(var="sample")%>%
pivot_longer(!sample,names_to="trait",values_to="gift") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
mutate(functionid = substr(trait, 1, 3)) %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_element ~ GIFT_db$Element[match(trait, GIFT_db$Code_element)],
TRUE ~ trait
)) %>%
mutate(functionid = case_when(
functionid %in% GIFT_db$Code_function ~ GIFT_db$Function[match(functionid, GIFT_db$Code_function)],
TRUE ~ functionid
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Element))) %>%
mutate(functionid=factor(functionid,levels=unique(GIFT_db$Function))) %>%
ggplot(aes(x=sample,y=trait,fill=gift)) +
geom_tile(colour="white", linewidth=0.2)+
scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
facet_grid(functionid ~ fox_behaviour, scales="free",space="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.text.y = element_text(size=8),
strip.text.y = element_text(angle = 0)
) +
labs(y="Traits",x="Samples",fill="GIFT")
GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(var="sample")%>%
pivot_longer(!sample,names_to="trait",values_to="gift") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
mutate(functionid = substr(trait, 1, 3)) %>%
mutate(trait = case_when(
trait %in% GIFT_db$Code_element ~ GIFT_db$Element[match(trait, GIFT_db$Code_element)],
TRUE ~ trait
)) %>%
mutate(functionid = case_when(
functionid %in% GIFT_db$Code_function ~ GIFT_db$Function[match(functionid, GIFT_db$Code_function)],
TRUE ~ functionid
)) %>%
mutate(trait=factor(trait,levels=unique(GIFT_db$Element))) %>%
mutate(functionid=factor(functionid,levels=unique(GIFT_db$Function))) %>%
ggplot(aes(x=sample,y=trait,fill=gift)) +
geom_tile(colour="white", linewidth=0.2)+
scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
facet_grid(functionid ~ gut_location, scales="free",space="free") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.text.y = element_text(size=10),
strip.text.y = element_text(angle = 0)
) +
labs(y="Traits",x="Samples",fill="GIFT")
9.0.2 Wilcoxon
9.0.2.1 Community elements differences
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(sample_metadata %>% select(sample,gut_location,fox_behaviour), by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,gut_location,fox_behaviour),
names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = kruskal.test(value ~ gut_location)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
rownames_to_column("Elements") %>%
left_join(uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(sample_metadata %>% select(sample,gut_location,fox_behaviour), by = join_by(sample == sample))
element_gift_filt %>%
select(-c(sample,fox_behaviour))%>%
group_by(gut_location) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element)) Elements F_colon M_colon D_ileum K_ileum Function
1 B0101 9.042021e-01 0.9452185000 8.483872e-01 0.9922576000 Nucleic acid biosynthesis_Inosinic acid (IMP)
2 B0102 8.980429e-01 0.9521216000 8.814970e-01 0.9943065000 Nucleic acid biosynthesis_Uridylic acid (UMP)
3 B0103 9.951756e-01 1.0000000000 9.995860e-01 1.0000000000 Nucleic acid biosynthesis_UDP/UTP
4 B0104 4.694174e-01 0.7673844000 5.957346e-01 0.9676218000 Nucleic acid biosynthesis_CDP/CTP
5 B0105 7.311438e-01 0.8748024000 7.400919e-01 0.9838109000 Nucleic acid biosynthesis_ADP/ATP
6 B0106 7.701193e-01 0.9008721000 8.581334e-01 0.9919054000 Nucleic acid biosynthesis_GDP/GTP
7 B0204 6.525033e-01 0.8419967000 8.029068e-01 0.9996184000 Amino acid biosynthesis_Serine
8 B0205 6.319491e-01 0.5940950000 6.627389e-01 0.6003995000 Amino acid biosynthesis_Threonine
9 B0206 8.270660e-01 0.9140888000 8.975307e-01 0.9990401000 Amino acid biosynthesis_Cysteine
10 B0207 6.737819e-01 0.7691794000 7.416252e-01 0.8596646000 Amino acid biosynthesis_Methionine
11 B0208 5.442286e-01 0.8057401000 7.529939e-01 0.9991327000 Amino acid biosynthesis_Valine
12 B0209 6.149349e-01 0.7086889000 7.212632e-01 0.7999369000 Amino acid biosynthesis_Isoleucine
13 B0210 5.010157e-01 0.7767003000 5.476730e-01 0.9975052000 Amino acid biosynthesis_Leucine
14 B0211 8.093479e-01 0.9126115000 9.190120e-01 0.9996184000 Amino acid biosynthesis_Lysine
15 B0212 6.776948e-01 0.6714719000 7.344924e-01 0.7053871000 Amino acid biosynthesis_Arginine
16 B0213 7.602218e-01 0.8909218000 7.917562e-01 0.9826544000 Amino acid biosynthesis_Proline
17 B0214 1.675417e-01 0.6677091000 3.346663e-01 0.9980025000 Amino acid biosynthesis_Glutamate
18 B0215 5.427452e-01 0.2392771700 5.219641e-01 0.0922133200 Amino acid biosynthesis_Histidine
19 B0216 2.297055e-01 0.1800818000 3.407278e-01 0.1741441000 Amino acid biosynthesis_Tryptophan
20 B0217 5.472435e-01 0.5034696000 4.691925e-01 0.4997109000 Amino acid biosynthesis_Phenylalanine
21 B0218 5.873705e-01 0.5161483000 5.094473e-01 0.4997109000 Amino acid biosynthesis_Tyrosine
22 B0219 1.003090e-01 0.0445315200 5.800346e-02 0.0000000000 Amino acid biosynthesis_GABA
23 B0220 1.016808e-01 0.0320511400 6.707939e-02 0.0000000000 Amino acid biosynthesis_Beta-alanine
24 B0221 5.496581e-01 0.6472659000 6.526429e-01 0.7491621000 Amino acid biosynthesis_Ornithine
25 B0303 3.466983e-01 0.3716875000 3.835668e-01 0.3997687000 Amino acid derivative biosynthesis_Ectoine
26 B0307 4.711377e-01 0.4760710000 4.454020e-01 0.4984230000 Amino acid derivative biosynthesis_Spermidine
27 B0309 4.463165e-02 0.0201724509 7.330752e-02 0.0005782321 Amino acid derivative biosynthesis_Putrescine
28 B0401 8.676715e-01 0.9596650000 9.260377e-01 0.9994218000 SCFA biosynthesis_Acetate
29 B0601 4.218585e-01 0.6429601000 5.411449e-01 0.7990212000 Organic anion biosynthesis_Succinate
30 B0602 7.492179e-01 0.7647096000 8.037116e-01 0.8003995000 Organic anion biosynthesis_Fumarate
31 B0603 5.372484e-01 0.4983080000 6.118297e-01 0.4983969000 Organic anion biosynthesis_Citrate
32 B0701 7.722910e-01 0.8843268000 8.451066e-01 1.0000000000 Vitamin biosynthesis_Thiamine (B1)
33 B0702 5.924357e-01 0.8247600000 7.562427e-01 0.9830361000 Vitamin biosynthesis_Riboflavin (B2)
34 B0704 4.176964e-01 0.7436086000 7.398774e-01 0.9990286000 Vitamin biosynthesis_Pantothenate (B5)
35 B0705 3.055590e-01 0.5240781000 4.311833e-01 0.6706592000 Vitamin biosynthesis_Pyridoxal-P (B6)
36 B0706 3.093418e-01 0.7369524000 4.730413e-01 0.9971351000 Vitamin biosynthesis_Biotin (B7)
37 B0707 6.927592e-01 0.8705010000 7.134942e-01 0.9930670000 Vitamin biosynthesis_Tetrahydrofolate (B9)
38 B0708 4.577213e-01 0.2544572000 2.930417e-01 0.1396510000 Vitamin biosynthesis_Cobalamin (B12)
39 B0709 5.759020e-02 0.0210456100 3.118951e-02 0.0000000000 Vitamin biosynthesis_Tocopherol/tocotorienol (E)
40 B0710 1.023563e-01 0.2501727000 1.276095e-01 0.3633975000 Vitamin biosynthesis_Phylloquinone (K1)
41 B0711 2.263965e-01 0.4743034000 3.552586e-01 0.6717029000 Vitamin biosynthesis_Menaquinone (K2)
42 B0712 1.006872e-01 0.1032767000 1.449259e-01 0.1169035000 Vitamin biosynthesis_Ubiquinone (Q10)
43 B0901 2.838361e-02 0.0095024940 7.457180e-03 0.0000000000 Metallophore biosynthesis_Staphyloferrin
44 B1026 6.711524e-04 0.0001491887 3.046125e-04 0.0000000000 Antibiotic biosynthesis_Prodigiosin
45 B1028 7.108709e-03 0.0045020230 6.407860e-03 0.0000000000 Antibiotic biosynthesis_Pyocyanin
46 B1041 1.362378e-04 0.0000000000 2.956816e-05 0.0000000000 Antibiotic biosynthesis_Validamycin A
47 D0101 3.670251e-01 0.4582496000 2.988741e-01 0.4828121000 Lipid degradation_Triglyceride
48 D0102 1.309692e-01 0.2029140000 1.325122e-01 0.2495006000 Lipid degradation_Fatty acid
49 D0309 6.354652e-01 0.5669771000 5.382669e-01 0.5078843000 Sugar degradation_Galactose
50 D0507 1.275300e-01 0.1495669000 1.738214e-01 0.1727479000 Amino acid degradation_Leucine
51 D0508 3.797487e-02 0.0184245300 3.756944e-02 0.0000000000 Amino acid degradation_Lysine
52 D0509 3.077201e-01 0.2725519000 2.734893e-01 0.2497897000 Amino acid degradation_Arginine
53 D0512 1.455705e-01 0.6568219000 3.220533e-01 0.9980025000 Amino acid degradation_Histidine
54 D0513 6.604162e-02 0.1082463800 7.564664e-02 0.1354670500 Amino acid degradation_Tryptophan
55 D0601 1.364125e-01 0.0386060100 9.505466e-02 0.0000000000 Nitrogen compound degradation_Nitrate
56 D0801 9.737248e-03 0.0033102250 1.973716e-02 0.0000000000 Xenobiotic degradation_Toluene
57 D0802 9.737248e-03 0.0033102250 1.973716e-02 0.0000000000 Xenobiotic degradation_Xylene
58 D0805 3.148056e-02 0.0104870100 1.012563e-02 0.0000000000 Xenobiotic degradation_Benzoate
59 D0807 1.322358e-01 0.0503759504 7.360588e-02 0.0001398254 Xenobiotic degradation_Catechol
60 D0809 5.113784e-03 0.0026462880 2.296552e-03 0.0000000000 Xenobiotic degradation_Biphenyl
61 D0811 1.669032e-02 0.0078700780 9.721016e-03 0.0000000000 Xenobiotic degradation_Benzoyl-CoA
62 D0812 8.995387e-05 0.0000000000 1.055236e-04 0.0000000000 Xenobiotic degradation_Naphthalene
63 D0816 1.013086e-01 0.0627671800 7.714537e-02 0.0391931200 Xenobiotic degradation_Phenylacetate
64 D0817 4.839446e-02 0.1183707700 7.698667e-02 0.1639595300 Xenobiotic degradation_Trans-cinnamate
65 D0901 3.219898e-01 0.7262781000 4.480076e-01 0.9988435000 Antibiotic degradation_Penicillin
66 D0902 3.985471e-02 0.0131677300 1.038850e-02 0.0000000000 Antibiotic degradation_Carbapenem
67 D0903 3.149326e-02 0.0124366700 1.023255e-02 0.0000000000 Antibiotic degradation_Cephalosporin
68 D0904 2.204813e-02 0.0219268639 7.383287e-04 0.0000000000 Antibiotic degradation_Oxacillin
69 D0905 2.706878e-01 0.1027584000 2.279679e-01 0.0000000000 Antibiotic degradation_Streptogramin
70 D0907 7.954688e-01 0.2816661860 7.049213e-01 0.0019975050 Antibiotic degradation_Tetracycline
71 D0908 4.457863e-01 0.1472018020 3.426434e-01 0.0019975050 Antibiotic degradation_Macrolide
72 D0910 2.304137e-01 0.0956614640 1.825602e-01 0.0011564640 Antibiotic degradation_Chloramphenicol
73 D0911 2.915679e-01 0.1096223000 2.437490e-01 0.0000000000 Antibiotic degradation_Lincosamide
74 D0912 3.863519e-02 0.0119847200 7.708563e-02 0.0000000000 Antibiotic degradation_Streptothricindifference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(gut_location) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) means_gift <- element_gift_filt %>%
select(-c(gut_location,fox_behaviour)) %>%
pivot_longer(!sample, names_to = "elements", values_to = "abundance") %>%
left_join(sample_metadata, by=join_by(sample==sample)) %>%
group_by(gut_location, elements) %>%
summarise(mean=mean(abundance))
# log_fold <- means_gift %>%
# group_by(elements) %>%
# summarise(
# logfc_high_low = log2(mean[gut_location == "high"] / mean[gut_location == "low"])
# )element_gift_names <- element_gift_filt %>%
select(-c(gut_location,fox_behaviour)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
select(-Elements)%>%
select(Function, everything())%>%
t()%>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(sample_metadata %>% select(sample,gut_location,fox_behaviour), by = join_by(sample == sample))
colNames <- names(element_gift_names %>% select(-c(sample,gut_location,fox_behaviour)))
for(i in colNames){
plt <- ggplot(element_gift_names, aes(x=gut_location, y=.data[[i]], color = gut_location, fill=gut_location)) +
geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
geom_jitter(width = 0.1, show.legend = TRUE) +
scale_color_manual(values=gut_location_colors)+
scale_fill_manual(values=gut_location_colors) +
theme_minimal() +
theme(
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank())
print(plt)
}
