.Rhistory

#extract samples of interest with specified strings - returns vector of strings and NAs
extract.lungsamp.sampctrls <- str_extract(rownames(otu_good), pattern="_BAL_|_Lung_|PBS|CLEAN|USED|Homog")
#select subsets based on extracted samp names - returns a logical vector (NA --> FALSE)
extract.lungsamp.sampctrls.good <- !is.na(extract.lungsamp.sampctrls)
#subset otu counts
otu.good.lungsamp.sampctrls <- otu_good[extract.lungsamp.sampctrls.good, ]
#subset otu.df for PCA
otu.df.lungsamp.sampctrls <- filter(otu_df, Sample_Type == "BAL" | Sample_Type == "Lung" | Sample_Type == "PBS" | Sample_Type == "Clean_Rinse" | Sample_Type == "Used_Rinse" | Sample_Type == "Homog_Ctrl")
#make a hellinger normalized df and pca object from the normalized df
make_hel_pca(otu.good.lungsamp.sampctrls)
# summarize pca object to allow for pulling the coordinates in the PC space
lungsamp_sampctrls_pca_summary <- summary(otu.good.lungsamp.sampctrls_pca)
# select PC1 and PC2 to plot
#"sites" refers to each individual sample, i.e. ecologic sites
pc_space <- lungsamp_sampctrls_pca_summary$sites %>% as.data.frame() %>% rownames_to_column() %>% dplyr::select(1:3)
colnames(pc_space) <- c("Sample_Name", "PC1", "PC2")
pc_loadings <- lungsamp_sampctrls_pca_summary$species # species as loading vectors
# attach your data from the pca object to your labels
centroid_df <- data.frame(Sample_Name = pc_space$Sample_Name, Sample_Type2 = otu.df.lungsamp.sampctrls[,"Sample_Type2"], PC1 = pc_space$PC1, PC2 = pc_space$PC2)
# calculate the mean x coordinate and y coordinate for each sample in PC1/PC2 space
centroids <- centroid_df %>%
group_by(Sample_Type2) %>%
summarize(x_coord=mean(PC1), y_coord=mean(PC2)) %>%
ungroup()
# to EACH SUBJECT in the pc object, add a column which contains the x and y
#coordinates for the centroid for the group that the subject belongs to
pc_plot_data <- full_join(centroid_df, centroids, by ="Sample_Type2")
# Now you can use your dataframe and mainpulate it how you want
pc_plot_data$Sample_Type2 <- factor(pc_plot_data$Sample_Type2,
levels = c("Lung", "BAL", "PBS", "Syringe_Rinse",  "Homog_Ctrl"),
labels = c("Lung", "BAL", "Saline", "Syringe Rinse",  "Homogenization\nControl"))
pca_pal <- c( "#EB3D0E", #scarlet - lung
"#0300A6", # blue - BAL
"#666B6E", #- saline
"#9A99E9", #- used rinse
"#F2C121" # saffron - water
)
# plot PCA
ggplot(pc_plot_data, aes(x = PC1, y = PC2)) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_segment(aes(x=x_coord, xend = PC1, y=y_coord, yend = PC2), color = "black") +
geom_point(show.legend = FALSE, aes(color = Sample_Type2), size = 3) +
# segment starting at centroid, ending at point
geom_label(aes(x=x_coord, y=y_coord, label=Sample_Type2, color=Sample_Type2), size=3, show.legend = FALSE) + # centroid information
theme_classic() +
theme(legend.title = element_blank(),
axis.text.y = element_text(size=15),
axis.text.x = element_text(size=15),
axis.title.y = element_text(size=12, face="bold", vjust = 0),
axis.title.x = element_text(size=12, face="bold", vjust = 0),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")
) +
scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, 0.5)) +
scale_y_continuous(limits = c(-1.1, 1.1), breaks = seq(-1,1,0.5)) +
labs(
x = "\nPrincipal Component 1 (13.7% explained)",
y = "Principal Component 2  (7.0% explained)\n"
) + scale_color_manual(values=pca_pal)
#prepare hellinger-transformed df to be filtered for pairwise comparisons
adonis.select.vec <- dplyr::select(otu.df.lungsamp.sampctrls, Sample:Sample_Type)
adonis.hel.df <- rownames_to_column(otu.good.lungsamp.sampctrls_hel, var = "Sample")
adonis.hel.df <- left_join(adonis.hel.df, adonis.select.vec, by = "Sample") %>%
dplyr::select(Sample, Sample_Type, everything())
#multiple comparisons - whole lung v. BAL v. pooled sampling controls
adonis.hel.df.wbn <- column_to_rownames(adonis.hel.df, var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis(adonis.hel.df.wbn~otu.df.lungsamp.sampctrls$RA_Groups, method="euclidean", permutations=10000)
#pairwise comparisons - whole lung v. bal
adonis.hel.bl <- dplyr::filter(adonis.hel.df, Sample_Type == "Lung" | Sample_Type == "BAL") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.bl <- dplyr::filter(otu.df.lungsamp.sampctrls,
Sample_Type == "Lung" | Sample_Type == "BAL")
adonis(adonis.hel.bl~adonis.otudf.bl$Sample_Type, method="euclidean", permutations=10000)
#pairwise comparison - whole lung v. sampling negative controls
adonis.hel.wn <- dplyr::filter(adonis.hel.df, Sample_Type != "BAL") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.wn <- dplyr::filter(otu.df.lungsamp.sampctrls, Sample_Type != "BAL")
adonis(adonis.hel.wn~adonis.otudf.wn$RA_Groups, method="euclidean", permutations=10000)
#pairwise comparison - BAL v. sampling negative controls
adonis.hel.bn <- dplyr::filter(adonis.hel.df, Sample_Type != "Lung") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.bn <- dplyr::filter(otu.df.lungsamp.sampctrls, Sample_Type != "Lung")
adonis(adonis.hel.bn~adonis.otudf.bn$Organ, method="euclidean", permutations=10000)
#extract samples of interest with specified strings - returns vector of strings and NAs
extract.lungsamp.isoctrls <- str_extract(rownames(otu_good), pattern="_BAL_|_Lung_|_IsoCtrl_|AE")
#select subsets based on extracted samp names - returns a logical vector (NA --> FALSE)
extract.lungsamp.isoctrls.good <- !is.na(extract.lungsamp.isoctrls)
#subset otu counts
otu.good.lungsamp.isoctrls <- otu_good[extract.lungsamp.isoctrls.good, ]
#subset otu.df for PCA
otu.df.lungsamp.isoctrls <- filter(otu_df, Sample_Type == "BAL" | Sample_Type == "Lung" | Sample_Type == "Iso_Ctrl" | Sample_Type == "AE")
#make a hellinger normalized df and pca object from the normalized df
make_hel_pca(otu.good.lungsamp.isoctrls)
# summarize pca object to allow for pulling the coordinates in the PC space
lungsamp_isoctrls_pca_summary <- summary(otu.good.lungsamp.isoctrls_pca)
# select PC1 and PC2 to plot
#"sites" refers to each individual sample, i.e. ecologic sites
pc_space <- lungsamp_isoctrls_pca_summary$sites %>% as.data.frame() %>% rownames_to_column() %>% dplyr::select(1:3)
colnames(pc_space) <- c("Sample_Name", "PC1", "PC2")
pc_loadings <- lungsamp_isoctrls_pca_summary$species  # species as loading vectors
# attach your data from the pca object to your labels
centroid_df <- data.frame(Sample_Name = pc_space$Sample_Name, Sample_Type2 = otu.df.lungsamp.isoctrls[,"Sample_Type2"], PC1 = pc_space$PC1, PC2 = pc_space$PC2)
# calculate the mean x coordinate and y coordinate for each sample in PC1/PC2 space
centroids <- centroid_df %>%
group_by(Sample_Type2) %>%
summarize(x_coord=mean(PC1), y_coord=mean(PC2)) %>%
ungroup()
# to EACH SUBJECT in the pc object, add a column which contains the x and y
#coordinates for the centroid for the group that the subject belongs to
pc_plot_data <- full_join(centroid_df, centroids, by ="Sample_Type2")
# Now you can use your dataframe and mainpulate it how you want
pc_plot_data$Sample_Type2 <- factor(pc_plot_data$Sample_Type2,
levels = c("Lung", "BAL", "AE", "Iso_Ctrl"),
labels = c("Lung", "BAL", "Elution\nBuffer", "Isolation\nControl"))
pca_pal <- c( "#EB3D0E", #scarlet - lung
"#0300A6", # blue - BAL
"#DBE6EC", #- AE
"#666B6E" #- iso ctrl
)
# plot PCA
ggplot(pc_plot_data, aes(x = PC1, y = PC2)) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_segment(aes(x=x_coord, xend = PC1, y=y_coord, yend = PC2), color = "black") +
geom_point(show.legend = FALSE, aes(color = Sample_Type2), size = 3) +
# segment starting at centroid, ending at point
geom_label(aes(x=x_coord, y=y_coord, label=Sample_Type2, color=Sample_Type2), size=3, show.legend = FALSE) + # centroid information
theme_classic() +
theme(legend.title = element_blank(),
axis.text.y = element_text(size=15),
axis.text.x = element_text(size=15),
axis.title.y = element_text(size=12, face="bold", vjust = 0),
axis.title.x = element_text(size=12, face="bold", vjust = 0),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")
) +
scale_color_manual(values = pca_pal) +
scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, 0.5)) +
scale_y_continuous(limits = c(-1.1, 1.1), breaks = seq(-1,1,0.5)) +
labs(
x = "\nPrincipal Component 1 (12.1% explained)",
y = "Principal Component 2  (9.0% explained)\n"
)
#prepare hellinger-transformed df to be filtered for pairwise comparisons
adonis.select.vec <- dplyr::select(otu.df.lungsamp.isoctrls, Sample:Sample_Type)
adonis.hel.df <- rownames_to_column(otu.good.lungsamp.isoctrls_hel, var = "Sample")
adonis.hel.df <- left_join(adonis.hel.df, adonis.select.vec, by = "Sample") %>%
dplyr::select(Sample, Sample_Type, everything())
#multiple comparisons - whole lung v. BAL v. pooled isolation controls
adonis(otu.good.lungsamp.isoctrls_hel~otu.df.lungsamp.isoctrls$Sample_Type,
method="euclidean", permutations=10000)
#pairwise comparison - whole lung v. iso ctrls
adonis.hel.il <- dplyr::filter(adonis.hel.df, Sample_Type != "BAL") %>%
column_to_rownames(var = "Sample") %>% dplyr::select(-Sample_Type)
adonis.otudf.il <- dplyr::filter(otu.df.lungsamp.isoctrls, Sample_Type != "BAL")
adonis(adonis.hel.il~adonis.otudf.il$RA_Groups, method="euclidean", permutations=10000)
#pairwise comparison - BAL v. iso ctrls
adonis.hel.ib <- dplyr::filter(adonis.hel.df, Sample_Type != "Lung") %>%
column_to_rownames(var = "Sample") %>% dplyr::select(-Sample_Type)
adonis.otudf.ib <- dplyr::filter(otu.df.lungsamp.isoctrls, Sample_Type != "Lung")
adonis(adonis.hel.ib~adonis.otudf.ib$RA_Groups, method="euclidean", permutations=10000)
#extract samples of interest with specified strings - returns vector of strings and NAs
extract.lungsamp.seqctrls <- str_extract(rownames(otu_good), pattern="_BAL_|_Lung_|Water|Blank")
#select subsets based on extracted samp names - returns a logical vector (NA --> FALSE)
extract.lungsamp.seqctrls.good <- !is.na(extract.lungsamp.seqctrls)
#subset otu counts
otu.good.lungsamp.seqctrls <- otu_good[extract.lungsamp.seqctrls.good, ]
#subset otu.df for PCA
otu.df.lungsamp.seqctrls <- filter(otu_df, Sample_Type == "BAL" | Sample_Type == "Lung" | Sample_Type == "Water" | Sample_Type == "Empty")
#make a hellinger normalized df and pca object from the normalized df
make_hel_pca(otu.good.lungsamp.seqctrls)
# summarize pca object to allow for pulling the coordinates in the PC space
lungsamp_seqctrls_pca_summary <- summary(otu.good.lungsamp.seqctrls_pca)
# select PC1 and PC2 to plot
#"sites" refers to each individual sample, i.e. ecologic sites
pc_space <- lungsamp_seqctrls_pca_summary$sites %>% as.data.frame() %>% rownames_to_column() %>% dplyr::select(1:3)
colnames(pc_space) <- c("Sample_Name", "PC1", "PC2")
pc_loadings <- lungsamp_seqctrls_pca_summary$species # species as loading vectors
# attach your data from the pca object to your labels
centroid_df <- data.frame(Sample_Name = pc_space$Sample_Name, Sample_Type2 = otu.df.lungsamp.seqctrls[,"Sample_Type2"], PC1 = pc_space$PC1, PC2 = pc_space$PC2)
# calculate the mean x coordinate and y coordinate for each sample in PC1/PC2 space
centroids <- centroid_df %>%
group_by(Sample_Type2) %>%
summarize(x_coord=mean(PC1), y_coord=mean(PC2)) %>%
ungroup()
# to EACH SUBJECT in the pc object, add a column which contains the x and y
#coordinates for the centroid for the group that the subject belongs to
pc_plot_data <- full_join(centroid_df, centroids, by ="Sample_Type2")
# Now you can use your dataframe and mainpulate it how you want
pc_plot_data$Sample_Type2 <- factor(pc_plot_data$Sample_Type2,
levels = c("Lung", "BAL", "Empty", "Water"),
labels = c("Lung", "BAL", "Empty\nWell", "Water"))
pca_pal <- c( "#EB3D0E", #scarlet - lung
"#0300A6", # blue - BAL
"#F2C121", #- Empty
"#666B6E" #- Water
)
ggplot(pc_plot_data, aes(x = PC1, y = PC2)) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_segment(aes(x=x_coord, xend = PC1, y=y_coord, yend = PC2), color = "black") +
geom_point(show.legend = FALSE, aes(color = Sample_Type2), size = 3) +
# segment starting at centroid, ending at point
geom_label(aes(x=x_coord, y=y_coord, label=Sample_Type2, color=Sample_Type2), size=3, show.legend = FALSE) + # centroid information
theme_classic() +
theme(legend.title = element_blank(),
axis.text.y = element_text(size=15),
axis.text.x = element_text(size=15),
axis.title.y = element_text(size=12, face="bold", vjust = 0),
axis.title.x = element_text(size=12, face="bold", vjust = 0),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")
) +
scale_color_manual(values = pca_pal) +
scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, 0.5)) +
scale_y_continuous(limits = c(-1.1, 1.1), breaks = seq(-1,1,0.5)) +
labs(
x = "\nPrincipal Component 1 (13.5% explained)",
y = "Principal Component 2  (6.6% explained)\n"
)
#prepare hellinger-transformed df to be filtered for pairwise comparisons
adonis.select.vec <- dplyr::select(otu.df.lungsamp.seqctrls, Sample:Sample_Type)
adonis.hel.df <- rownames_to_column(otu.good.lungsamp.seqctrls_hel, var = "Sample")
adonis.hel.df <- left_join(adonis.hel.df, adonis.select.vec, by = "Sample") %>%
dplyr::select(Sample, Sample_Type, everything())
#multiple comparisons - whole lung v. BAL v. pooled sequencing controls
adonis(otu.good.lungsamp.seqctrls_hel~otu.df.lungsamp.seqctrls$Sample_Type,
method="euclidean",
permutations=10000)
#pairwise comparison - whole lung v. seq ctrls
adonis.hel.sl <- dplyr::filter(adonis.hel.df, Sample_Type != "BAL") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.sl <- dplyr::filter(otu.df.lungsamp.seqctrls, Sample_Type != "BAL")
adonis(adonis.hel.sl~adonis.otudf.sl$RA_Groups, method="euclidean", permutations=10000)
#pairwise comparison - BAL v. seq ctrls
adonis.hel.sb <- dplyr::filter(adonis.hel.df, Sample_Type != "Lung") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.sb <- dplyr::filter(otu.df.lungsamp.seqctrls, Sample_Type != "Lung")
adonis(adonis.hel.sb~adonis.otudf.sb$RA_Groups, method="euclidean", permutations=10000)
#extract samples of interest with specified strings - returns vector of strings and NAs
extract.tong.lungsamp <- str_extract(rownames(otu_good), pattern="Tong|_BAL_|_Lung_")
#select subsets based on extracted samp names - returns a logical vector (NA --> FALSE)
extract.tong.lungsamp.good <- !is.na(extract.tong.lungsamp)
#subset otu counts
otu.good.tong.lungsamp <- otu_good[extract.tong.lungsamp.good, ]
#subset otu.df for PCA
otu.df.tong.lungsamp <- filter(otu_df, Sample_Type == "Tongue" | Sample_Type == "BAL" | Sample_Type == "Lung")
otu.df.tong.lungsamp <- mutate(otu.df.tong.lungsamp, PCA_Groups = ifelse(
test = Organ == "Lung",
yes =  as.character(RA_Groups),
no = as.character(Organ))) %>%
dplyr::select(Sample:Gene_16S_copies_per_mL, PCA_Groups, everything())
otu.df.tong.lungsamp$PCA_Groups <- as.factor(otu.df.tong.lungsamp$PCA_Groups)
#make a hellinger normalized df and pca object from the normalized df
make_hel_pca(otu.good.tong.lungsamp)
tong_lungsamp_pca_summary <- summary(otu.good.tong.lungsamp_pca) # summarize pca object to allow for pulling the coordinates in the PC space
# select PC1 and PC2 for plotting
#"sites" refers to each individual sample, i.e. ecologic sites
pc_space <- tong_lungsamp_pca_summary$sites %>% as.data.frame() %>% rownames_to_column() %>% dplyr::select(1:3)
colnames(pc_space) <- c("Sample_Name", "PC1", "PC2")
pc_loadings <- tong_lungsamp_pca_summary$species # species as loading vectors
# attach your data from the pca object to your labels
centroid_df <- data.frame(Sample_Name = pc_space$Sample_Name, PCA_Groups = otu.df.tong.lungsamp[,"PCA_Groups"], PC1 = pc_space$PC1, PC2 = pc_space$PC2)
# calculate the mean x coordinate and y coordinate for each sample in PC1/PC2 space
centroids <- centroid_df %>%
group_by(PCA_Groups) %>%
summarize(x_coord=mean(PC1), y_coord=mean(PC2)) %>%
ungroup()
# to EACH SUBJECT in the pc object, add a column which contains the x and y
#coordinates for the centroid for the group that the subject belongs to
pc_plot_data <- full_join(centroid_df, centroids, by ="PCA_Groups")
# Now you can use your dataframe and mainpulate it how you want
pc_plot_data$PCA_Groups <- factor(pc_plot_data$PCA_Groups,
levels = c("Lung", "BAL", "Tongue"))
pca_pal <- c( "#EB3D0E", #scarlet - lung
"#0300A6", # blue - BAL
"#FB9986" #vivid tangerine - tongue
)
ggplot(pc_plot_data, aes(x = PC1, y = PC2)) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_segment(aes(x=x_coord, xend = PC1, y=y_coord, yend = PC2), color = "black") +
geom_point(show.legend = FALSE, aes(color = PCA_Groups), size = 2) +
# segment starting at centroid, ending at point
geom_label(aes(x=x_coord, y=y_coord, label=PCA_Groups, color=PCA_Groups), size=3, show.legend = FALSE) + # centroid information
theme_classic() +
theme(legend.title = element_blank(),
axis.text.y = element_text(size=15),
axis.text.x = element_text(size=15),
axis.title.y = element_text(size=12, face="bold", vjust = 0),
axis.title.x = element_text(size=12, face="bold", vjust = 0),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")
) +
scale_color_manual(values = pca_pal) +
scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, 0.5)) +
scale_y_continuous(limits = c(-1.1, 1.1), breaks = seq(-1,1,0.5)) +
labs(
x = "\nPrincipal Component 1 (10.5% explained)",
y = "Principal Component 2  (4.6% explained)\n"
)
#multiple comparisons - whole lung v. BAL v. tongue
adonis(otu.good.tong.lungsamp_hel~otu.df.tong.lungsamp$Sample_Type,
method="euclidean",
permutations=10000)
#prepare hellinger-transformed df to be filtered for pairwise comparisons
adonis.select.vec <- dplyr::select(otu.df.tong.lungsamp, Sample:Sample_Type)
adonis.hel.df <- rownames_to_column(otu.good.tong.lungsamp_hel, var = "Sample")
adonis.hel.df <- left_join(adonis.hel.df, adonis.select.vec, by = "Sample") %>%
dplyr::select(Sample, Sample_Type, everything())
#pairwise comparison - whole lung v. tongue
adonis.hel.tl <- dplyr::filter(adonis.hel.df, Sample_Type == "Lung" | Sample_Type == "Tongue") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.tl <- dplyr::filter(otu.df.tong.lungsamp,
Sample_Type == "Lung" | Sample_Type == "Tongue")
adonis(adonis.hel.tl~adonis.otudf.tl$Sample_Type, method="euclidean", permutations=10000)
#pairwise comparison - BAL v. tongue
adonis.hel.tb <- dplyr::filter(adonis.hel.df, Sample_Type == "BAL" | Sample_Type == "Tongue") %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(-Sample_Type)
adonis.otudf.tb <- dplyr::filter(otu.df.tong.lungsamp, Sample_Type == "BAL" | Sample_Type == "Tongue")
adonis(adonis.hel.tb~adonis.otudf.tb$Sample_Type, method="euclidean", permutations=10000)
#extract samples of interest with specified strings - returns vector of strings and NAs
extract.tong.lungsamp.allnegctrls <- str_extract(rownames(otu_good), pattern="Tong|_BAL_|_Lung_|PBS|CLEAN|USED|Homog|AE|Iso|Water|Blank")
#select subsets based on extracted samp names - returns a logical vector (NA --> FALSE)
extract.tong.lungsamp.allnegctrls.good <- !is.na(extract.tong.lungsamp.allnegctrls)
#subset otu counts
otu.good.tong.lungsamp.allnegctrls <- otu_good[extract.tong.lungsamp.allnegctrls.good, ]
#subset otu.df for PCA
otu.df.tong.lungsamp.allnegctrls <- filter(otu_df, Sample_Type == "Tongue" | Sample_Type == "BAL" | Sample_Type == "Lung" | Sample_Type == "PBS" | Sample_Type == "Clean_Rinse" | Sample_Type == "Used_Rinse" | Sample_Type == "Homog_Ctrl" | Sample_Type == "AE" | Sample_Type == "Iso_Ctrl" | Sample_Type == "Water" | Sample_Type == "Empty")
otu.df.tong.lungsamp.allnegctrls <- mutate(otu.df.tong.lungsamp.allnegctrls, PCA_Groups = ifelse(
test = Organ == "Lung",
yes =  as.character(RA_Groups),
no = as.character(Organ))) %>%
dplyr::select(Sample:Gene_16S_copies_per_mL, PCA_Groups, everything())
otu.df.tong.lungsamp.allnegctrls$PCA_Groups <- as.factor(otu.df.tong.lungsamp.allnegctrls$PCA_Groups)
#make a hellinger normalized df and pca object from the normalized df
make_hel_pca(otu.good.tong.lungsamp.allnegctrls)
# summarize pca object to allow for pulling the coordinates in the PC space
tong_lungsamp_relevnegctrls_pca_summary <- summary(otu.good.tong.lungsamp.allnegctrls_pca)
# select PC1 and PC2 for plotting
#"sites" refers to each individual sample, i.e. ecologic sites
pc_space <- tong_lungsamp_relevnegctrls_pca_summary$sites %>% as.data.frame() %>% rownames_to_column() %>% dplyr::select(1:3)
colnames(pc_space) <- c("Sample_Name", "PC1", "PC2")
pc_loadings <- tong_lungsamp_relevnegctrls_pca_summary$species # species as loading vectors
# attach your data from the pca object to your labels
centroid_df <- data.frame(Sample_Name = pc_space$Sample_Name, PCA_Groups = otu.df.tong.lungsamp.allnegctrls[,"PCA_Groups"], PC1 = pc_space$PC1, PC2 = pc_space$PC2)
# calculate the mean x coordinate and y coordinate for each sample in PC1/PC2 space
centroids <- centroid_df %>%
group_by(PCA_Groups) %>%
summarize(x_coord=mean(PC1), y_coord=mean(PC2)) %>%
ungroup()
# to EACH SUBJECT in the pc object, add a column which contains the x and y
#coordinates for the centroid for the group that the subject belongs to
pc_plot_data <- full_join(centroid_df, centroids, by ="PCA_Groups")
# Now you can use your dataframe and mainpulate it how you want
pc_plot_data$PCA_Groups <- factor(pc_plot_data$PCA_Groups,
levels = c("Lung", "BAL", "Tongue", "Sampling Control", "Isolation Control", "Sequencing Control"))
pca_pal <- c( "#EB3D0E", #scarlet - lung
"#0300A6", # blue - BAL
"#FB9986", #vivid tangerine - tongue
"#0BDBA7", #sea foam green - sampling ctrls
"#F2C121", #- isolation ctrls
"#666B6E" #- sequencing ctrls
)
# plot PCA
ggplot(pc_plot_data, aes(x = PC1, y = PC2)) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_segment(aes(x=x_coord, xend = PC1, y=y_coord, yend = PC2), color = "black") +
geom_point(show.legend = FALSE, aes(color = PCA_Groups), size = 2) +
# segment starting at centroid, ending at point
geom_label(aes(x=x_coord, y=y_coord, label=PCA_Groups, color=PCA_Groups), size=3, show.legend = FALSE) + # centroid information
theme_classic() +
theme(legend.title = element_blank(),
axis.text.y = element_text(size=15),
axis.text.x = element_text(size=15),
axis.title.y = element_text(size=12, face="bold", vjust = 0),
axis.title.x = element_text(size=12, face="bold", vjust = 0),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")
) +
scale_color_manual(values = pca_pal) +
scale_x_continuous(limits = c(-0.75, 0.75), breaks =  c(-0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75)) +
scale_y_continuous(limits = c(-0.75, 0.75), breaks = c(-0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75)) +
labs(
x = "\nPrincipal Component 1 (11.6% explained)",
y = "Principal Component 2  (6.3% explained)\n"
)
#Rarefaction curves can also be used to explore α-diversity. First, the data needs to be prepared.
otu_raw_rare <- otu_raw
rownames(otu_raw_rare) <-  str_remove(string = rownames(otu_raw_rare), pattern = "_S\\d+_L001_R1_001") %>% str_remove(pattern = "_S\\d+_L001_R1")
otu_raw_rare <- dplyr::select(otu_raw_rare, -label, -numOtus)
tvl.rownames <- vars_select(rownames(otu_raw_rare), !contains("Mock")) %>% as.data.frame()
colnames(tvl.rownames) <- "Sample"
div.rare.bl <- rownames_to_column(otu_raw_rare, var = "Sample")
#filter out low read sample
div.rare.bl <- left_join(tvl.rownames, div.rare.bl) %>% column_to_rownames(var = "Sample")
# Generate the rarefaction data, combine with metadata, transform for plotting, and aggregate
rarefy_agg <- rarefy(div.rare.bl, sample = 1000) %>% # subsample
as.data.frame() %>%
rownames_to_column("Sample")
colnames(rarefy_agg) <- c("Sample", "Unique_Otus_per_1k_reads")
rarefy_agg <- left_join(rarefy_agg, dplyr::select(otu_df, Sample, Organ, Sample_Type2), by = "Sample")
rarefy_agg <- dplyr::filter(rarefy_agg, !is.na(Organ))
rarefy_agg_summary  <- rarefy_agg %>%
group_by(Sample_Type2) %>%
summarize(Mean_Species = mean(Unique_Otus_per_1k_reads),
SEM = sqrt(var(Unique_Otus_per_1k_reads)/length(Unique_Otus_per_1k_reads)))
rarefy_agg$Sample_Type2 <- factor(rarefy_agg$Sample_Type2,
levels = c("Empty", "Water",  "AE", "Iso_Ctrl","Syringe_Rinse", "PBS", "Homog_Ctrl","BAL", "Lung",  "Tongue", "Cecum"))
rarefy_agg_summary$Sample_Type2 <- factor(rarefy_agg_summary$Sample_Type2,
levels = c("Empty", "Water",  "AE", "Iso_Ctrl","Syringe_Rinse", "PBS", "Homog_Ctrl","BAL", "Lung",  "Tongue", "Cecum"))
rich_pal <- c( "#EB3D0E", #scarlet - lung
"#0300A6", # blue - BAL
"#FB9986", #vivid tangerine - tongue
"#0BDBA7", #sea foam green - sampling ctrls
"#F2C121", #- isolation ctrls
"#666B6E" #- sequencing ctrls
)
ggplot(rarefy_agg_summary, aes(x = Sample_Type2, y = Mean_Species, fill = Sample_Type2)) + # This brings in the aggregated dataframe and assigns columns to various plot elements
geom_col(color="black", show.legend = FALSE) + # Add bars representing the mean of specimens' diversity
geom_point(data = rarefy_agg, aes(x = Sample_Type2, y = Unique_Otus_per_1k_reads ), width = .2, show.legend = FALSE) + # Use specimen level data, map columns to scatterplot attributes
geom_errorbar(aes(ymin = Mean_Species - SEM, ymax = Mean_Species + SEM), width = 0.5, size=0.7) + # Add an errorbar layer
theme_classic() + # A minimalist theme
theme(axis.title.y = element_text(face = "bold", size = 12),
axis.text.y = element_text(size = 15),
axis.text.x = element_text(size = 10),
axis.title.x = element_text(face = "bold", size = 12),
plot.margin=unit(c(0.5,0.5,0.5,1), "cm")) +
scale_y_continuous(limits = c(0,150), breaks = seq(0,150,25), expand = c(0,0)) +
labs(y = "Richness\n(No. of Unique OTUs/1000 reads)\n", x = "\nSample Type")
otu_df <- mutate(otu_df, Alpha_Div_Groups =  ifelse(test = Organ != "Lung",
yes = as.character(Organ),
no = as.character(Sample_Type))) %>%
dplyr::select(Sample:Gene_16S_copies_per_mL, Alpha_Div_Groups, everything())
otu_df  <- left_join(otu_df, dplyr::select(rarefy_agg, Sample, Unique_Otus_per_1k_reads)) %>%
dplyr::select(Sample:Gene_16S_copies_per_mL, Unique_Otus_per_1k_reads, everything())
tukey_otu_df <- filter(otu_df, Sample_Type!="Mock" & Sample_Type!="Cecum" & Sample_Type!="Tongue")
TukeyHSD(aov(tukey_otu_df[,"Unique_Otus_per_1k_reads"] ~ tukey_otu_df[, "Alpha_Div_Groups"]))
#calculate Shannon diversity indices for all samples
div.shan <- diversity(otu_good, "shannon") %>% data.frame(Shannon = .) %>% rownames_to_column("Sample") %>% left_join(otu_df[,1:6], by = "Sample")
otu_df  <- mutate(otu_df, Shannon = diversity(otu_good, "shannon")) %>% dplyr::select(Sample:Alpha_Div_Groups, Shannon, everything())
# Create summary statistics
agg_shan_all <- div.shan %>%
group_by(Sample_Type2) %>%
summarize(Mean_Sh = mean(Shannon), SEM_Sh = sqrt(var(Shannon)/length(Shannon)))
agg_shan_all$Sample_Type2 <- factor(agg_shan_all$Sample_Type2,
levels = c("Mock", "Empty", "Water",  "AE", "Iso_Ctrl","Syringe_Rinse", "PBS", "Homog_Ctrl","BAL", "Lung",  "Tongue", "Cecum"))
agg_shan_all <- filter(agg_shan_all, Sample_Type2 != "Mock")
div.shan$Sample_Type2 <- factor(div.shan$Sample_Type2, levels = c("Mock", "Empty", "Water",  "AE", "Iso_Ctrl","Syringe_Rinse", "PBS", "Homog_Ctrl","BAL", "Lung",  "Tongue", "Cecum"))
div.shan <- filter(div.shan, Sample_Type2 != "Mock")
ggplot(agg_shan_all, aes(x = Sample_Type2, y = Mean_Sh, fill = Sample_Type2)) + # This brings in the aggregated dataframe and assigns columns to various plot elements
geom_col(color="black", show.legend = FALSE) + # Add bars representing the mean of specimens' Shannon diversity by vendor
geom_point(data = div.shan, aes(x = Sample_Type2, y = Shannon), width = .2, show.legend = FALSE) + # Use specimen level data, map columns to scatterplot attributes
geom_errorbar(aes(ymin = Mean_Sh - SEM_Sh, ymax = Mean_Sh + SEM_Sh), width = 0.5, size=0.7) + # Add an errorbar layer
theme_classic() + # A minimalist theme
theme(axis.title.y = element_text(face = "bold", size = 15),
axis.text.y = element_text(size = 15),
axis.text.x = element_text(size = 10),
axis.title.x = element_text(face = "bold", size = 15),
plot.margin=unit(c(0.5,0.5,0.5,1), "cm")) +
scale_y_continuous(limits = c(-0.01,4.5), breaks = seq(0,4,1), expand = c(0,0)) +
labs(y = "Shannon Diversity Index\n", x = "\nSample Type")
tukey_otu_df <- dplyr::filter(otu_df, Sample_Type != "Mock" & Sample_Type != "Cecum" & Sample_Type != "Tongue")
TukeyHSD(aov(tukey_otu_df[,"Shannon"] ~ tukey_otu_df[, "Alpha_Div_Groups"]))
bray_df_lbtce <- dplyr::filter(otu_df, Sample_Type2 == "Lung" | Sample_Type2 == "Tongue" | Sample_Type2 == "Cecum" | Sample_Type2 == "BAL" | Sample_Type2 == "Empty") %>%
droplevels() %>%
column_to_rownames("Sample")  %>% # This turns the specimen column into the rownames.
dplyr::select(-c(Sample_Type:Alpha_Div_Groups, Shannon))
#identify all unique comparisons within sample types
sample_df <- rownames(bray_df_lbtce) %>% as.data.frame()
colnames(sample_df) <- c("Sample")
sample_df <- mutate(sample_df, Sample_Type = factor(str_extract(Sample, pattern="Lung|BAL|Tongue|Cecum|Blank")))
dim(sample_df) #[1] 87  2
comb_all <- combinations(n = 87, r = 2, v = sample_df$Sample) %>% as.data.frame()
colnames(comb_all) <- c("Sample", "Comparison")
comb_all_mut <- mutate(comb_all, Sample_Type = factor(str_extract(Sample, pattern="Lung|BAL|Tongue|Cecum|Blank")),
Comparison_Type = factor(str_extract(Comparison, pattern="Lung|BAL|Tongue|Cecum|Blank")))
comb_uniq <- dplyr::filter(comb_all_mut, (Sample_Type == Comparison_Type))
#calculate distance, convert to symmetric matrix, and calculate row means (one mean per sample)
bray_dist_lbtce <- vegdist(bray_df_lbtce, method = "bray")
bray_dist_lbtce <- as.matrix(bray_dist_lbtce) %>% as.data.frame() %>% rownames_to_column(var="Sample")
bray_dist_lbtce <- dplyr::select(bray_dist_lbtce, Sample, everything())
bray_dist_lbtce_long <- pivot_longer(bray_dist_lbtce, -Sample, names_to = "Comparison", values_to = "BC_Index")
bray_dist_lbtce_long <- dplyr::filter(bray_dist_lbtce_long, Sample != Comparison)
bray_dist_lbtce_long_mut_filt_uniq <- left_join(comb_uniq, bray_dist_lbtce_long)
bray_dist_lbtce_summary <- group_by(bray_dist_lbtce_long_mut_filt_uniq, Sample) %>%
summarize(Mean_BC = mean(BC_Index), SEM = sem(BC_Index))
bray_dist_lbtce_summary <- ungroup(bray_dist_lbtce_summary) %>%
mutate(Sample_Type = factor(str_extract(Sample,pattern="Lung|BAL|Tongue|Cecum|Blank")))
bray_dist_lbtce_summary$Sample_Type <- factor(bray_dist_lbtce_summary$Sample_Type,
levels = c("Blank", "BAL","Lung", "Tongue", "Cecum"))
pdf("Fig3_081820_D0.pdf", useDingbats = FALSE, width=8, height = 5)
ggplot(bray_dist_lbtce_summary, aes(x=Sample_Type, y = Mean_BC)) +
geom_boxplot(aes(fill=Sample_Type), show.legend = FALSE) +
geom_point() +
scale_y_continuous(limits = c(0,1.01), expand = c(0,0)) +
theme_classic() +
theme(axis.title.y = element_text(face = "bold", size = 15),
axis.text.y = element_text(size = 15),
axis.text.x = element_text(size = 10),
axis.title.x = element_text(face = "bold", size = 15),
plot.margin=unit(c(0.5,0.5,0.5,1), "cm")) +
labs(y = "Bray-Curtis Dissimilarity\n", x = "\nSample Type") +
scale_fill_manual(values=c("#DBE6EC","#0300A6", "#EB3D0E", "#0BDBA7", "#0BDBA7"))
dev.off()
View(otu_df)
View(ddpcr)