UYÊN BÉO Ú Ụ
Hoang et al
2024-08-07
Custom script
indices_normality <- function(rich, nrow, ncol) {
### p-value < 0.05 means data failed normality test
par(mfrow = c(nrow, ncol))
for (i in names(rich)) {
shap <- shapiro.test(rich[, i])
qqnorm(rich[, i], main = i, sub = shap$p.value)
qqline(rich[, i])
}
par(mfrow = c(1, 1))
}#ARGs Abundance
library(pheatmap)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
#Import file input
log_abun_data <- read.csv("/Users/qcvp/R/mẫu plume báo trong nước/log10.txt", header = T, sep = "\t")
log_abun_data2 <- as.matrix(log_abun_data[,-1])
rownames(log_abun_data2) <- log_abun_data[,1]
log_abun_data3 <- as.matrix(log_abun_data2[,-c(1,8)])
log_abun_data4 <- as.matrix(log_abun_data2[,c(1,8)])
#Heatmap with Intl
pheatmap::pheatmap(log_abun_data3, cluster_rows = TRUE,cluster_cols = TRUE) 
#Heatmap without Intl
pheatmap::pheatmap(log_abun_data4, cluster_rows = TRUE,cluster_cols = FALSE) 
PCOA
Plot
library(vegan)## Loading required package: permute## Loading required package: latticelibrary(ape)##
## Attaching package: 'ape'## The following object is masked from 'package:dplyr':
##
## whererel.bray <- vegdist(log_abun_data2, method = "bray")
pcoa.sub <- pcoa(rel.bray)
pcoa_coord <- pcoa.sub$vectors[,1:2]
trans_data <- read.table("/Users/qcvp/R/mẫu plume báo trong nước/Transect.txt", header = T, sep = "\t")
rel_pcoa <- cbind(pcoa_coord, trans_data)
a1 <- ggplot() +
geom_hline(yintercept=0, colour="lightgrey", linetype = 2) +
geom_vline(xintercept=0, colour="lightgrey", linetype = 2) +
geom_point(data = rel_pcoa, aes(x=Axis.1, y=Axis.2, color = Station.Substrate, shape = Transect), alpha = 0.7, size = 5) +
geom_point(data = rel_pcoa, aes(x=Axis.1, y=Axis.2), colour = "white", size = 0.5) +
xlab(paste("PCo1 (", round(pcoa.sub$values$Relative_eig[1]*100, 1), "%)")) +
ylab(paste("PCo2 (", round(pcoa.sub$values$Relative_eig[2]*100, 1), "%)")) +
scale_shape_manual(values = c(15, 16, 17,18)) +
theme_bw() +
coord_equal() +
labs(color = "Sample types") +
scale_colour_manual(breaks =c ("Sediment", "Water"),
values = c("brown", "blue")) +
theme(axis.title.x = element_text(size=15), # remove x-axis labels
axis.title.y = element_text(size=15), # remove y-axis labels
panel.background = element_blank(),
panel.grid.major = element_blank(), #remove major-grid labels
panel.grid.minor = element_blank(), #remove minor-grid labels
plot.background = element_blank())
a1
##Permanova test
adonis2(rel.bray ~ Transect, data = rel_pcoa,permutations = 1000, method ="bray")## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = rel.bray ~ Transect, data = rel_pcoa, permutations = 1000, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Model 3 0.06970 0.13062 1.0016 0.4196
## Residual 20 0.46392 0.86938
## Total 23 0.53362 1.00000adonis2(rel.bray ~ Station.Substrate, data = rel_pcoa,permutations = 1000, method ="bray")## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = rel.bray ~ Station.Substrate, data = rel_pcoa, permutations = 1000, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.34551 0.64748 40.408 0.000999 ***
## Residual 22 0.18811 0.35252
## Total 23 0.53362 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Barplot of abundances
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithslibrary(stringr)
# Load data
Abs_abun_arg <- read.table("/Users/qcvp/R/mẫu plume báo trong nước/log10.txt", header = TRUE, sep = "\t")
# Reshape the data to long format
Abs_abun_arg_long <- melt(Abs_abun_arg, id.vars = "ID", variable.name = "ARG", value.name = "Abundance")
Abs_abun_arg_long <- Abs_abun_arg_long %>%
mutate(Type = str_sub(ID, 4, 4)) %>%
mutate(Transect = str_sub(ID, 1, 2))
#Calculate errors
Abs_abun_summary <- Abs_abun_arg_long %>%
group_by(Transect, ARG, Type) %>%
summarise(
Mean = mean(Abundance, na.rm = TRUE),
SD = sd(Abundance, na.rm = TRUE))## `summarise()` has grouped output by 'Transect', 'ARG'. You can override using
## the `.groups` argument.Abs_abun_summary## # A tibble: 64 × 5
## # Groups: Transect, ARG [32]
## Transect ARG Type Mean SD
## <chr> <fct> <chr> <dbl> <dbl>
## 1 BD X16S S 7.02 0.117
## 2 BD X16S W 7.7 0.384
## 3 BD Sul1 S 6.68 0.633
## 4 BD Sul1 W 3.18 0.261
## 5 BD Sul2 S 7.78 0.706
## 6 BD Sul2 W 4.64 0.52
## 7 BD tetQ S 6.41 0.439
## 8 BD tetQ W 4.19 0.473
## 9 BD tetX S 7.46 0.868
## 10 BD tetX W 4.49 0.307
## # ℹ 54 more rows#Only Sed
Abs_abun_arg_sed <-Abs_abun_summary[Abs_abun_summary$Type == "S",]
Abs_abun_arg_wat <-Abs_abun_summary[Abs_abun_summary$Type== "W",]
# Plot using ggplot2
# Plot for sediment data
s <- ggplot(Abs_abun_arg_sed, aes(x = Transect, y = Mean, fill = Transect)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymin = Mean-SD, ymax = Mean+SD), width = 0.2, position = position_dodge(0.9)) +
theme_minimal() +
theme(legend.text = element_text(size = 15),
legend.position = "none",
plot.title = element_text(size = 20, face = "bold"),
legend.title = element_text(size = 18),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.text.y = element_text(size = 15),
strip.text = element_text(size = 15),
axis.text.x = element_text(size = 15)) +
labs(x = "ARGs", y = "Mean Absolute Abundance (Log10)", title = "Sediment") +
facet_grid(~ ARG, scales = "free")
s
# Plot for water data
w <- ggplot(Abs_abun_arg_wat, aes(x = Transect, y = Mean, fill = Transect)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymin = Mean-SD, ymax = Mean+SD), width = 0.2, position = position_dodge(0.9)) +
theme_minimal() +
theme(legend.text = element_text(size = 15),
legend.position = "none",
plot.title = element_text(size = 20, face = "bold"),
legend.title = element_text(size = 18),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.text.y = element_text(size = 15),
strip.text = element_text(size = 15),
axis.text.x = element_text(size = 15)) +
labs(x = "ARGs", y = "Mean Absolute Abundance (Log10)", title = "Water") +
facet_grid(~ ARG, scales = "free")
w
#Statistical tests ## Load the data
# Transform the data to long format for plotting
Abs_abun_arg_long <- Abs_abun_arg %>%
gather(key = "Gene", value = "Abundance", -ID)Checking distribution and add Group column
#Add a column to identify sample groups
Abs_abun_arg_long <- Abs_abun_arg_long %>%
mutate(Type = str_sub(ID, 4, 4)) %>%
mutate(Transect = str_sub(ID, 1, 2))
#Create a gene list
genes <- c("Sul1", "Sul2", "tetQ", "tetM", "tetX", "VanA", "Intl")
# Checking sample distribution for each Gene
for (gene in genes) {
# Filter data for the current gene
gene_log10 <- Abs_abun_arg_long[Abs_abun_arg_long$Gene == gene, ]
# Checking distribution
gene_log10 |>
dplyr::select(Abundance) |>
indices_normality(nrow = 1, ncol =1)}
Test between water and sediment for each gene
# Loop through each gene
for (gene in genes) {
# Posthoc Dunn test for each gene
test_wilcox <- function(gene) {
test <- ggpubr::compare_means(Abundance ~ Type,
data = filter(Abs_abun_arg_long, Gene == gene),
method = "wilcox.test")
return(test)
}
wilcox_test_log10 <- lapply(unique(Abs_abun_arg_long$Gene), test_wilcox)
}
# Display results
wilcox_test_log10## [[1]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.00510 0.0051 0.0051 ** Wilcoxon
##
## [[2]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.000000740 0.00000074 7.4e-07 **** Wilcoxon
##
## [[3]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.000000740 0.00000074 7.4e-07 **** Wilcoxon
##
## [[4]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.000476 0.00048 0.00048 *** Wilcoxon
##
## [[5]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.264 0.26 0.26 ns Wilcoxon
##
## [[6]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.432 0.43 0.43 ns Wilcoxon
##
## [[7]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.0000359 0.000036 3.6e-05 **** Wilcoxon
##
## [[8]]
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Abundance S W 0.000000740 0.00000074 7.4e-07 **** WilcoxonTest between Transect
# Loop through each gene
for (gene in genes) {
# Calculate statistical significance based on Kruskal-Wallis test
kruskal_test <- function(gene) {
KW_result <- kruskal.test(Abundance ~ Transect, data = filter(Abs_abun_arg_long, Gene == gene))
return(KW_result) }
KW_result_log10 <- lapply(unique(Abs_abun_arg_long$Gene), kruskal_test)
# Posthoc Dunn test for each gene
test_dunn <- function(gene) {
dunn_result <- FSA::dunnTest(Abundance ~ Transect,
data = filter(Abs_abun_arg_long, Gene == gene),
method = "bh")
return(dunn_result) }
dunn_result_log10 <- lapply(unique(Abs_abun_arg_long$Gene), test_dunn)
}## Registered S3 methods overwritten by 'FSA':
## method from
## confint.boot car
## hist.boot car## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.
## Warning: Transect was coerced to a factor.# Kruskal-Wallis
KW_result_log10## [[1]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 3.1097, df = 3, p-value = 0.375
##
##
## [[2]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 1.02, df = 3, p-value = 0.7964
##
##
## [[3]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 3.1267, df = 3, p-value = 0.3725
##
##
## [[4]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 3.1564, df = 3, p-value = 0.3681
##
##
## [[5]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 5.1867, df = 3, p-value = 0.1586
##
##
## [[6]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 9.3344, df = 3, p-value = 0.02516
##
##
## [[7]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 4.5499, df = 3, p-value = 0.2079
##
##
## [[8]]
##
## Kruskal-Wallis rank sum test
##
## data: Abundance by Transect
## Kruskal-Wallis chi-squared = 0.46, df = 3, p-value = 0.9276# Dunn test
dunn_result_log10## [[1]]## Dunn (1964) Kruskal-Wallis multiple comparison## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR -1.3475123 0.1778153 0.5334459
## 2 BD - HM 0.2245854 0.8223018 0.8223018
## 3 CR - HM 1.5720977 0.1159279 0.6955673
## 4 BD - NP -0.7554236 0.4499949 0.6749923
## 5 CR - NP 0.5920888 0.5537912 0.6645494
## 6 HM - NP -0.9800090 0.3270817 0.6541634
##
## [[2]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR 0.3265986 0.7439715 0.8927658
## 2 BD - HM 0.9389711 0.3477456 1.0000000
## 3 CR - HM 0.6123724 0.5402914 1.0000000
## 4 BD - NP 0.6940221 0.4876684 1.0000000
## 5 CR - NP 0.3674235 0.7133032 1.0000000
## 6 HM - NP -0.2449490 0.8064959 0.8064959
##
## [[3]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR -0.04082483 0.9674355 0.9674355
## 2 BD - HM 1.51051867 0.1309111 0.3927334
## 3 CR - HM 1.55134350 0.1208194 0.7249163
## 4 BD - NP 0.48989795 0.6242061 0.7490473
## 5 CR - NP 0.53072278 0.5956109 0.8934163
## 6 HM - NP -1.02062073 0.3074342 0.6148683
##
## [[4]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR 0.4695876 0.6386497 0.7663796
## 2 BD - HM 1.3475123 0.1778153 0.5334459
## 3 CR - HM 0.8779247 0.3799846 0.5699769
## 4 BD - NP 1.5312640 0.1257042 0.7542249
## 5 CR - NP 1.0616764 0.2883826 0.5767653
## 6 HM - NP 0.1837517 0.8542083 0.8542083
##
## [[5]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR 2.0592384 0.03947141 0.2368284
## 2 BD - HM 1.8720349 0.06120178 0.1836053
## 3 CR - HM -0.1872035 0.85150109 0.8515011
## 4 BD - NP 1.3104244 0.19005229 0.3801046
## 5 CR - NP -0.7488140 0.45396934 0.6809540
## 6 HM - NP -0.5616105 0.57438145 0.6892577
##
## [[6]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR -0.9669810 0.333553529 0.50033029
## 2 BD - HM -0.6789441 0.497173275 0.59660793
## 3 CR - HM 0.2880369 0.773318496 0.77331850
## 4 BD - NP 1.8105176 0.070215553 0.14043111
## 5 CR - NP 2.7774986 0.005477907 0.03286744
## 6 HM - NP 2.4894617 0.012793669 0.03838101
##
## [[7]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR 0.2860849 0.77481307 0.7748131
## 2 BD - HM 1.8799867 0.06010989 0.3606593
## 3 CR - HM 1.5939018 0.11095804 0.3328741
## 4 BD - NP 1.2669476 0.20517404 0.4103481
## 5 CR - NP 0.9808626 0.32666049 0.4899907
## 6 HM - NP -0.6130391 0.53985046 0.6478205
##
## [[8]]## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Benjamini-Hochberg method.## Comparison Z P.unadj P.adj
## 1 BD - CR 0.1224745 0.9025233 0.9025233
## 2 BD - HM 0.6123724 0.5402914 1.0000000
## 3 CR - HM 0.4898979 0.6242061 1.0000000
## 4 BD - NP 0.4082483 0.6830914 1.0000000
## 5 CR - NP 0.2857738 0.7750514 1.0000000
## 6 HM - NP -0.2041241 0.8382565 1.0000000ARG resistance heatmap
ENtero
library(reshape)##
## Attaching package: 'reshape'## The following objects are masked from 'package:reshape2':
##
## colsplit, melt, recast## The following objects are masked from 'package:tidyr':
##
## expand, smiths## The following object is masked from 'package:dplyr':
##
## renamelibrary(ggplot2)
ksddata <- read.csv("/Users/qcvp/R/mẫu plume báo trong nước/entero.txt", header = T, sep = "\t")
for (x in 3:23) {
ksddata[,x] <- gsub( 'R', '1', ksddata[,x])
ksddata[,x] <- gsub( 'S', '0', ksddata[,x])
ksddata[,x] <- gsub( 'I', '0.5', ksddata[,x])
ksddata[,x] <- gsub( 'NA', '2', ksddata[,x])}
ksddata2 <- as.matrix(ksddata[,-1])
rownames(ksddata2) <- ksddata2[,1]
ksddata2 <- as.matrix(ksddata2[,-1])
data_long <- melt(ksddata2)## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE#Add group of antitiobics
data_long$group <- NA
data_long$group[data_long$X2 %in% c('TCC','PTZ', 'AUG', 'PIP', 'TIC')] <- 'Penicillin'
data_long$group[data_long$X2 %in% c('FEP', 'CAZ' , 'CTX', 'FOX')] <- 'Cephalosporins'
data_long$group[data_long$X2 %in% c('IPM', 'ETP')] <- 'Carbapanems'
data_long$group[data_long$X2 %in% c('CN', 'AK', 'TOB')] <- 'Aminoglycosides'
data_long$group[data_long$X2 %in% c('CIP', 'LEV')] <- 'Fluoroquiolones'
data_long$group[data_long$X2 %in% c('ATM')] <- 'Monobactams'
data_long$group[data_long$X2 %in% c('SXT')] <- 'Miscellaneous Antibiotics'
data_long$group[data_long$X2 %in% c('FF', 'TGC', 'T')] <- 'Other'
#Plotting
data_long$group <- factor(data_long$group, levels = c("Penicillin", "Cephalosporins", "Carbapanems", "Aminoglycosides", "Fluoroquiolones", "Monobactams", "Miscellaneous Antibiotics", "Other"))
data_long$X2 <- factor(data_long$X2, levels = c("PTZ", "TCC", "CAZ", "CTX", "FEP", "IPM", "AK", "CN", "TOB", "CIP", "LEV", "ATM", "SXT", "AUG", "PIP", "TIC", "FOX", "ETP", "FF", "TGC", "T"))
kd_plot <- ggplot(data_long, aes(x = X2, y = X1, fill = value)) +
geom_tile(colour = "black") +
scale_fill_manual(values = c('1' = 'red', '0.5' = 'yellow', '0' = 'green', '2' = 'gray')) +
theme_minimal() +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = "none",
axis.text.x = element_text(size=13, angle = 45, hjust = 1),
axis.text.y = element_text(size=13),
legend.title = element_blank(),
legend.text = element_blank(),
panel.spacing = unit(0, "lines"),
strip.text = element_text(size=13))+
facet_grid(~group, scales = "free", space = "free")
kd_plot
Gram duonwg
library(reshape)
library(ggplot2)
ksddata <- read.csv("/Users/qcvp/R/mẫu plume báo trong nước/gramduong.txt", header = T, sep = "\t")
for (x in 3:18) {
ksddata[,x] <- gsub( 'R', '1', ksddata[,x])
ksddata[,x] <- gsub( 'S', '0', ksddata[,x])
ksddata[,x] <- gsub( 'I', '0.5', ksddata[,x])
ksddata[,x] <- gsub( 'NA', '2', ksddata[,x])}
ksddata2 <- as.matrix(ksddata[,-1])
rownames(ksddata2) <- ksddata2[,1]
ksddata2 <- as.matrix(ksddata2[,-1])
data_long <- melt(ksddata2)## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE#Add group of antitiobics
data_long$group <- NA
data_long$group[data_long$X2 %in% c( 'AUG', 'TIC')] <- 'Penicillin'
data_long$group[data_long$X2 %in% c('E')] <- 'Macrolides'
data_long$group[data_long$X2 %in% c('IPM')] <- 'Carbapanems'
data_long$group[data_long$X2 %in% c('CN', 'K')] <- 'Aminoglycosides'
data_long$group[data_long$X2 %in% c('CIP', 'LEV', 'OFX')] <- 'Fluoroquiolones'
data_long$group[data_long$X2 %in% c('V')] <- 'Glycopeptide'
data_long$group[data_long$X2 %in% c('SXT')] <- 'Miscellaneous Antibiotics'
data_long$group[data_long$X2 %in% c('L')] <- 'Oxazolidinones'
data_long$group[data_long$X2 %in% c('DA', 'T', 'RA', 'FF')] <- 'Other'
#Plotting
data_long$group <- factor(data_long$group, levels = c("Penicillin", "Macrolides", "Carbapanems" ,"Oxazolidinones", "Aminoglycosides", "Fluoroquiolones", "Glycopeptide", "Miscellaneous Antibiotics", NA , "Other"))
data_long$X2 <- factor(data_long$X2, levels = c( "AUG", "TIC", "E", "IPM", "CN", "K", "CIP", "LEV", "OFX", "V", "SXT", "L", "FF", "T", "DA", "RA"))
kd_plot <- ggplot(data_long, aes(x = X2, y = X1, fill = value)) +
geom_tile(colour = "black") +
scale_fill_manual(values = c('1' = 'red', '0.5' = 'yellow', '0' = 'green', '2' = 'gray')) +
theme_minimal() +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = "none",
axis.text.x = element_text(size=13, angle = 45, hjust = 1),
axis.text.y = element_text(size=13),
legend.title = element_blank(),
legend.text = element_blank(),
panel.spacing = unit(0, "lines"),
strip.text = element_text(size=13))+
facet_grid(~group, scales = "free", space = "free")
kd_plot
#MAR plot ##Preparing data
library(tidyr)
MAR_data <- read.table("/Users/qcvp/R/mẫu plume báo trong nước/mar.txt", sep = "\t", header = T)
names(MAR_data)## [1] "Sample" "Taxonomy" "AUG"
## [4] "TCC" "TIC" "PTZ"
## [7] "PIP" "E" "IMP"
## [10] "CN" "TOB" "K"
## [13] "CIP" "LEV" "OFX"
## [16] "ATM" "V" "L"
## [19] "DA" "R" "T"
## [22] "SXT" "FF" "CAZ"
## [25] "CTX" "FOX" "FEP"
## [28] "ETP" "AK" "TGC"
## [31] "Transect" "Station.Substrate"row.names(MAR_data) <- MAR_data$SampleID
MAR_data2 <- MAR_data[,-1]
MAR_data2 <- MAR_data2 %>% separate(Taxonomy, c("Genus", "Species"))##Plotting
library(ggplot2)
# Adding color for plotting
custom_colors <- c( "#CBD588", "#5F7FC7","#DA5724", "#508578", "#CD9BCD",
"#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
"#8569D5", "#5E738F","#D1A33D", "#8A7C64", "#599861", "black", "red", "green")
# MAR calculate
MAR_data2$MAR <- apply(MAR_data2[,3:30], 1, function(x) (length(which(x == "R"))/length(which(x %in% c("R", "S", "I")))))
# Plotting
box3 = ggplot(MAR_data2, aes(x = Genus , y = MAR, color = Genus)) +
geom_jitter(position = position_jitter(width = .20), alpha = 0.5, size = 3) + theme_bw() +
geom_boxplot(aes(colour = Genus), alpha=0.1, outlier.colour = NA) +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size = 20),
axis.text.y = element_text(size=20),
axis.text.x = element_blank()) +
ylab("MAR index") + scale_fill_manual(values = custom_colors) + geom_hline(yintercept = 0.2, linetype="dashed") +
facet_wrap(~Transect, nrow=1, scales = "free_x")
box3
#Bacterial composition plot
library(dplyr)
library(reshape)
library(stringr)
library(ggplot2)
MAR_data <- MAR_data %>%
mutate(Abundance = 1/nrow(MAR_data))
#Plotting
total <-ggplot(MAR_data, aes(x = Station.Substrate,y = Abundance, fill = Taxonomy)) +
geom_bar(stat = "identity", position ="stack") + scale_color_brewer("Set1") +
theme_bw() +
ggtitle("Total bacterial isolates") +
theme(legend.text = element_text(face = "italic")) +
theme(axis.title.x = element_blank(),
#legend.position = "none",
legend.text = element_text(size = 15),
legend.title = element_text(size = 18),
axis.title.y = element_text(size = 18),
axis.text.y = element_text(size=15),
axis.text.x = element_text(angle=0, colour = "black", size=18),
axis.ticks.x = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(), #remove major-grid labels
panel.grid.minor = element_blank())
total