Lab 4
Diego Zapata
2024-05-01
STEP 1
#a
babybact <- read_table("lab4/backhed_data_stats202.txt")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double(),
## id = col_character(),
## age = col_character()
## )
## ℹ Use `spec()` for the full column specifications.#b
dim(babybact) #nrows = 391, # ncol = 697## [1] 391 697#c- Understand your data! What fields are present? What values are in each field?
Each row represents a microbiome sample from a baby of varying age. Columns represent the relative abundance of unique microbiota taxa. There is one column that denotes the baby’s age (newborn - 1 year old).
#d - Create a data frame d2 consisting of only the relative abundances of all the species (e.g., exclude columns 1 and 2). This will be used for the analyses below.
d2 <- babybact[,3:697]#e - Visualize your data! Make a histogram using hist() of the relative abundances of the common species Bacteroides_vulgatus_57955. Make sure to label your axes! Compare this distribution to that of any randomly chosen species. What do these distributions tell you?
#plotting Bacteroides vulgatus
par(mfrow = c(3,1))
hist(d2$Bacteroides_vulgatus_57955, main = "Relative abundance of Bacteroides vulgatus", xlab = "relative abundance", breaks = 20)
sum(d2$Bacteroides_vulgatus_57955 == 0) / (nrow(d2)) ## [1] 0.05626598sum(d2$Bacteroides_vulgatus_57955 > 0.2) / (nrow(d2)) ## [1] 0.1508951#Plotting two other random species
hist(d2$Salmonella_bongori_55351, main = "Relative abundance of Salmonella bongori", xlab = "relative abundance", breaks = 20)
hist(d2$Escherichia_coli_57907, main = "Relative abundance of Escherichia coli", xlab = "relative abundance", breaks = 20)
While the bacteria seems to be common across samples (ie only 5% of samples do not have this species), It mostly occurs at lower relative abundances. There are some samples however where there are much higher proportions of the bacteria. For B. vulgatus, ~15% of all samples have a relative abundance greater than 0.2! Other species have similar right-tailed distributions, but the ones I have selected have a lower and more narrow range of relative abundances compared to B. vulgatus.
#f - It turns out that the data is not normally distributed, and yet this is an assumption of PCA. We will transform it using a common statistical technique for microbiome data knows as the central log ratio (this is equivalent to scaling the data, but specific for compositional data in which features sum to one). To apply clr, type: d3 = clr(d2).
d3 <- clr(d2)STEP 2
#a - Perform PCA using PCA=prcomp() using the new dataframe.
PCA <- prcomp(d3)#b- Plot PC1 vs PC2 using autoplot(). You can plot and color the points by age as follows:
autoplot(PCA,
x=1,
y=2,
data=babybact,
colour = 'age')
#c - What does each point in the PCA plot represent? Which PCs separate which age groups?
Each point represents one gut microbiome from a baby (ie an individual sample). PC1 seems to separate samples by age group most, with newborn babies increasingly loading to the left of PC1 with age.
#d. Plot the variance explained by each PC using fviz_eig() using choice=”variance”. How much overall variance is explained by the first four PCs?
fviz_eig(PCA, choice = "variance")
sum <- summary(PCA)
sum$importance[2,]## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 0.12640 0.04658 0.03961 0.02915 0.02337 0.01943 0.01789 0.01655 0.01486 0.01417
## PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20
## 0.01352 0.01301 0.01168 0.01112 0.01067 0.01029 0.01013 0.00960 0.00906 0.00885
## PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30
## 0.00862 0.00826 0.00818 0.00809 0.00790 0.00774 0.00748 0.00717 0.00705 0.00691
## PC31 PC32 PC33 PC34 PC35 PC36 PC37 PC38 PC39 PC40
## 0.00664 0.00646 0.00635 0.00627 0.00607 0.00599 0.00589 0.00564 0.00555 0.00544
## PC41 PC42 PC43 PC44 PC45 PC46 PC47 PC48 PC49 PC50
## 0.00538 0.00534 0.00522 0.00512 0.00507 0.00504 0.00493 0.00481 0.00472 0.00467
## PC51 PC52 PC53 PC54 PC55 PC56 PC57 PC58 PC59 PC60
## 0.00459 0.00447 0.00442 0.00431 0.00424 0.00417 0.00412 0.00407 0.00397 0.00395
## PC61 PC62 PC63 PC64 PC65 PC66 PC67 PC68 PC69 PC70
## 0.00390 0.00382 0.00380 0.00374 0.00370 0.00364 0.00360 0.00353 0.00350 0.00344
## PC71 PC72 PC73 PC74 PC75 PC76 PC77 PC78 PC79 PC80
## 0.00336 0.00334 0.00333 0.00326 0.00323 0.00319 0.00316 0.00312 0.00308 0.00302
## PC81 PC82 PC83 PC84 PC85 PC86 PC87 PC88 PC89 PC90
## 0.00301 0.00294 0.00291 0.00290 0.00285 0.00284 0.00280 0.00276 0.00270 0.00264
## PC91 PC92 PC93 PC94 PC95 PC96 PC97 PC98 PC99 PC100
## 0.00262 0.00257 0.00254 0.00251 0.00249 0.00248 0.00244 0.00237 0.00235 0.00234
## PC101 PC102 PC103 PC104 PC105 PC106 PC107 PC108 PC109 PC110
## 0.00234 0.00232 0.00228 0.00225 0.00222 0.00220 0.00215 0.00214 0.00213 0.00211
## PC111 PC112 PC113 PC114 PC115 PC116 PC117 PC118 PC119 PC120
## 0.00208 0.00204 0.00203 0.00202 0.00199 0.00196 0.00194 0.00191 0.00189 0.00188
## PC121 PC122 PC123 PC124 PC125 PC126 PC127 PC128 PC129 PC130
## 0.00186 0.00183 0.00182 0.00182 0.00179 0.00177 0.00173 0.00172 0.00171 0.00169
## PC131 PC132 PC133 PC134 PC135 PC136 PC137 PC138 PC139 PC140
## 0.00166 0.00164 0.00162 0.00159 0.00157 0.00156 0.00153 0.00151 0.00150 0.00149
## PC141 PC142 PC143 PC144 PC145 PC146 PC147 PC148 PC149 PC150
## 0.00148 0.00146 0.00146 0.00145 0.00144 0.00141 0.00140 0.00139 0.00136 0.00135
## PC151 PC152 PC153 PC154 PC155 PC156 PC157 PC158 PC159 PC160
## 0.00135 0.00132 0.00131 0.00130 0.00128 0.00125 0.00124 0.00123 0.00122 0.00121
## PC161 PC162 PC163 PC164 PC165 PC166 PC167 PC168 PC169 PC170
## 0.00118 0.00118 0.00116 0.00114 0.00114 0.00112 0.00111 0.00109 0.00107 0.00107
## PC171 PC172 PC173 PC174 PC175 PC176 PC177 PC178 PC179 PC180
## 0.00106 0.00104 0.00103 0.00102 0.00101 0.00099 0.00098 0.00097 0.00095 0.00095
## PC181 PC182 PC183 PC184 PC185 PC186 PC187 PC188 PC189 PC190
## 0.00094 0.00093 0.00092 0.00091 0.00091 0.00090 0.00088 0.00088 0.00087 0.00085
## PC191 PC192 PC193 PC194 PC195 PC196 PC197 PC198 PC199 PC200
## 0.00084 0.00083 0.00082 0.00082 0.00081 0.00080 0.00079 0.00078 0.00077 0.00076
## PC201 PC202 PC203 PC204 PC205 PC206 PC207 PC208 PC209 PC210
## 0.00076 0.00075 0.00074 0.00073 0.00072 0.00071 0.00071 0.00069 0.00069 0.00067
## PC211 PC212 PC213 PC214 PC215 PC216 PC217 PC218 PC219 PC220
## 0.00067 0.00066 0.00066 0.00065 0.00064 0.00063 0.00062 0.00062 0.00061 0.00061
## PC221 PC222 PC223 PC224 PC225 PC226 PC227 PC228 PC229 PC230
## 0.00061 0.00060 0.00059 0.00058 0.00057 0.00056 0.00056 0.00055 0.00054 0.00054
## PC231 PC232 PC233 PC234 PC235 PC236 PC237 PC238 PC239 PC240
## 0.00053 0.00052 0.00052 0.00052 0.00051 0.00050 0.00049 0.00049 0.00049 0.00047
## PC241 PC242 PC243 PC244 PC245 PC246 PC247 PC248 PC249 PC250
## 0.00047 0.00047 0.00046 0.00045 0.00045 0.00043 0.00043 0.00043 0.00042 0.00041
## PC251 PC252 PC253 PC254 PC255 PC256 PC257 PC258 PC259 PC260
## 0.00041 0.00040 0.00040 0.00039 0.00039 0.00038 0.00038 0.00037 0.00037 0.00036
## PC261 PC262 PC263 PC264 PC265 PC266 PC267 PC268 PC269 PC270
## 0.00036 0.00035 0.00035 0.00035 0.00034 0.00034 0.00033 0.00033 0.00032 0.00032
## PC271 PC272 PC273 PC274 PC275 PC276 PC277 PC278 PC279 PC280
## 0.00031 0.00031 0.00031 0.00030 0.00029 0.00029 0.00029 0.00028 0.00028 0.00028
## PC281 PC282 PC283 PC284 PC285 PC286 PC287 PC288 PC289 PC290
## 0.00027 0.00026 0.00026 0.00025 0.00025 0.00025 0.00025 0.00024 0.00023 0.00023
## PC291 PC292 PC293 PC294 PC295 PC296 PC297 PC298 PC299 PC300
## 0.00023 0.00022 0.00022 0.00021 0.00021 0.00021 0.00021 0.00020 0.00020 0.00020
## PC301 PC302 PC303 PC304 PC305 PC306 PC307 PC308 PC309 PC310
## 0.00020 0.00019 0.00019 0.00019 0.00018 0.00018 0.00017 0.00017 0.00017 0.00017
## PC311 PC312 PC313 PC314 PC315 PC316 PC317 PC318 PC319 PC320
## 0.00016 0.00016 0.00016 0.00015 0.00015 0.00015 0.00015 0.00014 0.00014 0.00014
## PC321 PC322 PC323 PC324 PC325 PC326 PC327 PC328 PC329 PC330
## 0.00014 0.00013 0.00013 0.00013 0.00013 0.00012 0.00012 0.00012 0.00011 0.00011
## PC331 PC332 PC333 PC334 PC335 PC336 PC337 PC338 PC339 PC340
## 0.00011 0.00011 0.00011 0.00011 0.00010 0.00010 0.00010 0.00010 0.00009 0.00009
## PC341 PC342 PC343 PC344 PC345 PC346 PC347 PC348 PC349 PC350
## 0.00009 0.00009 0.00009 0.00009 0.00008 0.00008 0.00008 0.00008 0.00007 0.00007
## PC351 PC352 PC353 PC354 PC355 PC356 PC357 PC358 PC359 PC360
## 0.00007 0.00007 0.00007 0.00007 0.00006 0.00006 0.00006 0.00006 0.00006 0.00006
## PC361 PC362 PC363 PC364 PC365 PC366 PC367 PC368 PC369 PC370
## 0.00005 0.00005 0.00005 0.00005 0.00005 0.00004 0.00004 0.00004 0.00004 0.00004
## PC371 PC372 PC373 PC374 PC375 PC376 PC377 PC378 PC379 PC380
## 0.00004 0.00004 0.00003 0.00003 0.00003 0.00003 0.00003 0.00003 0.00003 0.00003
## PC381 PC382 PC383 PC384 PC385 PC386 PC387 PC388 PC389 PC390
## 0.00003 0.00002 0.00002 0.00002 0.00002 0.00002 0.00002 0.00001 0.00001 0.00001
## PC391
## 0.00000var_ex <- sum$importance[2,]
PC1to4 <- sum(var_ex[1:4])
PC1to4 #24.174% of variance is explained by the first four PCs.## [1] 0.24174#e. Plot the loadings for the PCA using fviz_pca_var(). Which species are the
#most informative? Note, you will use this information below in the k-means exercise.
fviz_pca_var(PCA)
It looks like the Escherichia spp./Bacteroides spp. and Ruminococcus_gnavis seem to load very much on PC1 and PC2 respectively.
###Part 2
#1. Perform K-means clustering on the data using kmeans(d2, K). Which K will you use? Note, we are not using the clr transformed data here because R is a bit buggy.
I will use k=3 since it looked like there were three major clusters according to the PCA, but for comparison, I will include K=2 and 4.
k2 <- kmeans(d2,2)
k3 <- kmeans(d2,3)
k4 <- kmeans(d2,4)#2. How well did the different ages get classified? Try making a table to assess this using table().
It seems the age groups did not get classified very well according to their microbiomes. There are equal numbers of each age class, but each k-means clustering failed to reproduce these evenly distributed clusters. In K=4, there was a large skew towards the first cluster. The same unevenness applies to k=2 and 3.
table(babybact$age)##
## 12M 4M B M
## 97 98 98 98table(k2$cluster)##
## 1 2
## 35 356table(k3$cluster)##
## 1 2 3
## 328 31 32table(k4$cluster)##
## 1 2 3 4
## 45 230 83 33#3 Plot the K-means clustering assignments using two species that seemed particularly informative from question 2e above. Color the points using the results of the K-means clustering (remember, you can pull out the cluster that each point is in using $cluster). Repeat for another pair of species.
# Extract unique cluster labels
unique_clusters <- unique(k2$cluster)
s=sort(unique(k2$cluster))
c=1:length(s)
names(c)=s
cs=c[babybact$age]
d2$cluster <- k2$cluster
plot(d2$Escherichia_fergusonii_56914,
d2$Bacteroides_ovatus_58035,
xlab = "Escherichia albertii_56276",
ylab = "Ruminococcus gnavus_57638",
col = d2$cluster)
legend("topright", legend = unique_clusters, col = unique_clusters, pch = 1, title = "Cluster")
unique_clusters <- unique(k2$cluster)
plot(d2$Barnesiella_intestinihominis_62208,
d2$Escherichia_albertii_56276,
xlab = "Barnesiella intestinihominis_62208",
ylab = "Escherichia albertii_56276",
col = k2$cluster)
legend("topright", legend = unique_clusters, col = unique_clusters, pch = 1, title = "Cluster")
#Step 4 - Now, make the same plots again as you just did, except color the points according to the actual ages specified in the original data. How well does the clustering perform in accurately labeling the points? Note, you might need to assign a color to each age. Here is some code to help you:
#creating colors for each age
strings=sort(unique(babybact$age))
colors=1:length(strings)
names(colors)=strings
colors=colors[babybact$age]
#plotting by age group instead
plot(d2$Escherichia_fergusonii_56914,
d2$Bacteroides_ovatus_58035,
xlab = "Escherichia albertii_56276",
ylab = "Ruminococcus gnavus_57638",
col = colors)
legend("topright", legend = strings, col = unique(colors), pch = 1, title = "Age Group")
plot(d2$Barnesiella_intestinihominis_62208,
d2$Escherichia_albertii_56276,
xlab = "Barnesiella intestinihominis_62208",
ylab = "Escherichia albertii_56276",
col = colors)
legend("topright", legend = strings, col = unique(colors), pch = 1, title = "Age Group")
Actually using the real age groups in data show a much more clear bifurcation between age group ~ microbiome abundances. For the first plot (Escherichia_fergusonii_56914 vs Bacteroides_ovatus_58035), there is a clearly more R. gnavus in newborn babies than in other age groups. E. albertii is largely more abundant in babies other than newborns. For the seconf plot (Barnesiella intestinihominis_62208 vs Escherichia albertii_56276), B. intestinihominis is most common in 4M old babies. Other age classes largely do not have this bacteria in significant amounts compared to 4m babies.