Chapter 10 Functional differences
10.1 Metagenomics
genome_counts_filt <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0))%>%
rownames_to_column(., "genome")
genome_gifts <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts <- genome_gifts[, colSums(genome_gifts != 0) > 0]
genome_counts_filt <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),]
#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
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)
unique_funct_db <- GIFT_db[c(3, 4, 5)] %>%
distinct(Code_function, .keep_all = TRUE)10.1.1 Functional capacity of the MAGs
GIFTs_elements %>%
reshape2::melt() %>%
dplyr::rename(Genome = Var1, Code_element = Var2, GIFT = value) %>%
inner_join(GIFT_db,by="Code_element") %>%
ggplot(., aes(x=Code_element, y=Genome, fill=GIFT, group=Code_function))+
geom_tile()+
scale_y_discrete(guide = guide_axis(check.overlap = TRUE))+
scale_x_discrete(guide = guide_axis(check.overlap = TRUE))+
scale_fill_gradientn(limits = c(0,1), colours=brewer.pal(7, "YlGnBu"))+
facet_grid(. ~ Code_function, scales = "free", space = "free")+
theme_grey(base_size=8)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text.x = element_text(angle = 90))
10.1.2 MCI
GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
left_join(sample_metadata%>% select(sample, Species), by = join_by(sample == sample)) %>%
left_join(unique_funct_db[c(1,2)], by = join_by(trait == Code_function)) %>%
group_by(Species, Domain) %>%
summarise(mean = mean(value), sd = sd(value)) %>%
mutate(MCI = str_c(round(mean, 3), "±", round(sd, 3))) %>%
select(-mean, -sd) %>%
pivot_wider(names_from = Domain,
values_from = MCI) %>%
tt()| Species | Biosynthesis | Degradation | Structure |
|---|---|---|---|
| Eb | 0.449±0.31 | 0.225±0.167 | 0.342±0.318 |
| Ha | 0.427±0.326 | 0.19±0.17 | 0.321±0.332 |
| Pk | 0.42±0.326 | 0.236±0.197 | 0.288±0.308 |
unique_funct_db <- GIFT_db[c(3, 4, 5)] %>%
distinct(Code_function, .keep_all = TRUE)
MCI <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
left_join(sample_metadata%>% select(sample, Species), by = join_by(sample == sample)) %>%
left_join(unique_funct_db[c(1,2)], by = join_by(trait == Code_function))
# Biosynthesis
MCI_sub <- MCI %>% filter(Domain=="Biosynthesis")
shapiro.test(MCI_sub$value)
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.91112, p-value = 9.188e-16
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.49824, df = 2, p-value = 0.7795
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.8995, p-value = 9.34e-16
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 5.0406, df = 2, p-value = 0.08044
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.80948, p-value = 7.71e-13
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 2.0268, df = 2, p-value = 0.36310.1.3 Differences in functional traits considering the three species
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(1,2)], by=join_by("sample"=="sample"))element_gift %>%
pivot_longer(-c(sample,Species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = summary(aov(value ~ Species))[[1]][["Pr(>F)"]][1]) %>%
filter(p_value < 0.05) %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))# A tibble: 19 × 3
trait p_value Function
<chr> <dbl> <chr>
1 B0104 0.0285 Nucleic acid biosynthesis_CDP/CTP
2 B0309 0.0216 Amino acid derivative biosynthesis_Putrescine
3 B0310 0.0302 Amino acid derivative biosynthesis_Tryptamine
4 D0205 0.00797 Polysaccharide degradation_Pectin
5 D0211 0.00451 Polysaccharide degradation_Alpha-mannan
6 D0306 0.0199 Sugar degradation_D-Xylose
7 D0308 0.0487 Sugar degradation_L-Rhamnose
8 D0310 0.00000333 Sugar degradation_NeuAc
9 D0517 0.0266 Amino acid degradation_Ornithine
10 D0604 0.0124 Nitrogen compound degradation_GlcNAc
11 D0612 0.0111 Nitrogen compound degradation_Hypotaurine
12 D0613 0.0389 Nitrogen compound degradation_Taurine
13 D0706 0.0460 Alcohol degradation_Ethylene glycol
14 D0801 0.0395 Xenobiotic degradation_Toluene
15 D0802 0.0395 Xenobiotic degradation_Xylene
16 D0905 0.0382 Antibiotic degradation_Streptogramin
17 D0911 0.0427 Antibiotic degradation_Lincosamide
18 S0104 0.0190 Cellular structure_Lipoteichoic acid
19 S0105 0.0298 Cellular structure_Lipopolysaccharide Pairwise comparison
anova_results <- element_gift %>%
pivot_longer(-c(sample, Species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(anova_p = summary(aov(value ~ Species))[[1]][["Pr(>F)"]][1] ) %>%
# mutate(p_adjust_ANOVA = p.adjust(anova_p, method = "BH")) %>%
filter(anova_p < 0.05)
element_gift %>%
pivot_longer(-c(sample, Species), names_to = "trait", values_to = "value") %>%
filter(trait %in% anova_results$trait) %>%
group_by(trait) %>%
summarise(
tukey_results = list(
as.data.frame(TukeyHSD(aov(value ~ Species))$Species) |>
tibble::rownames_to_column("contrast")
)
) %>%
unnest(tukey_results) %>%
filter(`p adj` < 0.05) %>%
left_join(., uniqueGIFT_db[c(1,3)], by = join_by(trait == Code_element))# A tibble: 20 × 7
trait contrast diff lwr upr `p adj` Function
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 B0104 Pk-Ha -0.138 -0.258 -0.0175 0.0213 Nucleic acid biosynthesis_CDP/CTP
2 B0309 Pk-Ha -0.156 -0.290 -0.0214 0.0196 Amino acid derivative biosynthesis_Putrescine
3 B0310 Ha-Eb -0.000528 -0.00101 -0.0000467 0.0285 Amino acid derivative biosynthesis_Tryptamine
4 D0205 Ha-Eb -0.0837 -0.155 -0.0119 0.0188 Polysaccharide degradation_Pectin
5 D0205 Pk-Ha 0.0646 0.00702 0.122 0.0246 Polysaccharide degradation_Pectin
6 D0211 Ha-Eb -0.0491 -0.0918 -0.00646 0.0205 Polysaccharide degradation_Alpha-mannan
7 D0211 Pk-Ha 0.0436 0.00939 0.0778 0.00935 Polysaccharide degradation_Alpha-mannan
8 D0306 Pk-Ha 0.201 0.0344 0.369 0.0146 Sugar degradation_D-Xylose
9 D0310 Pk-Eb 0.178 0.0209 0.335 0.0229 Sugar degradation_NeuAc
10 D0310 Pk-Ha 0.305 0.176 0.434 0.00000193 Sugar degradation_NeuAc
11 D0517 Pk-Ha -0.192 -0.364 -0.0210 0.0245 Amino acid degradation_Ornithine
12 D0604 Pk-Ha 0.329 0.0719 0.586 0.00906 Nitrogen compound degradation_GlcNAc
13 D0612 Pk-Ha 0.0700 0.0160 0.124 0.00814 Nitrogen compound degradation_Hypotaurine
14 D0613 Pk-Ha -0.229 -0.441 -0.0166 0.0319 Nitrogen compound degradation_Taurine
15 D0801 Pk-Ha 0.0319 0.00245 0.0614 0.0310 Xenobiotic degradation_Toluene
16 D0802 Pk-Ha 0.0319 0.00245 0.0614 0.0310 Xenobiotic degradation_Xylene
17 D0905 Pk-Ha 0.121 0.00456 0.237 0.0400 Antibiotic degradation_Streptogramin
18 D0911 Pk-Ha 0.140 0.00512 0.275 0.0403 Antibiotic degradation_Lincosamide
19 S0104 Pk-Ha 0.110 0.0187 0.201 0.0146 Cellular structure_Lipoteichoic acid
20 S0105 Pk-Ha -0.228 -0.430 -0.0260 0.0235 Cellular structure_Lipopolysaccharide 10.2 Metabarcoding Standard
10.2.1 Prepare for PicrusT
write_tsv(genome_counts_filt,"resources/amplicon/picrust/counts.tsv")
asvs <- genome_metadata$genome
fasta_all <- readDNAStringSet("resources/amplicon/ASVs_seqtab_overlap10_nochim0.fa", format="fasta")
fasta_filtered <- fasta_all[asvs]
writeXStringSet(fasta_filtered,
filepath = "resources/amplicon/picrust/sequences.fa",
format = "fasta")10.2.3 Analyse functions
amplicon_ko <- read_tsv("resources/amplicon/picrust/results/combined_KO_predicted.tsv.xz")
amplicon_ec <- read_tsv("resources/amplicon/picrust/results/combined_EC_predicted.tsv.xz")
tree_bacteria <- read.tree("resources/amplicon/picrust/results/bac_reduced.tre")
save(amplicon_ko,amplicon_ec,tree_bacteria, file="resources/functional_files.Rdata")amplicon_ko_long <- amplicon_ko %>%
pivot_longer(!sequence,names_to="kegg",values_to="value") %>%
filter(value>0) %>%
mutate(kegg = stringr::str_remove(kegg, "^ko:"))
amplicon_ec_long <- amplicon_ec %>%
pivot_longer(!sequence,names_to="ec",values_to="value") %>%
filter(value>0) %>%
mutate(ec = str_c("[EC:",ec,"]"))
# Distil KO and EC separately (for later merging)
amplicon_ko_distil <- distill(amplicon_ko_long,GIFT_db,genomecol=1,annotcol=2, verbosity=T)
amplicon_ec_distil <- distill(amplicon_ec_long,GIFT_db,genomecol=1,annotcol=2, verbosity=T)
# Save distilled tables
amplicon_ko_distil %>%
as.data.frame() %>%
rownames_to_column("genome") %>%
write_tsv("resources/amplicon/picrust/results/gifts_ko.tsv.xz")
amplicon_ec_distil %>%
as.data.frame() %>%
rownames_to_column("genome") %>%
write_tsv("resources/amplicon/picrust/results/gifts_ec.tsv.xz")# Merge distilled tables
amplicon_ko_distil <- read_tsv("resources/amplicon/picrust/results/gifts_ko.tsv.xz") %>% column_to_rownames(var="genome")
amplicon_ec_distil <- read_tsv("resources/amplicon/picrust/results/gifts_ec.tsv.xz") %>% column_to_rownames(var="genome")
amplicon_distil <- pmax(amplicon_ko_distil, amplicon_ec_distil)
amplicon_elements <- to.elements(amplicon_distil,GIFT_db)
#Aggregate element-level GIFTs into the function level
amplicon_functions <- to.functions(amplicon_elements,GIFT_db)
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
amplicon_elem_community_standard <- to.community(amplicon_elements,genome_counts_row,GIFT_db)
amplicon_func_community_standard <- to.community(amplicon_functions,genome_counts_row,GIFT_db)10.2.4 MCI
amplicon_func_community_standard %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
left_join(sample_metadata%>% select(sample, Species), by = join_by(sample == sample)) %>%
left_join(unique_funct_db[c(1,2)], by = join_by(trait == Code_function)) %>%
group_by(Species, Domain) %>%
summarise(mean = mean(value), sd = sd(value)) %>%
mutate(MCI = str_c(round(mean, 3), "±", round(sd, 3))) %>%
select(-mean, -sd) %>%
pivot_wider(names_from = Domain,
values_from = MCI) %>%
tt()| Species | Biosynthesis | Degradation | Structure |
|---|---|---|---|
| Eb | 0.703±0.517 | 0.335±0.22 | 0.416±0.218 |
| Ha | 0.704±0.533 | 0.305±0.195 | 0.435±0.21 |
| Pk | 0.702±0.528 | 0.316±0.217 | 0.427±0.205 |
MCI <- amplicon_func_community_standard %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
left_join(sample_metadata%>% select(sample, Species), by = join_by(sample == sample)) %>%
left_join(unique_funct_db[c(1,2)], by = join_by(trait == Code_function))
#Biosynthesis
MCI_sub <- MCI %>% filter(Domain=="Biosynthesis")
shapiro.test(MCI_sub$value)
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.90918, p-value = 6.041e-16
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.10533, df = 2, p-value = 0.9487
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.96435, p-value = 4.1e-09
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.99844, df = 2, p-value = 0.607
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.95961, p-value = 0.0001922
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.18124, df = 2, p-value = 0.9134amplicon_func_community_standard %>%
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 ~ Species, 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")
10.2.5 Differences in functional traits considering the three species
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift_MS <- amplicon_elem_community_standard %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c("sample","Species")], by=join_by("sample"=="sample"))element_gift_MS %>%
pivot_longer(-c(sample,Species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = summary(aov(value ~ Species))[[1]][["Pr(>F)"]][1]) %>%
filter(p_value < 0.05) %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))# A tibble: 8 × 3
trait p_value Function
<chr> <dbl> <chr>
1 B0210 0.000138 Amino acid biosynthesis_Leucine
2 B0218 0.0328 Amino acid biosynthesis_Tyrosine
3 B1006 0.000469 Antibiotic biosynthesis_Carbapenem-3-carboxylate
4 D0306 0.0154 Sugar degradation_D-Xylose
5 D0508 0.0367 Amino acid degradation_Lysine
6 D0604 0.0384 Nitrogen compound degradation_GlcNAc
7 D0609 0.0341 Nitrogen compound degradation_L-carnitine
8 D0910 0.00808 Antibiotic degradation_Chloramphenicol Pairwise comparison
anova_results_MS <- element_gift_MS %>%
pivot_longer(-c(sample, Species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(anova_p = summary(aov(value ~ Species))[[1]][["Pr(>F)"]][1] ) %>%
# mutate(p_adjust_ANOVA = p.adjust(anova_p, method = "BH")) %>%
filter(anova_p < 0.05)
element_gift_MS %>%
pivot_longer(-c(sample, Species), names_to = "trait", values_to = "value") %>%
filter(trait %in% anova_results_MS$trait) %>%
group_by(trait) %>%
summarise(
tukey_results = list(
as.data.frame(TukeyHSD(aov(value ~ Species))$Species) |>
tibble::rownames_to_column("contrast")
)
) %>%
unnest(tukey_results) %>%
filter(`p adj` < 0.05) %>%
left_join(., uniqueGIFT_db[c(1,3)], by = join_by(trait == Code_element))# A tibble: 12 × 7
trait contrast diff lwr upr `p adj` Function
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 B0210 Ha-Eb -0.226 -0.344 -0.107 0.0000869 Amino acid biosynthesis_Leucine
2 B0210 Pk-Eb -0.131 -0.247 -0.0154 0.0230 Amino acid biosynthesis_Leucine
3 B0218 Ha-Eb 0.146 0.0148 0.276 0.0259 Amino acid biosynthesis_Tyrosine
4 B1006 Ha-Eb -0.00141 -0.00230 -0.000513 0.00114 Antibiotic biosynthesis_Carbapenem-3-carboxylate
5 B1006 Pk-Eb -0.00143 -0.00230 -0.000555 0.000723 Antibiotic biosynthesis_Carbapenem-3-carboxylate
6 D0306 Ha-Eb 0.158 0.0101 0.305 0.0338 Sugar degradation_D-Xylose
7 D0306 Pk-Ha -0.123 -0.242 -0.00506 0.0392 Sugar degradation_D-Xylose
8 D0508 Ha-Eb -0.0721 -0.139 -0.00561 0.0308 Amino acid degradation_Lysine
9 D0604 Ha-Eb 0.175 0.0121 0.338 0.0327 Nitrogen compound degradation_GlcNAc
10 D0609 Pk-Eb 0.0798 0.00570 0.154 0.0321 Nitrogen compound degradation_L-carnitine
11 D0910 Pk-Eb -0.218 -0.436 -0.00116 0.0485 Antibiotic degradation_Chloramphenicol
12 D0910 Pk-Ha -0.217 -0.396 -0.0387 0.0136 Antibiotic degradation_Chloramphenicol 10.2.6 Metabarcoding selected
load("resources/amplicon/selected_filtering.Rdata")
genome_counts_row <- genome_counts_filt_amplicon %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
amplicon_elements_filtered <- amplicon_elements[rownames(genome_counts_row),]
#Aggregate element-level GIFTs into the function level
amplicon_functions <- to.functions(amplicon_elements_filtered,GIFT_db)
amplicon_elem_community_selected <- to.community(amplicon_elements_filtered,genome_counts_row,GIFT_db)
amplicon_func_community_selected <- to.community(amplicon_functions,genome_counts_row,GIFT_db)amplicon_elem_community_selected %>%
as.data.frame() %>%
rownames_to_column(var="sample")%>%
pivot_longer(!sample,names_to="trait",values_to="gift") %>%
left_join(sample_metadata_amplicon, 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 ~ Species, 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")
10.2.7 MCI
amplicon_func_community_selected %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
left_join(sample_metadata%>% select(sample, Species), by = join_by(sample == sample)) %>%
left_join(unique_funct_db[c(1,2)], by = join_by(trait == Code_function)) %>%
group_by(Species, Domain) %>%
summarise(mean = mean(value), sd = sd(value)) %>%
mutate(MCI = str_c(round(mean, 3), "±", round(sd, 3))) %>%
select(-mean, -sd) %>%
pivot_wider(names_from = Domain,
values_from = MCI)# A tibble: 3 × 4
# Groups: Species [3]
Species Biosynthesis Degradation Structure
<chr> <chr> <chr> <chr>
1 Eb 0.556±0.328 0.308±0.215 0.468±0.294
2 Ha 0.562±0.339 0.303±0.218 0.485±0.278
3 Pk 0.553±0.353 0.317±0.23 0.449±0.273unique_funct_db <- GIFT_db[c(3, 4, 5)] %>%
distinct(Code_function, .keep_all = TRUE)
MCI <- amplicon_func_community_selected %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
pivot_longer(-sample, names_to = "trait", values_to = "value") %>%
left_join(sample_metadata%>% select(sample, Species), by = join_by(sample == sample)) %>%
left_join(unique_funct_db[c(1,2)], by = join_by(trait == Code_function))
#Biosynthesis
MCI_sub <- MCI %>% filter(Domain=="Biosynthesis")
shapiro.test(MCI_sub$value)
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.87076, p-value < 2.2e-16
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.020347, df = 2, p-value = 0.9899
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.94832, p-value = 7.699e-11
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.25591, df = 2, p-value = 0.8799
Shapiro-Wilk normality test
data: MCI_sub$value
W = 0.90306, p-value = 5.54e-08
Kruskal-Wallis rank sum test
data: value by Species
Kruskal-Wallis chi-squared = 0.37772, df = 2, p-value = 0.827910.2.8 Differences in functional traits considering the three species
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift_MSel <- amplicon_elem_community_selected %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c("sample","Species")], by=join_by("sample"=="sample"))element_gift_MSel %>%
pivot_longer(-c(sample,Species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = summary(aov(value ~ Species))[[1]][["Pr(>F)"]][1]) %>%
filter(p_value < 0.05) %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))# A tibble: 10 × 3
trait p_value Function
<chr> <dbl> <chr>
1 B0402 0.0250 SCFA biosynthesis_Butyrate
2 B0901 0.0268 Metallophore biosynthesis_Staphyloferrin
3 B1011 0.00223 Antibiotic biosynthesis_Erythromycin
4 D0502 0.00893 Amino acid degradation_Threonine
5 D0505 0.0206 Amino acid degradation_Valine
6 D0507 0.00164 Amino acid degradation_Leucine
7 D0518 0.0424 Amino acid degradation_GABA
8 D0607 0.0446 Nitrogen compound degradation_Creatinine
9 D0808 0.0389 Xenobiotic degradation_Cumate
10 S0103 0.0427 Cellular structure_Teichoic acid Pairwise comparison
anova_results_MSel <- element_gift_MSel %>%
pivot_longer(-c(sample, Species), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(anova_p = summary(aov(value ~ Species))[[1]][["Pr(>F)"]][1] ) %>%
# mutate(p_adjust_ANOVA = p.adjust(anova_p, method = "BH")) %>%
filter(anova_p < 0.05)
element_gift_MSel %>%
pivot_longer(-c(sample, Species), names_to = "trait", values_to = "value") %>%
filter(trait %in% anova_results_MSel$trait) %>%
group_by(trait) %>%
summarise(
tukey_results = list(
as.data.frame(TukeyHSD(aov(value ~ Species))$Species) |>
tibble::rownames_to_column("contrast")
)
) %>%
unnest(tukey_results) %>%
filter(`p adj` < 0.05) %>%
left_join(., uniqueGIFT_db[c(1,3)], by = join_by(trait == Code_element))# A tibble: 11 × 7
trait contrast diff lwr upr `p adj` Function
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 B0402 Pk-Ha 0.0862 0.00801 0.164 0.0278 SCFA biosynthesis_Butyrate
2 B0901 Ha-Eb -0.00592 -0.0116 -0.000255 0.0388 Metallophore biosynthesis_Staphyloferrin
3 B0901 Pk-Eb -0.00608 -0.0119 -0.000295 0.0375 Metallophore biosynthesis_Staphyloferrin
4 B1011 Ha-Eb -0.102 -0.171 -0.0335 0.00226 Antibiotic biosynthesis_Erythromycin
5 B1011 Pk-Eb -0.0902 -0.160 -0.0200 0.00884 Antibiotic biosynthesis_Erythromycin
6 D0502 Ha-Eb 0.288 0.0678 0.508 0.00764 Amino acid degradation_Threonine
7 D0505 Pk-Eb 0.121 0.0163 0.225 0.0200 Amino acid degradation_Valine
8 D0507 Ha-Eb 0.104 0.0219 0.187 0.0100 Amino acid degradation_Leucine
9 D0507 Pk-Ha -0.0981 -0.168 -0.0278 0.00426 Amino acid degradation_Leucine
10 D0808 Pk-Eb -0.000696 -0.00138 -0.00000735 0.0471 Xenobiotic degradation_Cumate
11 S0103 Pk-Eb 0.103 0.00630 0.200 0.0346 Cellular structure_Teichoic acid