DIABIMMUNE Cohort
Visualizing antibiotics exposures in the DIABIMMUNE ABx Dataset
Load-in required data
Data were processed by Nephele
phyloYassour <- import_biom("../data/Yassour/taxa.biom")
metaYassour <- read.csv("../data/Yassour/YassourM_2016_metadata.csv")
treeYassour <- read_tree("../data/Yassour/phylo/rooted_tree.nwk")
# One of the metadata rows is missing and should be removed
metaYass <- sample_data(metaYassour)
metaYass$SampleID <- metaYass$X.SampleID
metaYass <- metaYass[!is.na(metaYass$SampleID)]
sample_names(metaYass) <- metaYass$SampleID
sample_data(phyloYassour) <- metaYass
yass <- phyloYassour
phy_tree(yass) <- treeYassourCalculate Beta diversity at each time point
Let’s make a function!
betaDivPlot <- function(x){
dd <- sample_data(yass)$study_condition
dd.col <- rainbow(length(dd))
names(dd.col) <- dd
subphylo <- prune_samples(sample_data(yass)$subject_id==x, x=yass)
bray <- distance(subphylo, method="bray") %>% as.matrix() %>% as.data.frame()
bray$sample_id <- colnames(bray)
phylosamp <- sample_data(yass) %>% data.frame()
joined <- inner_join (phylosamp, bray)
minsamp <- joined %>% filter(infant_age==min(infant_age)) %>% select(sample_id) %>% .[1,1]
joinsub <- joined %>% select(sample_id, infant_age, study_condition, braycurtis=all_of(minsamp)) %>% arrange(infant_age)
joinsub$study_condition[is.na(joinsub$study_condition)] <- "no abx"
plot <- ggplot(joinsub, aes(x=infant_age, y=braycurtis))+geom_line()+geom_point(aes(color=factor(study_condition))) +
scale_color_manual("Legend", values = dd.col)+ theme_minimal() +theme(legend.position="none") +
theme(axis.text.x=element_blank(), #remove x axis labels
axis.ticks.x=element_blank(), #remove x axis ticks
axis.text.y=element_blank(), #remove y axis labels
axis.ticks.y=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank()#remove y axis ticks
)
return(plot)
}Test it out
betaDivPlot("E010481")Joining, by = "sample_id"
subjects <- yass %>% sample_data() %>% data.frame() %>% arrange(desc(study_condition)) %>% .$subject_id %>% unique()
listPlots <- lapply(subjects, betaDivPlot)Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"n <- length(listPlots)
nCol <- floor(sqrt(n))
do.call("grid.arrange", c(listPlots, ncol=nCol))
Let’s try it while removing early days which might be introducing a different effect
yassleft <- subset_samples(yass, infant_age>=180)betaDivPlotLeft <- function(x){
dd <- sample_data(yassleft)$study_condition
dd.col <- rainbow(length(dd))
names(dd.col) <- dd
subphylo <- prune_samples(sample_data(yassleft)$subject_id==x, x=yassleft)
bray <- distance(subphylo, method="bray") %>% as.matrix() %>% as.data.frame()
bray$sample_id <- colnames(bray)
phylosamp <- sample_data(yassleft) %>% data.frame()
joined <- inner_join (phylosamp, bray)
minsamp <- joined %>% filter(infant_age==min(infant_age)) %>% select(sample_id) %>% .[1,1]
joinsub <- joined %>% select(sample_id, infant_age, study_condition, braycurtis=all_of(minsamp)) %>% arrange(infant_age)
joinsub$study_condition[is.na(joinsub$study_condition)] <- "no abx"
plot <- ggplot(joinsub, aes(x=infant_age, y=braycurtis))+geom_line()+geom_point(aes(color=factor(study_condition))) +
scale_color_manual("Legend", values = dd.col)+ theme_minimal() +theme(legend.position="none") +
theme(axis.text.x=element_blank(), #remove x axis labels
axis.ticks.x=element_blank(), #remove x axis ticks
axis.text.y=element_blank(), #remove y axis labels
axis.ticks.y=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank()#remove y axis ticks
)
return(plot)
}subjectsleft <- yassleft %>% sample_data() %>% data.frame() %>% arrange(desc(study_condition)) %>% .$subject_id %>% unique()
listPlots <- lapply(subjectsleft, betaDivPlotLeft)Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"n <- length(listPlots)
nCol <- floor(sqrt(n))
do.call("grid.arrange", c(listPlots, ncol=nCol))
That didn’t seem to do much. Let’s average across study conditions and see if there’s a difference
Calculate dissimilarity over time for each subject
Determine study condition and add to sample_data()
library(microViz)
microViz version 0.9.3 - Copyright (C) 2021 David Barnett
* Website: https://david-barnett.github.io/microViz/
* Useful? For citation info, run: citation('microViz')
* Silence: suppressPackageStartupMessages(library(microViz))studyconds <- sample_data(yass) %>% group_by(subject_id) %>% summarize(subject_id=dplyr::first(subject_id), abxexpose = n_distinct(study_condition)>1)
yass <- ps_join(yass, studyconds)Joining, by = "subject_id"betaValues <- function(x){
dd <- sample_data(yass)$study_condition
dd.col <- rainbow(length(dd))
names(dd.col) <- dd
subphylo <- prune_samples(sample_data(yass)$subject_id==x, x=yass)
bray <- distance(subphylo, method="bray") %>% as.matrix() %>% as.data.frame()
bray$sample_id <- colnames(bray)
phylosamp <- sample_data(yass) %>% data.frame()
joined <- inner_join (phylosamp, bray)
minsamp <- joined %>% filter(infant_age==max(infant_age)) %>% select(sample_id) %>% .[1,1]
joinsub <- joined %>% select(subject_id, sample_id, infant_age, agecat, study_condition, abxexpose, braycurtis=all_of(minsamp)) %>% arrange(infant_age)
joinsub$study_condition[is.na(joinsub$study_condition)] <- "no abx"
return(joinsub)
}sample_data(yass)$agecat[sample_data(yass)$infant_age<100] <- "0-99 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=100 & sample_data(yass)$infant_age < 200] <- "100-199 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=200 & sample_data(yass)$infant_age < 300] <- "200-299 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=300 & sample_data(yass)$infant_age < 400] <- "300-399 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=400 & sample_data(yass)$infant_age < 500] <- "400-499 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=500 & sample_data(yass)$infant_age < 600] <- "500-599 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=600 & sample_data(yass)$infant_age < 700] <- "600-699 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=700 & sample_data(yass)$infant_age < 800] <- "700-799 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=800 & sample_data(yass)$infant_age < 900] <- "800-899 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=900 & sample_data(yass)$infant_age < 1000] <- "900-999 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=1000 & sample_data(yass)$infant_age < 1100] <- "x1000-1099 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=1100 & sample_data(yass)$infant_age < 1200] <- "x1100-1199 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=1200 & sample_data(yass)$infant_age < 1300] <- "x1200-1299 days"listValues <- lapply(subjects, betaValues)Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"
Joining, by = "sample_id"dfValues <- as.data.frame(do.call(rbind, listValues)) %>% arrange(infant_age)
plotValues <- dfValues %>% group_by(abxexpose, agecat) %>% summarize(bc=mean(braycurtis))`summarise()` has grouped output by 'abxexpose'. You can override using the
`.groups` argument.ggplot(data=plotValues, aes(x=agecat, y=bc)) + geom_line(aes(group=abxexpose, color=abxexpose))+ scale_x_discrete(labels = function(x) str_wrap(x, width = 0.01))
Calculate average dissimilarity for each time category
Select a time point then calculate average within participant dissimilarity within that timepoint
t1cases <- prune_samples(sample_data(yass)$agecat=="800-899 days" & sample_data(yass)$abxexpose==TRUE, x=yass)
brayt1cases <- distance(t1cases, method="bray") %>% as.matrix()
brayt1cases %>% mean()[1] 0.8447845avgBetaDivPlot <- function(x){
dd <- sample_data(yass)$study_condition
dd.col <- rainbow(length(dd))
names(dd.col) <- dd
subphylo <- prune_samples(sample_data(yass)$subject_id==x, x=yass)
bray <- distance(subphylo, method="bray") %>% as.matrix() %>% as.data.frame()
bray$sample_id <- colnames(bray)
phylosamp <- sample_data(yass) %>% data.frame()
joined <- inner_join (phylosamp, bray)
minsamp <- joined %>% filter(infant_age==min(infant_age)) %>% select(sample_id) %>% .[1,1]
joinsub <- joined %>% select(sample_id, infant_age, study_condition, braycurtis=all_of(minsamp)) %>% arrange(infant_age)
joinsub$study_condition[is.na(joinsub$study_condition)] <- "no abx"
return(joinsub)
}avgBetaDivPlot("E010481")Joining, by = "sample_id" sample_id infant_age
1 SRR8115708 63
2 SRR8117230 63
3 SRR8119454 78
4 SRR8119450 121
5 SRR8119453 150
6 SRR8119182 178
7 SRR8119448 207
8 SRR8119187 239
9 SRR8119353 272
10 SRR8119455 295
11 SRR8119449 327
12 SRR8116363 366
13 SRR8117118 366
14 SRR8119451 420
15 SRR8118993 463
16 SRR8118997 482
17 SRR8119181 518
18 SRR8119350 547
19 SRR8119000 577
20 SRR8119188 612
21 SRR8119716 628
22 SRR8119186 668
23 SRR8119185 693
24 SRR8116361 719
25 SRR8117120 719
26 SRR8119696 778
27 SRR8119698 819
28 SRR8119452 846
29 SRR8119722 877
30 SRR8115717 907
31 SRR8117233 907
32 SRR8115716 942
33 SRR8117232 942
34 SRR8119697 971
35 SRR8115709 1003
36 SRR8117231 1003
37 SRR8116362 1034
38 SRR8117117 1034
39 SRR8116368 1067
40 SRR8117119 1067
41 SRR8116364 1089
42 SRR8117234 1089
study_condition braycurtis
1 Amoxicillin 0.0000000
2 no abx 0.0000000
3 no abx 0.9649608
4 no abx 0.7099739
5 Amoxicillin and clavulanic acid 0.9690648
6 Amoxicillin and clavulanic acid 0.9377355
7 Azithromycin 0.9241796
8 no abx 0.9442477
9 Cefalexin, Amoxicillin and clavulanic acid 0.9551169
10 no abx 0.9650799
11 no abx 0.9879075
12 no abx 0.9977106
13 no abx 0.9977106
14 no abx 0.9699916
15 no abx 0.8925375
16 no abx 0.9806153
17 Amoxicillin and clavulanic acid 0.9791326
18 no abx 0.9347861
19 no abx 0.9869364
20 no abx 0.9804419
21 no abx 0.9728718
22 no abx 0.9658010
23 Amoxicillin 0.9634368
24 no abx 0.9713640
25 no abx 0.9713640
26 Trimetoprime and sulfadiazine 1.0000000
27 no abx 0.9424588
28 Amoxicillin and clavulanic acid 0.8894468
29 Trimetoprime and sulfadiazine, Trimetoprime and sulfadiazine 1.0000000
30 Amoxicillin 0.8821137
31 no abx 0.8821137
32 no abx 0.9289515
33 no abx 0.9289515
34 no abx 0.9307430
35 no abx 0.9356681
36 no abx 0.9356681
37 Amoxicillin 1.0000000
38 no abx 1.0000000
39 no abx 0.9139071
40 no abx 0.9139071
41 no abx 1.0000000
42 no abx 1.0000000PCOA
Simplify antibiotics exposures
sample_data(yass)$abx <- "none"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin"] <- "penicillins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefalexin"] <- "cephalosporins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin and clavulanic acid"] <- "penicillins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Azithromycin"] <- "macrolides"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Amoxicillin and clavulanic acid"] <- "penicillins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin and clavulanic acid, Amoxicillin and clavulanic acid"] <- "penicillins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Trimetoprime and sulfadiazine, Amoxicillin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Trimetoprime and sulfadiazine"] <- "sulfonamides"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Cefaclor, Azithromycin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Azithromycin, Cefalexin, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Azithromycin, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Cefalexin, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefaclor"] <- "cephalosporins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin and clavulanic acid, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Phenoxymethylpenicillin"] <- "penicillins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Clarithromycin"] <- "macrolides"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Trimetoprime and sulfadiazine, Trimetoprime and sulfadiazine"] <- "sulfonamides"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Trimetoprime and sulfadiazine, Amoxicillin and clavulanic acid"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin and clavulanic acid, Azithromycin, Amoxicillin and clavulanic acid"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Clarithromycin, Cefalexin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefalexin, Amoxicillin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Ceftriaxone, Ceftriaxone"] <- "cephalosporins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Azithromycin, Amoxicillin and clavulanic acid" ] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Penicillin G, Netilmicin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Clarithromycin, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Ceftriaxone"] <- "cephalosporins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin and clavulanic acid, Trimetoprime and sulfadiazine, Cefalexin, Amoxicillin and clavulanic acid"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Amoxicillin and clavulanic acid, Amoxicillin, Amoxicillin and clavulanic acid, Trimetoprime and sulfadiazine, Azithromycin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Azithromycin, Azithromycin"] <- "macrolides"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Amoxicillin"] <- "penicillins"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefalexin, Amoxicillin and clavulanic acid"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefalexin, Cefalexin, Azithromycin, Amoxicillin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Trimetoprime and sulfadiazine, Azithromycin, Cefaclor"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Systemic antibiotic NAS"] <- "other"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin and clavulanic acid, Azithromycin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefalexin, Cefalexin, Trimetoprime and sulfadiazine"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Systemic antibiotic NAS" ] <- "other"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Trimetoprime and sulfadiazine, Trimetoprime and sulfadiazine, Amoxicillin, Amoxicillin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Amoxicillin, Trimetoprime and sulfadiazine, Amoxicillin"] <- "multi"
sample_data(yass)$abx[sample_data(yass)$study_condition=="Cefaclor, Cefalexin"] <- "other"#ordu = ordinate(yass, "PCoA", "unifrac", weighted=TRUE)
#plot_ordination(yass, ordu, color="abx", shape="human")
#plot_ordination(yass, ordu, color="infant_age", shape="human")Individual PCOA plots
pcoaPlot <- function(x){
dd <- sample_data(yass)$study_condition
dd.col <- rainbow(length(dd))
names(dd.col) <- dd
subphylo <- prune_samples(sample_data(yass)$subject_id==x, x=yass)
ordu <- ordinate(subphylo, "PCoA", "unifrac", weighted=TRUE)
plot_ordination(yass, ordu, color="abx", label="infant_age")
}subjects <- yass %>% sample_data() %>% data.frame() %>% arrange(desc(study_condition)) %>% .$subject_id %>% unique()
lapply(subjects, pcoaPlot)[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
[[19]]
[[20]]
[[21]]
[[22]]
[[23]]
[[24]]
[[25]]
[[26]]
[[27]]
[[28]]
[[29]]
[[30]]
[[31]]
[[32]]
[[33]]
[[34]]
[[35]]
[[36]]
[[37]]
[[38]]
Group by age
sample_data(yass)$agecat[sample_data(yass)$infant_age<100] <- "0-99 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=100 & sample_data(yass)$infant_age < 200] <- "100-199 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=200 & sample_data(yass)$infant_age < 300] <- "200-299 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=300 & sample_data(yass)$infant_age < 400] <- "300-399 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=400 & sample_data(yass)$infant_age < 500] <- "400-499 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=500 & sample_data(yass)$infant_age < 600] <- "500-599 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=600 & sample_data(yass)$infant_age < 700] <- "600-699 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=700 & sample_data(yass)$infant_age < 800] <- "700-799 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=800 & sample_data(yass)$infant_age < 900] <- "800-899 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=900 & sample_data(yass)$infant_age < 1000] <- "900-999 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=1000 & sample_data(yass)$infant_age < 1100] <- "x1000-1099 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=1100 & sample_data(yass)$infant_age < 1200] <- "x1100-1199 days"
sample_data(yass)$agecat[sample_data(yass)$infant_age>=1200 & sample_data(yass)$infant_age < 1300] <- "x1200-1299 days"Function to plot those abx ordinations but stratified by age category
pcoaPlot <- function(x){
subphylo <- prune_samples(sample_data(yass)$agecat==x, x=yass)
ordu = ordinate(subphylo, "PCoA", "unifrac", weighted=TRUE)
plot = plot_ordination(subphylo, ordu, color="abx", shape="human")
return(plot)
}agecats <- yass %>% sample_data() %>% .$agecat %>% unique()
listPlots <- lapply(agecats, pcoaPlot)Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.
Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"): Could not obtain coordinates from the provided `ordination`.
Please check your ordination method, and whether it is supported by `scores` or listed by phyloseq-package.Warning in plot_ordination(subphylo, ordu, color = "abx", shape = "human"):
Shape variable was not found in the available data you provided.No shape mapped.n <- length(listPlots)
nCol <- floor(sqrt(n))
do.call("grid.arrange", c(listPlots, ncol=nCol))