I need to organize and comment the code, but could host this at github. Or could provide as supplemental. Can only share patient data and supplementary tables, not raw microbiome, 16S or liver RNAseq. Raw Data available ate EGA. Enriched data available upon request (ps).
init
library(ggplot2)
library(ggsci)
library(ggpubr)
library(tidyverse)
library(lmerTest)
library(vegan)
library(DESeq2)
library(DT)
library(RColorBrewer)
library(reshape2)
library(rio)
library(biomaRt)
library(enrichR)
library(phyloseq)
library(gghalves)
library(ggpmisc)load data
ethanol_data <- rio::import_list("./Data/Patient_data.xlsx")
# ps data objects
ps.SI.16S <- readRDS("./Data/ps.BARIA.SI.16S.RDS")
ps.BARIA <- readRDS("./Data/BARIA.metagenomics.RDS")
ps.Myco <- readRDS("Data/BARIA.mycobiome.RDS")
ps.intervention <- readRDS("./Data/ps.ETHANASH.RDS")
ps.RNAseq <- readRDS("./Data/ps.BARIA.RNASeq.liver.RDS")
ps.BARIA.ko <- readRDS("./Data/BARIA.metagenomics.ko.RDS")
ps.RNAseq <- prune_samples(!sample_names(ps.RNAseq) %in% c("SB20","SB4"), ps.RNAseq)
#### global plotting settings ####
colvec <- readRDS("Data/named.colvec.RDS")
shape_vec <- c(19,17,15)
names(shape_vec) <- names(colvec)[1:3]
theme_set(theme_classic2())
logscale=F
legendpos="bottom"
stat_compare_mean_method = "t.test"fix some data
ethanol_data[["BARIA"]]$NAFLD <- factor(ethanol_data[["BARIA"]]$NAFLD, levels=names(colvec))
ethanol_data[["Antwerp"]]$NAFLD <- factor(ethanol_data[["Antwerp"]]$NAFLD, levels=names(colvec))
ethanol_data[["ETHANASH"]]$NAFLD <- factor(ethanol_data[["ETHANASH"]]$NAFLD, levels=names(colvec))
ps.BARIA.ko@sam_data$NAFLD <- factor(ps.BARIA.ko@sam_data$NAFLD, levels=names(colvec))transform data
ps.BARIA@sam_data$NAFLD <- factor(ps.BARIA@sam_data$NAFLD, levels=names(colvec)[1:3])
ps.BARIA.compositional <- transform_sample_counts(ps.BARIA, function(x) x/sum(x))
ps.BARIA.fam <- microbiome::aggregate_taxa(ps.BARIA.compositional, "Family")top13.fams <- names(sort(taxa_sums(ps.BARIA.fam), decreasing = T)[2:14])
colvec.fam <- c(RColorBrewer::brewer.pal(name = "Paired", n=12),"#000000")
names(colvec.fam) <- top13.famsS1 - Portal
df.portal.merged <- ethanol_data$BARIA[,colnames(ethanol_data$BARIA) %in% colnames(ethanol_data$Antwerp)]
df.portal.merged$Cohort <- "Prospective"
df.portal.merged.A <- ethanol_data$Antwerp[,colnames(ethanol_data$Antwerp) %in% colnames(ethanol_data$BARIA)]
df.portal.merged.A$Cohort <- "Validation"
df.portal.merged <- rbind(df.portal.merged, df.portal.merged.A[,colnames(df.portal.merged)])
df.portal.merged$Ethanol <- df.portal.merged$Portal
df.portal.merged <- df.portal.merged[!is.na(df.portal.merged$Ethanol),]
df.portal.merged$Steatosis <- factor(df.portal.merged$Steatosis)
df.portal.merged$Ballooning <- factor(df.portal.merged$Ballooning)
df.portal.merged$Inflammation <- factor(df.portal.merged$Inflammation)
df.portal.merged$Fibrosis <- factor(df.portal.merged$Fibrosis)
df.portal.merged$NAS <- factor(df.portal.merged$NAS)
df.portal.merged$SAFA <- factor(df.portal.merged$SAFA)
#levels(df.portal.merged$Fibrosis)[levels(df.portal.merged$Fibrosis) %in% c("0","1")] <- "0/1"
levels(df.portal.merged$Fibrosis)[levels(df.portal.merged$Fibrosis) %in% c("2","3","4","5")] <- "2+"
levels(df.portal.merged$Steatosis)[levels(df.portal.merged$Steatosis) %in% c("2","3")] <- "2/3"
levels(df.portal.merged$Ballooning)[levels(df.portal.merged$Ballooning) %in% c("1","2")] <- "1/2"
levels(df.portal.merged$Inflammation)[levels(df.portal.merged$Inflammation) %in% c("1","2")] <- "1/2"
levels(df.portal.merged$NAS)[!levels(df.portal.merged$NAS) %in% c("0","1")] <- "2+"
levels(df.portal.merged$SAFA)[!levels(df.portal.merged$SAFA) %in% c("0","1")] <- "2+"
pseudo_porta <- sort(unique(df.portal.merged$Ethanol))[2]/2
portal_plots <- function(df, mc, x){
p <- ggplot(df[!is.na(df$Ethanol) & !is.na(get(x, df)),], aes_string(x=x, y="Ethanol + pseudo_porta")) +
geom_boxplot(outlier.shape = NA) +
ggbeeswarm::geom_beeswarm() +
scale_color_manual(values = colvec) +
scale_y_log10() +
labs(x=NULL, y="Ethanol (mM)") +
#theme(axis.text.x = element_text(angle = 20, hjust = 0)) +
stat_compare_means(comparisons = mc, method = stat_compare_mean_method,label = "p.signif") +
#theme(axis.text.x = element_text(angle = 0)) +
NULL
p = p + facet_wrap(~Cohort)
return(p)
}
mc <- list(c(names(colvec)[1],names(colvec)[2]),c(names(colvec)[1],names(colvec)[3]))
pS2.A <- portal_plots(df = df.portal.merged, mc=mc, x = "NAFLD")+
theme(axis.text.x = element_text(angle = 30, vjust = 0.7)) +
NULL; plot(pS2.A)
mc.S <- list(c("0","1"),c("0","2/3"))
pS2.B <- portal_plots(df = df.portal.merged, mc=mc.S, x = "Steatosis") +
NULL; plot(pS2.B)
mc.I <- list(c("0","1/2"))
pS2.C <- portal_plots(df = df.portal.merged, mc=mc.I, x = "Inflammation") +
NULL; plot(pS2.C)
mc.B <- list(c("0","1/2"))
pS2.D <- portal_plots(df = df.portal.merged, mc=mc.B, x = "Ballooning") +
NULL; plot(pS2.D)
mc.F <- list(c("0","1"),c("0","2+"))
pS2.E <- portal_plots(df = df.portal.merged, x = "Fibrosis", mc = mc.F) +
NULL; plot(pS2.E)
mc.F <- list(c("0","1"),c("0","2+"))
pS2.F <- portal_plots(df = df.portal.merged, x = "NAS", mc = mc.F) +
NULL; plot(pS2.F)
mc.F <- list(c("0","1"),c("0","2+"))
pS2.G <- portal_plots(df = df.portal.merged, x = "SAFA", mc = mc.F) +
NULL; plot(pS2.G)
S2.panel <- ggarrange(plotlist = list(pS2.B + labs(title="Steatosis"),
pS2.C + labs(title="Inflammation"),
pS2.D + labs(title="Ballooning"),
pS2.E + labs(title="Fibrosis"),
pS2.A + labs(title="NAFLD"),
pS2.G + labs(title="SAF-A"),
pS2.F + labs(title="NAS")),
labels = c("A", "B", "C", "D","E","F","G"),
ncol = 4, nrow = 2); plot(S2.panel)
ggsave(plot = S2.panel, filename = "./Figures/Figure.S1.Portal.NALF.strat.pdf", height = 8, width = 12)
ggsave(plot = S2.panel, filename = "./Figures/Figure.S1.Portal.NALF.strat.png", height = 8, width = 12)Figure 1
theme_set(theme_classic2())
logscale=F
legendpos="bottom"
stat_compare_mean_method = "t.test"
#### fix some data ####
ethanol_data[["BARIA"]]$NAFLD <- factor(ethanol_data[["BARIA"]]$NAFLD, levels=names(colvec))
ethanol_data[["Antwerp"]]$NAFLD <- factor(ethanol_data[["Antwerp"]]$NAFLD, levels=names(colvec))
ethanol_data[["ETHANASH"]]$NAFLD <- factor(ethanol_data[["ETHANASH"]]$NAFLD, levels=names(colvec))
#### Numbers in paper:
# portal means
aggregate(ethanol_data$BARIA$Portal, by=list(ethanol_data$BARIA$NAFLD), function(x) mean(x, na.rm = T))## Group.1 x
## 1 No Steatosis 2.460714
## 2 NAFL 27.408333
## 3 NASH 33.746000aggregate(ethanol_data$Antwerp$Portal, by=list(ethanol_data$Antwerp$NAFLD), function(x) mean(x, na.rm = T))## Group.1 x
## 1 No Steatosis 2.107430
## 2 NAFL 8.068738
## 3 NASH 20.972855# portal medians
aggregate(ethanol_data$BARIA$Portal, by=list(ethanol_data$BARIA$NAFLD), function(x) median(x, na.rm = T))## Group.1 x
## 1 No Steatosis 0.735
## 2 NAFL 23.780
## 3 NASH 22.130aggregate(ethanol_data$Antwerp$Portal, by=list(ethanol_data$Antwerp$NAFLD), function(x) median(x, na.rm = T))## Group.1 x
## 1 No Steatosis 1.721820
## 2 NAFL 1.755293
## 3 NASH 11.642723# post prandial increase
sum(ethanol_data$BARIA$Post_Prandial>ethanol_data$BARIA$Fasted, na.rm = T)## [1] 101sum(!ethanol_data$BARIA$Post_Prandial>ethanol_data$BARIA$Fasted, na.rm = T)## [1] 8df.MMT <- ethanol_data$BARIA[!is.na(ethanol_data$BARIA$Post_Prandial),] %>% pivot_longer(values_to = "Ethanol",cols = c("Post_Prandial","Fasted"), names_to = "Time")
summary(lmer(Ethanol~Time + (1 | Subject_ID), df.MMT))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Ethanol ~ Time + (1 | Subject_ID)
## Data: df.MMT
##
## REML criterion at convergence: -815
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.65330 -0.49371 -0.02662 0.33395 2.82552
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject_ID (Intercept) 0.0008217 0.02867
## Residual 0.0007062 0.02657
## Number of obs: 218, groups: Subject_ID, 109
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.072e-02 3.744e-03 1.675e+02 13.547 < 2e-16 ***
## TimePost_Prandial 3.396e-02 3.600e-03 1.080e+02 9.433 9.08e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TmPst_Prndl -0.481lmer_res <- summary(lmer(Ethanol~Time*NAFLD + (1 | Subject_ID), df.MMT));lmer_res## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Ethanol ~ Time * NAFLD + (1 | Subject_ID)
## Data: df.MMT
##
## REML criterion at convergence: -847.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.57606 -0.45503 -0.04023 0.36054 2.85873
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject_ID (Intercept) 0.0004845 0.02201
## Residual 0.0006057 0.02461
## Number of obs: 218, groups: Subject_ID, 109
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.011e-02 5.035e-03 1.770e+02 7.965 1.93e-13 ***
## TimePost_Prandial 1.958e-02 5.308e-03 1.060e+02 3.689 0.000357 ***
## NAFLDNAFL 1.284e-02 6.669e-03 1.770e+02 1.925 0.055874 .
## NAFLDNASH 4.727e-02 1.210e-02 1.770e+02 3.906 0.000133 ***
## TimePost_Prandial:NAFLDNAFL 1.898e-02 7.030e-03 1.060e+02 2.700 0.008069 **
## TimePost_Prandial:NAFLDNASH 5.386e-02 1.276e-02 1.060e+02 4.221 5.14e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TmPs_P NAFLDNAF NAFLDNAS TP_P:NAFLDNAF
## TmPst_Prndl -0.527
## NAFLDNAFL -0.755 0.398
## NAFLDNASH -0.416 0.219 0.314
## TP_P:NAFLDNAF 0.398 -0.755 -0.527 -0.166
## TP_P:NAFLDNAS 0.219 -0.416 -0.166 -0.527 0.314coefficients(lmer_res)["TimePost_Prandial","Estimate"] + coefficients(lmer_res)["TimePost_Prandial:NAFLDNAFL","Estimate"]## [1] 0.03856491coefficients(lmer_res)["TimePost_Prandial","Estimate"] + coefficients(lmer_res)["TimePost_Prandial:NAFLDNASH","Estimate"]## [1] 0.07344111as.numeric(signif(coefficients(lmer_res)[,5], digits =1))+1## [1] 1.00000 1.00040 1.06000 1.00010 1.00800 1.00005# median/IQR increase
summary(ethanol_data$BARIA$Post_Prandial - ethanol_data$BARIA$Fasted)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.04922 0.01177 0.02300 0.03396 0.04436 0.13990 37#### Prospective ####
mc <- list( c(names(colvec)[1], names(colvec)[2]), c(names(colvec)[1], names(colvec)[3]) )
p1 <- ggplot(ethanol_data[["BARIA"]], aes(x=NAFLD, y=Portal, color=NAFLD)) +
geom_boxplot() +
ggbeeswarm::geom_beeswarm() +
scale_color_manual(values = colvec) +
labs(x=NULL, y="Ethanol (mM)") +
#scale_y_log10() +
#stat_compare_means(method = "aov") +
#stat_compare_means(comparisons = mc, method = "wilcox.test") +
stat_compare_means(comparisons = mc, method = stat_compare_mean_method,label = "p.signif", label.y = c(110,120)) +
coord_cartesian(ylim = c(0.1, 130)) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
NULL; plot(p1)
#### Validation ####
mc <- list( c(names(colvec)[1], names(colvec)[2]), c(names(colvec)[1], names(colvec)[3]) )
p2 <- ggplot(ethanol_data[["Antwerp"]], aes(x=NAFLD, y=Portal, color=NAFLD)) +
stat_compare_means(comparisons = mc, label = "p.signif", label.y = c(110,120)) +
geom_boxplot() +
ggbeeswarm::geom_beeswarm() +
scale_color_manual(values = colvec) +
labs(x=NULL, y="Ethanol (mM)") +
#scale_y_log10() +
#stat_compare_means(method = "aov") +
coord_cartesian(ylim = c(0.1, 130)) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
NULL; plot(p2)
#### MMT ####
df.MMT <- ethanol_data$BARIA %>% pivot_longer(values_to = "Ethanol",cols = c("Post_Prandial","Fasted"), names_to = "Time")
df.MMT$Dummy <- NA
df.MMT$Dummy[df.MMT$NAFLD==names(colvec)[1]] <- 1
df.MMT$Dummy[df.MMT$NAFLD==names(colvec)[2]] <- 3
df.MMT$Dummy[df.MMT$NAFLD==names(colvec)[3]] <- 5
df.MMT$Dummy[df.MMT$Time=="Post_Prandial"] <- df.MMT$Dummy[df.MMT$Time=="Post_Prandial"]+1
df.MMT$Dummy <- factor(df.MMT$Dummy)
mc3 <- list(c("2","4"),c("2","6"),c("1","3"),c("1","5"))
p3 <- ggplot(df.MMT, aes(x=Dummy, y=Ethanol, color=NAFLD, shape=Time)) +
geom_boxplot(outlier.shape = NA) +
ggbeeswarm::geom_beeswarm() +
scale_color_manual(values = colvec) +
#facet_wrap(~Group) +
#scale_y_log10() +
labs(x=NULL, y="Ethanol (mM)") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
geom_line(aes(group=Subject_ID), color="black", alpha=0.2) +
stat_compare_means(comparisons = mc3, method = "wilcox",label = "p.signif", label.y = c(0.24,0.26,0.20,0.22)) +
#stat_compare_means(comparisons = mc4, method = stat_compare_mean_method,label = "p.signif", label.y = c(-0.7,-0.6)) +
# stat_compare_means(comparisons = mc4, method = stat_compare_mean_method,label = "p.signif") +
coord_cartesian(ylim = c(min(df.MMT$Ethanol, na.rm = T), 0.28)) +
theme(strip.text.x = element_blank()) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
NULL; plot(p3)
#### ETHANASH ####
range_limit <- 0.4
ethanol_data$ETHANASH_Ethanol$Group <- factor(ethanol_data$ETHANASH_Ethanol$Group, levels=names(colvec))
ethanol_data$ETHANASH_Ethanol$Time2 <- "Fasted"
ethanol_data$ETHANASH_Ethanol$Time2[ethanol_data$ETHANASH_Ethanol$Time>0] <- "Post prandial"
df.ETHANASH_MMT <- ethanol_data$ETHANASH_Ethanol[ethanol_data$ETHANASH_Ethanol$Invervention %in% c("MMT"),]
df.ETHANASH_MMT$Ethanol[df.ETHANASH_MMT$Ethanol<0] <- 0
df.ETHANASH_MMT$Ethanol2 <- df.ETHANASH_MMT$Ethanol
df.ETHANASH_MMT$Ethanol2[df.ETHANASH_MMT$Ethanol>range_limit] <- range_limit
df.ETHANASH_MMT$out_of_bounds <- df.ETHANASH_MMT$Ethanol>range_limit
df.ETHANASH_MMT$out_of_bounds[df.ETHANASH_MMT$Time>0] <- "Post Prandial"
df.ETHANASH_MMT$out_of_bounds[df.ETHANASH_MMT$Time<=0] <- "Fasted"
df.ETHANASH_MMT$out_of_bounds[df.ETHANASH_MMT$Ethanol>range_limit] <- "Out of range"
p4 <- ggplot(df.ETHANASH_MMT, aes(x=Time, y=Ethanol, color=Group, group=Group)) +
geom_smooth(aes(group=Group), show.legend = T, se = T, formula = y ~ splines::ns(x), size=1.5) +
#geom_smooth(aes(group=Group)) +
coord_cartesian(ylim = c(-0.01,range_limit+0.05)) +
ggbeeswarm::geom_beeswarm(aes(y=Ethanol2, shape=out_of_bounds), size=1.5, show.legend = F, cex = 0.75) +
#geom_boxplot() +
scale_color_manual(values = colvec) +
scale_fill_manual(values = colvec) +
geom_point(aes(y=Ethanol2, shape=out_of_bounds)) +
scale_shape_manual(values=c(16,3,17))+
#scale_y_log10() +
#facet_wrap(~Group) +
labs(x="Time post prandial (min)", y="Ethanol (mM)") +
theme(legend.position = "none") +
NULL; plot(p4)
# median_data_MMT <- aggregate(df.ETHANASH_MMT$Ethanol, by=list(df.ETHANASH_MMT$Group, df.ETHANASH_MMT$Time, df.ETHANASH_MMT$NAFLD), FUN=median)
#
# df.ETHANASH_MMT$Time_Group <- interaction( df.ETHANASH_MMT$Time, df.ETHANASH_MMT$Group)
#
# p4 <- ggplot(df.ETHANASH_MMT, aes(x=Time, y=Ethanol, color=Group, group=Group)) +
# #geom_smooth(aes(group=Group), show.legend = T, se = T, formula = y ~ splines::ns(x), size=1.5) +
# #geom_smooth(aes(group=Group)) +
# geom_boxplot(aes(group=Time_Group,color=Group)) +
# geom_line(data = median_data_MMT, aes(group=Group.1, x=Group.2, y=x, color=Group.1), inherit.aes = F) +
# coord_cartesian(ylim = c(-0.01,range_limit+0.05)) +
# #ggbeeswarm::geom_beeswarm(aes(y=Ethanol2, shape=out_of_bounds), size=1.5, show.legend = F, cex = 0.75) +
# #geom_boxplot() +
# scale_color_manual(values = colvec) +
# scale_fill_manual(values = colvec) +
# #geom_point(aes(y=Ethanol2, shape=out_of_bounds)) +
# scale_shape_manual(values=c(16,3,17))+
# #scale_y_log10() +
# #facet_wrap(~Group) +
# labs(x="Time post prandial (min)", y="Ethanol (mM)") +
# theme(legend.position = "none") +
# NULL; plot(p4)
#### ETHANASH FMP ####
range_limit_FMP <- 3.1
df.ETHANASH_FMP <- ethanol_data$ETHANASH_Ethanol[!ethanol_data$ETHANASH_Ethanol$Invervention %in% c("MMT"),]
df.ETHANASH_FMP$Ethanol[df.ETHANASH_FMP$Ethanol<0] <- 0
df.ETHANASH_FMP$Ethanol2 <- df.ETHANASH_FMP$Ethanol
df.ETHANASH_FMP$Ethanol2[df.ETHANASH_FMP$Ethanol>range_limit_FMP] <- range_limit_FMP
df.ETHANASH_FMP$out_of_bounds <- df.ETHANASH_FMP$Ethanol>range_limit_FMP
df.ETHANASH_FMP$out_of_bounds[df.ETHANASH_FMP$Time>0] <- "Post Prandial"
df.ETHANASH_FMP$out_of_bounds[df.ETHANASH_FMP$Time<=0] <- "Fasted"
df.ETHANASH_FMP$out_of_bounds[df.ETHANASH_FMP$Ethanol>range_limit_FMP] <- "Out of Range"
df.ETHANASH_FMP$Group[is.na(df.ETHANASH_FMP$Group)] <- "NASH + 4-methylpyrazole + AB"
median_data <- aggregate(df.ETHANASH_FMP$Ethanol, by=list(df.ETHANASH_FMP$Group, df.ETHANASH_FMP$Time, df.ETHANASH_FMP$NAFLD), FUN=median)
df.ETHANASH_FMP$Time_Group <- interaction( df.ETHANASH_FMP$Time, df.ETHANASH_FMP$Group)
p5 <- ggplot(df.ETHANASH_FMP, aes(x=Time, y=Ethanol, color=Group, group=Group)) +
geom_smooth(aes(group=Group), show.legend = T, se = T, formula = y ~ splines::ns(x), size=1.5) +
#geom_boxplot(aes(group=Time_Group,color=Group)) +
#geom_line(data = median_data, aes(group=Group.1, x=Group.2, y=x, color=Group.1), inherit.aes = F) +
coord_cartesian(ylim = c(-0.01,range_limit_FMP)) +
ggbeeswarm::geom_beeswarm(aes(y=Ethanol2, shape=out_of_bounds), size=1.5, show.legend = F, cex = 0.75) +
scale_color_manual(values = colvec) +
scale_fill_manual(values = colvec) +
scale_shape_manual(values=c(16,3,17))+
labs(x="Time post prandial (min)", y="Ethanol (mM)") +
theme(legend.position = "none") +
NULL; plot(p5)
#### shared legend ####
df.leg <- ethanol_data$ETHANASH_Ethanol
df.leg$Time <- as.character(df.leg$Time)
df.leg$Time[!df.leg$Time %in% c("-30","0")] <- "Post prandial"
df.leg$Time[df.leg$Time %in% c("-30","0")] <- "Fasted"
df.leg$Time <- factor(df.leg$Time, levels=c("Fasted","Post prandial","Post Prandial (Out of Bounds)"))
#levels(df.leg$Time)[!levels(df.leg$Time) %in% c("Fasted", "Post prandial")] <- "Post prandial"
#levels(df.leg$Time) <- c(as.character(df.leg$Time))
df.leg <- df.leg[df.leg$Group %in% names(colvec),]
df.leg$Time[1] <- "Post Prandial (Out of Bounds)"
p.legend <- ggplot(df.leg, aes(x=Group, y=Ethanol, color=Group, shape=Time)) +
geom_point(show.legend = T) +
scale_color_manual(values = colvec) +
scale_shape_manual(values=c(16, 17, 3)) +
#scale_y_log10() +
theme(legend.position = legendpos) +
theme(legend.text=element_text(size=10)) +
theme(legend.box="vertical") +
theme(legend.margin=margin()) +
theme(legend.title = element_text(size=12),
legend.key.size = unit(0.06, 'cm'),
legend.key.height= unit(0.06, 'line'),
legend.key.width = unit(0.06, 'line')) +
#guides(shape = guide_legend(override.aes = list(size = 3))) +
#guides(color = guide_legend(override.aes = list(size = 3))) +
# guides(color = guide_legend(order = 1),
# shape = guide_legend(order = 2)) +
NULL;plot(p.legend)
#### aggregate ####
p.m1 <- ggarrange(p1 + labs(y="Ethanol (mM)", title="Prospective portal"),
p2 + labs(y="Ethanol (mM)", title="Validation portal"),
p3 + labs(y="Ethanol (mM)", title="Prospective MMT"), common.legend = T, legend = "none",
nrow = 1,
labels = c("A","B","C")); plot(p.m1)
p.m2 <- ggarrange(p4 + labs(y="Ethanol (mM)", title="Intervention MMT"),
p5 + labs(y="Ethanol (mM)", title="Intervention + AB/4MP"), common.legend = T, legend = "none",
nrow = 1,
labels = c("D","E")); plot(p.m2)
p.m <- ggarrange(p.m1, p.m2, common.legend = T,
nrow = 2,
legend.grob = get_legend(p.legend),
legend = legendpos); plot(p.m)
ggsave(plot = p.m, filename = "Figures/Figure1.Ethanol.Concentrations.pdf", height = 6, width = 9)
ggsave(plot = p.m, filename = "Figures/Figure1.Ethanol.Concentrations.png", height = 6, width = 9)S3 - MMT
df.MMT <- ethanol_data$BARIA[,] %>% pivot_longer(values_to = "Ethanol",cols = c("Post_Prandial","Fasted"), names_to = "Time")
df.MMT$Steatosis <- factor(df.MMT$Steatosis)
df.MMT$Ballooning <- factor(df.MMT$Ballooning)
df.MMT$Inflammation <- factor(df.MMT$Inflammation)
df.MMT$Fibrosis <- factor(df.MMT$Fibrosis)
df.MMT$NAS <- factor(df.MMT$NAS)
df.MMT$SAFA <- factor(df.MMT$SAFA)
#levels(df.MMT$Fibrosis)[levels(df.MMT$Fibrosis) %in% c("0","1")] <- "0/1"
levels(df.MMT$Fibrosis)[levels(df.MMT$Fibrosis) %in% c("2","3","4","5")] <- "2+"
levels(df.MMT$Steatosis)[levels(df.MMT$Steatosis) %in% c("2","3")] <- "2/3"
levels(df.MMT$Ballooning)[levels(df.MMT$Ballooning) %in% c("1","2")] <- "1/2"
levels(df.MMT$NAS)[!levels(df.MMT$NAS) %in% c("0","1")] <- "2+"
levels(df.MMT$SAFA)[!levels(df.MMT$SAFA) %in% c("0","1")] <- "2+"
# add pseudo level for log_y plotting
pseudo_mm <- sort(unique(df.MMT$Ethanol))[2]/2
MMT_plot <- function(df, mc, x){
p <- ggplot(df[!is.na(df$Ethanol) & !is.na(get(x, df)),], aes_string(x=x, y="Ethanol + pseudo_mm")) +
geom_boxplot(outlier.shape = NA) +
ggbeeswarm::geom_beeswarm() +
scale_color_manual(values = colvec) +
facet_wrap(~Time) +
labs(x=NULL, y="Ethanol (mM)") +
stat_compare_means(comparisons = mc, method = stat_compare_mean_method,label = "p.signif") +
scale_y_log10() +
NULL
return(p)
}
mc.N <- list(c(levels(df.MMT$NAFLD)[1],levels(df.MMT$NAFLD)[2]),
c(levels(df.MMT$NAFLD)[1],levels(df.MMT$NAFLD)[3]))
pS1.A <- MMT_plot(df = df.MMT, x = "NAFLD", mc = mc.N)+
theme(axis.text.x = element_text(angle = 30, vjust = 0.7)) +
NULL; plot(pS1.A)
mc.S <- list(c("0","1"),c("0","2/3"))
pS1.B <- MMT_plot(df = df.MMT, x = "Steatosis", mc = mc.S) +
NULL; plot(pS1.B)
mc.I <- list(c("0","1"),c("0","2"))
pS1.C <- MMT_plot(df = df.MMT, x = "Inflammation", mc = mc.I) +
NULL; plot(pS1.C)
mc.B <- list(c("0","1/2"))
pS1.D <- MMT_plot(df = df.MMT, x = "Ballooning", mc = mc.B) +
NULL; plot(pS1.D)
mc.F <- list(c("0","1"),c("0","2+"))
pS1.E <- MMT_plot(df = df.MMT, x = "Fibrosis", mc = mc.F) +
NULL; plot(pS1.E)
pS1.F <- MMT_plot(df = df.MMT, x = "NAS", mc = mc.F) +
NULL; plot(pS1.E)
pS1.G <- MMT_plot(df = df.MMT, x = "SAFA", mc = mc.F) +
NULL; plot(pS1.E)
S1.panel <- ggarrange(plotlist = list(pS1.B + labs(title="Steatosis"),
pS1.C + labs(title="Inflammation"),
pS1.D + labs(title="Ballooning"),
pS1.E + labs(title="Fibrosis"),
pS1.A + labs(title="NAFLD"),
pS1.G + labs(title="SAF-A"),
pS1.F + labs(title="NAS")
),
labels = c("A", "B", "C", "D","E","F","G"),
ncol = 4, nrow = 2); plot(S1.panel)
ggsave(plot = S1.panel, filename = "./Figures/Figure.S2.BARIA.MMT.NALF.strat.pdf", height = 8, width = 12)
ggsave(plot = S1.panel, filename = "./Figures/Figure.S2.BARIA.MMT.NALF.strat..png", height = 8, width = 12)S3 - Insuline
pointsize=2
df.BARIA.insuline.long <- ethanol_data$BARIA %>%
pivot_longer(cols = c("Insuline_120","Insuline_0"), values_to = "Insuline", names_to = "Insuline_Time")%>%
pivot_longer(cols = c("Portal","Fasted","Post_Prandial"), values_to = "Ethanol", names_to = "Plasma_Type")
df.BARIA.insuline.long$Plasma_Type <- factor(df.BARIA.insuline.long$Plasma_Type) %>%
recode_factor(Portal="Portal", Fasted="Fasted", Post.prandial="Post prandial")
df.BARIA.insuline.long$Insuline_Time <- factor(df.BARIA.insuline.long$Insuline_Time) %>%
recode_factor(Insuline_0="Insuline (Fasted)",Insuline_120="Insuline (Post Prandial)", )
Insuline_labs <- c("Insulin Fasted", "Insulin Post Prandial")
p.insuline <- ggplot(df.BARIA.insuline.long, aes(x=Ethanol, y=Insuline, color=NAFLD, shape=NAFLD)) +
geom_point(size=pointsize) +
scale_y_log10() +
scale_x_log10() +
geom_smooth(method = "lm", aes(group=NA), show.legend = F) +
stat_cor(method="s", aes(group=NA), show.legend = F) +
facet_grid(Insuline_Time~Plasma_Type, scales="free",
labeller = labeller(Insuline = Insuline_labs)) +
scale_color_manual(values = colvec[names(colvec) %in% df.BARIA.insuline.long$NAFLD]) +
labs(x="Ethanol (mM)", y="Insuline (pmol/L)") +
NULL; plot(p.insuline)
ggsave(plot = p.insuline, filename = "./Figures/Figure.S3.Insulin.pdf", height = 6, width = 9)
ggsave(plot = p.insuline, filename = "./Figures/Figure.S3.Insulin.png", height = 6, width = 9)S4 - Portal / MMT correlations
corgroup="NAFLD"
df.BARIA.wide <- ethanol_data$BARIA
pS5.A <- ggplot(df.BARIA.wide, aes(x=Portal+sort(unique(df.BARIA.wide$Portal))[2]/2,
y=Fasted+sort(unique(df.BARIA.wide$Fasted))[2]/2,
color=NAFLD,
label=Subject_ID, shape=NAFLD)) +
geom_point() +
geom_smooth(aes_string(group=corgroup), show.legend = F, se = F,method = "lm") +
geom_smooth(aes_string(group=NA), show.legend = F, se = F,method = "lm", col="black") +
scale_color_manual(values = colvec[names(colvec) %in% df.BARIA.wide$NAFLD]) +
stat_cor(aes_string(group=corgroup), show.legend = F) +
stat_cor(aes_string(group=NA), show.legend = F, label.y = 0.12) +
scale_shape_manual(values = shape_vec) +
labs(x="Portal Ethanol (mM)", y="Fasted peripheral Ethanol (mM)") +
NULL; plot(pS5.A)
pS5.B <- ggplot(df.BARIA.wide, aes(x=Portal+sort(unique(Portal))[2]/2,
y=Post_Prandial+sort(unique(Post_Prandial))[2]/2,
color=NAFLD, shape=NAFLD, label=Subject_ID)) +
ggbeeswarm::geom_beeswarm(size=1) +
geom_point() +
geom_smooth(aes_string(group=corgroup), show.legend = F, se = F,method = "lm") +
geom_smooth(aes_string(group=NA), show.legend = F, se = F,method = "lm", col="black") +
scale_color_manual(values = colvec[names(colvec) %in% df.BARIA.wide$NAFLD]) +
stat_cor(aes_string(group=corgroup), show.legend = F, label.y = c(0.18,0.2,0.22)) +
stat_cor(aes_string(group=NA), show.legend = F, label.y = 0.24) +
scale_shape_manual(values = shape_vec) +
labs(x="Portal Ethanol (mM)", y="Post Prandial peripheral Ethanol (mM)") +
#scale_x_log10() +
NULL; plot(pS5.B)
pS5.C <- ggplot(df.BARIA.wide, aes(x=Fasted+sort(unique(Fasted))[2]/2,
y=Post_Prandial+sort(unique(Post_Prandial))[2]/2,
color=NAFLD, shape=NAFLD,
label=Subject_ID)) +
ggbeeswarm::geom_beeswarm(size=1) +
geom_point() +
geom_smooth(aes_string(group=corgroup), show.legend = F, se = F,method = "lm") +
geom_smooth(aes_string(group=NA), show.legend = F, se = F,method = "lm", col="black") +
scale_color_manual(values = colvec[names(colvec) %in% df.BARIA.wide$NAFLD]) +
stat_cor(aes_string(group=corgroup), show.legend = F) +
stat_cor(aes_string(group=NA), show.legend = F, label.y = 0.24) +
scale_shape_manual(values = shape_vec) +
labs(x="Fasted Ethanol (mM)", y="Post Prandial peripheral Ethanol (mM)") +
NULL; plot(pS5.C)
S5.panels <- ggarrange(plotlist = list(pS5.A, pS5.B),labels = c("A","B"),
ncol = 2, nrow = 1,
common.legend = T, legend = "bottom"); plot(S5.panels)
ggsave(plot = S5.panels, filename = "./Figures/Figure.S4.Ant.Ethanol.correlations.pdf", height = 4, width = 9)
ggsave(plot = S5.panels, filename = "./Figures/Figure.S4.Ant.Ethanol.correlations.png", height = 4, width = 9)
#Caption:
# Correlations between ethanol concentrations in portal and fasted plasma (A),
# portal and post prandial plasma (B) and between fasted and post prandial plasma. ST.1 - RNAseq table
colnames(ps.RNAseq@sam_data) <- make.names(colnames(ps.RNAseq@sam_data))
# ps.RNAseq@sam_data$peripheral.ethanol.umol.l <- ps.RNAseq@sam_data$peripheral.ethanol.umol.l/1000
#
# ps.RNAseq@sam_data$post.mmt.ETOH.umol.l <- ps.RNAseq@sam_data$post.mmt.ETOH.umol.l/1000
Enrich.databases <- c("GO_Molecular_Function_2021",
"GO_Biological_Process_2021",
"GO_Cellular_Component_2021",
"KEGG_2021_Human",
"WikiPathways_2019_Human")
deseq.models <- list()
alpha = 0.05
# subsets
ps.porta <- prune_samples(!is.na(ps.RNAseq@sam_data$portal.ethanol.mmol.l),ps.RNAseq)
ps.porta@sam_data$Ethanol <- ps.porta@sam_data$portal.ethanol.mmol.l
ps.fasted <- prune_samples(!is.na(ps.RNAseq@sam_data$peripheral.ethanol.umol.l),ps.RNAseq)
ps.fasted@sam_data$Ethanol <- ps.fasted@sam_data$peripheral.ethanol.umol.l
ps.postp <- prune_samples(!is.na(ps.RNAseq@sam_data$post.mmt.ETOH.umol.l),ps.RNAseq)
ps.postp@sam_data$Ethanol <- ps.postp@sam_data$post.mmt.ETOH.umol.l
ps.nafl <- prune_samples(!is.na(ps.RNAseq@sam_data$NAFLD),ps.RNAseq)
ps.nafl@sam_data$NAFLD <- factor(ps.nafl@sam_data$NAFLD)
# DESEQ models
deseq.models[["Ethanol.Porta"]] <- phyloseq_to_deseq2(ps.porta, ~ Batch + Ethanol)
deseq.models[["Ethanol.Porta"]] <- DESeq(deseq.models[["Ethanol.Porta"]], test="Wald", fitType="parametric")
deseq.models[["Ethanol.fasted"]] <- phyloseq_to_deseq2(ps.fasted, ~ Batch + Ethanol)
deseq.models[["Ethanol.fasted"]] <- DESeq(deseq.models[["Ethanol.fasted"]], test="Wald", fitType="parametric")
deseq.models[["Ethanol.PostP"]] <- phyloseq_to_deseq2(ps.postp, ~ Batch + Ethanol)
deseq.models[["Ethanol.PostP"]] <- DESeq(deseq.models[["Ethanol.PostP"]], test="Wald", fitType="parametric")
deseq.models[["Ethanol.NAFL"]] <- phyloseq_to_deseq2(ps.nafl, ~ Batch + NAFLD)
deseq.models[["Ethanol.NAFL"]] <- DESeq(deseq.models[["Ethanol.NAFL"]], test="Wald", fitType="parametric")
# Results
res.porta = results(deseq.models[["Ethanol.Porta"]], cooksCutoff = FALSE, name = "Ethanol")
res.fasted = results(deseq.models[["Ethanol.fasted"]], cooksCutoff = FALSE, name = "Ethanol")
res.postp = results(deseq.models[["Ethanol.PostP"]], cooksCutoff = FALSE, name = "Ethanol")
res.NAFL = results(deseq.models[["Ethanol.NAFL"]], cooksCutoff = FALSE, name = "NAFLD_1_vs_0")
res.NASH = results(deseq.models[["Ethanol.NAFL"]], cooksCutoff = FALSE, name = "NAFLD_2_vs_0")
# Sigtabs
sigtab.porta = res.porta[which(res.porta$padj < alpha), ]
sigtab.porta$gene <- rownames(sigtab.porta)
sigtab.fasted = res.fasted[which(res.fasted$padj < alpha), ]
sigtab.fasted$gene <- rownames(sigtab.fasted)
sigtab.postp = res.postp[which(res.postp$padj < alpha), ]
sigtab.postp$gene <- rownames(sigtab.postp)
sigtab.NAFL = res.NAFL[which(res.NAFL$padj < alpha), ]
sigtab.NAFL$gene <- rownames(sigtab.NAFL)
sigtab.NASH = res.NASH[which(res.NASH$padj < alpha),]
sigtab.NASH$gene <- rownames(sigtab.NASH)
# EnrichR
Ethanol.Porta.enrich <- enrichr(genes = rownames(sigtab.porta), databases = Enrich.databases)## Uploading data to Enrichr... Done.
## Querying GO_Molecular_Function_2021... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying GO_Cellular_Component_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying WikiPathways_2019_Human... Done.
## Parsing results... Done.Ethanol.Porta.enrich.sig <- do.call(rbind,lapply(Ethanol.Porta.enrich, function(x) x[x$Adjusted.P.value<alpha,]))
Ethanol.fasted.enrich <- enrichr(genes = rownames(sigtab.fasted), databases = Enrich.databases)## Uploading data to Enrichr... Done.
## Querying GO_Molecular_Function_2021... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying GO_Cellular_Component_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying WikiPathways_2019_Human... Done.
## Parsing results... Done.Ethanol.fasted.enrich.sig <- do.call(rbind,lapply(Ethanol.fasted.enrich, function(x) x[x$Adjusted.P.value<alpha,]))
Ethanol.PostP.enrich <- enrichr(genes = rownames(sigtab.postp), databases = Enrich.databases)## Uploading data to Enrichr... Done.
## Querying GO_Molecular_Function_2021... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying GO_Cellular_Component_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying WikiPathways_2019_Human... Done.
## Parsing results... Done.Ethanol.PostP.enrich.sig <- do.call(rbind,lapply(Ethanol.PostP.enrich, function(x) x[x$Adjusted.P.value<alpha,]))
Ethanol.NAFL.enrich <- enrichr(genes = rownames(sigtab.NAFL), databases = Enrich.databases)## Uploading data to Enrichr... Done.
## Querying GO_Molecular_Function_2021... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying GO_Cellular_Component_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying WikiPathways_2019_Human... Done.
## Parsing results... Done.Ethanol.NAFL.enrich.sig <- do.call(rbind,lapply(Ethanol.NAFL.enrich, function(x) x[x$Adjusted.P.value<alpha,]))
Ethanol.NASH.enrich <- enrichr(genes = rownames(sigtab.NASH), databases = Enrich.databases)## Uploading data to Enrichr... Done.
## Querying GO_Molecular_Function_2021... Done.
## Querying GO_Biological_Process_2021... Done.
## Querying GO_Cellular_Component_2021... Done.
## Querying KEGG_2021_Human... Done.
## Querying WikiPathways_2019_Human... Done.
## Parsing results... Done.Ethanol.NASH.enrich.sig <- do.call(rbind,lapply(Ethanol.NASH.enrich, function(x) x[x$Adjusted.P.value<alpha,]))
# Export results
rio::export(list(Porta = data.frame(res.porta),
Fasted = data.frame(res.fasted),
PostPrandial = data.frame(res.postp),
NAFL = data.frame(res.NAFL),
NASH = data.frame(res.NASH),
Porta_EnrichR = Ethanol.Porta.enrich.sig,
Fasted_EnrichR = Ethanol.fasted.enrich.sig,
PostPrandial_EnrichR = Ethanol.PostP.enrich.sig,
NAFL_EnrichR = Ethanol.NAFL.enrich.sig,
NASH_EnrichR = Ethanol.NASH.enrich.sig),
"./Tables/Supplementary.Table.Liver.RNASeq.results.xlsx", row.names=TRUE)S5 - RNAseq MA plots
DESEQ_res <- rio::import_list("./Tables/Supplementary.Table.Liver.RNASeq.results.xlsx", rowNames =T,readxl = FALSE)
alpha=0.05
MA_function <- function(df){
df$Highlight_1 <- df$padj<alpha & !is.na(df$padj)
#df$Highlight_1 <- df$padj<alpha & !is.na(df$padj) & abs(df$log2FoldChange)>fc
df$Highlight_2 <- rownames(df) %in% c("ADH1A","CYP2E1")
df$gene <- rownames(df)
p <- ggplot(data.frame(df), aes(y=log2FoldChange, x=log(baseMean))) +
theme_bw() +
geom_point() +
ggrepel::geom_text_repel(data=subset(df, df$Highlight_1), aes(label=gene), color="blue") +
geom_point(data=subset(df, df$Highlight_1), color="blue") +
geom_point(data=subset(df, df$Highlight_2), color="red") +
ggrepel::geom_text_repel(data=subset(df, df$Highlight_2), aes(label=gene), color="red") +
NULL
return(p)
}
p.MA1 <- MA_function(DESEQ_res$Fasted)
p.MA2 <- MA_function(DESEQ_res$PostPrandial)
p.MA3 <- MA_function(DESEQ_res$Porta)
p.MA4 <- MA_function(DESEQ_res$NAFL)
p.MA5 <- MA_function(DESEQ_res$NASH)
p.merge <- ggarrange(plotlist = list(p.MA1,p.MA2,p.MA3,p.MA4,p.MA5), labels = c("A","B","C","D","E"));p.merge
ggsave(plot = p.merge, filename = "./Figures/Figure.S5.RNASeq_MAPlots.png", height = 14, width = 21)
ggsave(plot = p.merge, filename = "./Figures/Figure.S5.RNASeq_MAPlots.pdf", height = 14, width = 21)S6 - Small Intestinte
sample_names(ps.SI.16S) <- paste0("SB_",ps.SI.16S@sam_data$Subject_ID)
ps.SI.16S@sam_data$Fasted <-ethanol_data$BARIA$Fasted[match(sample_names(ps.SI.16S),ethanol_data$BARIA$Subject_ID)]
ps.SI.16S@sam_data$Post_Prandial <- ethanol_data$BARIA$Post_Prandial[match(sample_names(ps.SI.16S),ethanol_data$BARIA$Subject_ID)]
ps.SI.16S@sam_data$NAFLD <- factor(ps.SI.16S@sam_data$NAFLD.categorisch)
levels(ps.SI.16S@sam_data$NAFLD)[levels(ps.SI.16S@sam_data$NAFLD)==0] <- names(colvec)[1]
levels(ps.SI.16S@sam_data$NAFLD)[levels(ps.SI.16S@sam_data$NAFLD)==1] <- names(colvec)[2]
levels(ps.SI.16S@sam_data$NAFLD)[levels(ps.SI.16S@sam_data$NAFLD)==2] <- names(colvec)[3]
ASV <- taxa_names(ps.SI.16S)[ps.SI.16S@tax_table[,6] %in% "Streptococcus"]
df <- cbind(ps.SI.16S@otu_table[,ASV], ps.SI.16S@sam_data)
df$Streptococcus <- df$ASV_155
df$steatosis <- as.factor(df$steatosis)
summary(aov(Streptococcus~steatosis, df))## Df Sum Sq Mean Sq F value Pr(>F)
## steatosis 1 78551 78551 13.17 0.00248 **
## Residuals 15 89484 5966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(Streptococcus~NAFLD, df))## Df Sum Sq Mean Sq F value Pr(>F)
## NAFLD 2 103040 51520 11.1 0.0013 **
## Residuals 14 64995 4643
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(aov(Streptococcus~NAFLD, df))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Streptococcus ~ NAFLD, data = df)
##
## $NAFLD
## diff lwr upr p adj
## NAFL-No Steatosis 120.0000 26.01120 213.9888 0.0125998
## NASH-No Steatosis 242.3333 96.72652 387.9402 0.0017823
## NASH-NAFL 122.3333 -17.07458 261.7412 0.0893913kruskal.test(df$Streptococcus,df$NAFLD)##
## Kruskal-Wallis rank sum test
##
## data: df$Streptococcus and df$NAFLD
## Kruskal-Wallis chi-squared = 8.0657, df = 2, p-value = 0.01772cor.test(df$Streptococcus,df$EtOH.port, method = "s")##
## Spearman's rank correlation rho
##
## data: df$Streptococcus and df$EtOH.port
## S = 134.94, p-value = 0.07753
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5281821cor.test(df$Streptococcus,df$Fasted, method = "s")##
## Spearman's rank correlation rho
##
## data: df$Streptococcus and df$Fasted
## S = 140.98, p-value = 0.09246
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5070548cor.test(df$Streptococcus,df$Post_Prandial, method = "s")##
## Spearman's rank correlation rho
##
## data: df$Streptococcus and df$Post_Prandial
## S = 24.568, p-value = 0.04963
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7075276p1 <- ggplot(df, aes(y=Streptococcus/500, x=EtOH.port, color=NAFLD, shape=NAFLD)) +
geom_point(size=3) +
geom_smooth(aes(group = 1), method = "lm", show.legend = F) +
scale_color_manual(values = colvec[1:3]) +
stat_cor(aes(group = 1), method = "s") +
scale_shape_manual(values = shape_vec) +
labs(x="Portal Ethanol (mM)", y="Streptococcus Relative Abundance") +
NULL; plot(p1)
p2 <- ggplot(df, aes(y=Streptococcus/500, x=NAFLD, color=NAFLD, shape=NAFLD)) +
geom_boxplot(show.legend = F) +
ggbeeswarm::geom_beeswarm(size=3) +
scale_color_manual(values = colvec[1:3]) +
stat_compare_means(show.legend = F) +
scale_shape_manual(values = shape_vec) +
labs(y="Streptococcus Relative Abundance", x=NULL) +
NULL; plot(p2)
p3 <- ggplot(df, aes(y=Streptococcus/500, x=Fasted, color=NAFLD, shape=NAFLD)) +
geom_point(size=3) +
geom_smooth(aes(group = 1), method = "lm", show.legend = F) +
scale_color_manual(values = colvec[1:3]) +
stat_cor(aes(group = 1), method = "s") +
scale_shape_manual(values = shape_vec) +
labs(x="Fasted Ethanol (mM)", y="Streptococcus Relative Abundance") +
NULL; plot(p3)
p4 <- ggplot(df, aes(y=Streptococcus/500, x=Post_Prandial, color=NAFLD, shape=NAFLD)) +
geom_point(size=3) +
geom_smooth(aes(group = 1), method = "lm", show.legend = F) +
scale_color_manual(values = colvec[1:3]) +
stat_cor(aes(group = 1), method = "s") +
scale_shape_manual(values = shape_vec) +
labs(x="Post prandial Ethanol (mM)", y="Streptococcus Relative Abundance") +
NULL; plot(p4)
p.SI <- ggarrange(plotlist = list(p2, p3, p4, p1), labels = c("A","B","C","D"),
common.legend = T, ncol = 2, nrow = 2,
legend = "right"); plot(p.SI)
ggsave(plot = p.SI, filename = "Figures/Figure.S6.Streptococcus_Small_intestine.pdf", height = 7, width = 7)
ggsave(plot = p.SI, filename = "Figures/Figure.S6.Streptococcus_Small_intestine.png", height = 7, width = 7)S8 - Prospective cohort alpha & beta NAFLD
pointsize=2
ord <- ordinate(ps.BARIA.compositional, "PCoA", "bray")
ps.sub <- prune_samples(!is.na(ps.BARIA.compositional@sam_data$NAFLD), ps.BARIA.compositional)
d <- vegdist(t(ps.sub@otu_table@.Data), "bray")
res <- adonis(d~ps.sub@sam_data$NAFLD)
annotate <- paste0("Adonis p.value:",res$aov.tab$`Pr(>F)`[1]);annotate## [1] "Adonis p.value:0.959"p.b1 <- plot_ordination(ps.sub, ord, color="NAFLD", shape="NAFLD") +
geom_point(size=pointsize) +
scale_color_manual(values=colvec[1:3]) +
annotate("text", x=0, y=-0.2, label=annotate) +
scale_shape_manual(values = shape_vec[1:3]) +
labs(title="NAFLD") +
NULL;p.b1
ps.sub <- prune_samples(!is.na(ps.BARIA.compositional@sam_data$Portal), ps.BARIA.compositional)
d <- vegdist(t(ps.sub@otu_table@.Data), "bray")
res <- adonis(d~ps.sub@sam_data$Portal)
annotate <- paste0("Adonis p.value:",res$aov.tab$`Pr(>F)`[1]);annotate## [1] "Adonis p.value:0.513"p.b2 <- plot_ordination(ps.sub, ord, color="Portal") +
scale_color_viridis_c(option="B", end=0.8) +
geom_point(size=pointsize) +
annotate("text", x=0, y=-0.2, label=annotate) +
guides(fill=guide_legend(title="Ethanol mM")) +
labs(title="Portal") +
NULL;p.b2
ps.sub <- prune_samples(!is.na(ps.BARIA.compositional@sam_data$Fasted), ps.BARIA.compositional)
d <- vegdist(t(ps.sub@otu_table@.Data), "bray")
res <- adonis(d~ps.sub@sam_data$Fasted)
annotate <- paste0("Adonis p.value:",res$aov.tab$`Pr(>F)`[1]);annotate## [1] "Adonis p.value:0.236"p.b3 <- plot_ordination(ps.sub, ord, color="Fasted") +
geom_point(size=pointsize) +
scale_color_viridis_c(option="B", end=0.8) +
annotate("text", x=0, y=-0.2, label=annotate) +
guides(fill=guide_legend(title="Ethanol mM")) +
labs(title="Fasted") +
NULL;p.b3
ps.sub <- prune_samples(!is.na(ps.BARIA.compositional@sam_data$Post.prandial), ps.BARIA.compositional)
d <- vegdist(t(ps.sub@otu_table@.Data), "bray")
res <- adonis(d~ps.sub@sam_data$Post.prandial)
annotate <- paste0("Adonis p.value:",res$aov.tab$`Pr(>F)`[1]);annotate## [1] "Adonis p.value:0.208"p.b4 <- plot_ordination(ps.BARIA.compositional, ord, color="Post.prandial") +
geom_point(size=pointsize) +
scale_color_viridis_c(option="B", end=0.8) +
annotate("text", x=0, y=-0.2, label=annotate) +
#guides(color=guide_legend(title="Ethanol mM")) +
labs(title="Post Prandial") +
NULL;p.b4
p.b2$labels$colour <- 'Ethanol (mMM)'
p.b3$labels$colour <- 'Ethanol (mMM)'
p.b4$labels$colour <- 'Ethanol (mMM)'
p.beta <- ggarrange(plotlist = list(p.b1, p.b2, p.b3, p.b4),
ncol = 2, nrow = 2,
labels = c("A","B","C","D"), common.legend = F);p.beta
ggsave(plot = p.beta, filename = "./Figures/Figure.S8.Beta.Diversity.pdf", height = 6, width = 9)
ggsave(plot = p.beta, filename = "./Figures/Figure.S8.Beta.Diversity.png", height = 6, width = 9)S7 - Prospective cohort alpha & beta ethanol
rich <- estimate_richness(ps.BARIA.compositional, measures = "Shannon")
df <- ps.BARIA.compositional@sam_data
df$Shannon <- as.vector(rich$Shannon)
pointsize=1.5
pa1 <- ggplot(df, aes(x=NAFLD, y=Shannon, color=NAFLD, shape=NAFLD)) +
gghalves::geom_half_violin(side = "left", show.legend = F) +
gghalves::geom_half_point(side = "right", size=pointsize) +
labs(x=NULL, y="Shannon Index") +
scale_color_manual(values = colvec[1:3]) +
scale_y_continuous(limits = c(1.5, 2.8)) +
stat_compare_means(show.legend = F) +
labs(title="NAFLD") +
NULL;pa1
pa2 <- ggplot(df, aes(x=Portal, y=Shannon)) +
geom_point(size=pointsize, aes(color=NAFLD, shape=NAFLD)) +
labs(x="Ethanol (mM)", y="Shannon Index", title = "Portal") +
scale_color_manual(values = colvec[1:3]) +
scale_y_continuous(limits = c(1.5, 2.8)) +
stat_cor() +
NULL
pa3 <- ggplot(df, aes(x=Fasted/1000, y=Shannon)) +
geom_point(size=pointsize, aes(color=NAFLD, shape=NAFLD)) +
labs(x="Ethanol (mM)", y="Shannon Index", title = "Fasted") +
scale_color_manual(values = colvec[1:3]) +
scale_y_continuous(limits = c(1.5, 2.8)) +
stat_cor() +
NULL; pa3
pa4 <- ggplot(df, aes(x=Post.prandial/1000, y=Shannon)) +
geom_point(size=pointsize, aes(color=NAFLD, shape=NAFLD)) +
labs(x="Ethanol (mM)", y="Shannon Index", title = "Post prandial") +
scale_color_manual(values = colvec[1:3]) +
scale_y_continuous(limits = c(1.5, 2.8)) +
stat_cor() +
NULL; pa4
p.alpha <- ggarrange(plotlist = list(pa1, pa2, pa3, pa4),
ncol = 2, nrow = 2,
labels = c("A","B","C","D"), common.legend = T, legend = "bottom");p.alpha
ggsave(plot = p.alpha, filename = "./Figures/Figure.S7.BARIA.Diversity.Ethanol.pdf", height = 6, width = 9)
ggsave(plot = p.alpha, filename = "./Figures/Figure.S7.BARIA.Diversity.Ethanol.png", height = 6, width = 9)S9 - Prospective Mycobiome
pointsize=1.5
ps.Saccho <- prune_taxa(taxa_names(ps.Myco)[grep("Saccharomycetaceae",taxa_names(ps.Myco))], ps.Myco)
df <- ps.Saccho@sam_data
df$Fungi <- sample_sums(ps.Saccho)/df$reads_total
colvec_myco <- colvec[1:3]
df$NAFLD <- droplevels(factor(df$NAFLD, levels=names(colvec_myco)))
mc <- list( c(names(colvec)[1], names(colvec)[2]), c(names(colvec)[1], names(colvec)[3]) )
pa1 <- ggplot(df, aes(x=NAFLD, y=Fungi, color=NAFLD, shape=NAFLD)) +
gghalves::geom_half_violin(side = "left", show.legend = F) +
gghalves::geom_half_point(side = "right", size=pointsize) +
labs(x=NULL, y="Fungi") +
scale_color_manual(values = colvec_myco) +
#scale_y_continuous(limits = c(1.5, 2.8)) +
scale_y_log10(limits=c(0.00000001, 0.006)) +
#scale_y_log10() +
#stat_compare_means(show.legend = F, comparisons = mc, method = "wilcox",label = "p.signif",
# label.y = c(-3, -2.6)) +
stat_compare_means(show.legend = F, comparisons = mc, method = "wilcox",label = "p.signif",
label.y = c(-3, -2.6)) +
labs(title="NAFLD") +
NULL;pa1
pa2 <- ggplot(df, aes(x=Portal, y=Fungi)) +
geom_point(size=pointsize, aes(color=NAFLD, shape=NAFLD)) +
labs(x="Ethanol (mM)", y="Fungi", title = "Portal") +
scale_color_manual(values = colvec_myco) +
#scale_y_continuous(limits = c(1.5, 2.8)) +
scale_y_log10(limits=c(0.00000001, 0.006)) +
stat_cor() +
NULL; pa2
pa3 <- ggplot(df, aes(x=Fasted/1000, y=Fungi)) +
geom_point(size=pointsize, aes(color=NAFLD, shape=NAFLD)) +
labs(x="Ethanol (mM)", y="Fungi", title = "Fasted") +
scale_color_manual(values = colvec_myco) +
#scale_y_continuous(limits = c(1.5, 2.8)) +
scale_y_log10(limits=c(0.00000001, 0.006)) +
stat_cor() +
NULL; pa3
pa4 <- ggplot(df, aes(x=`Post prandial`/1000, y=Fungi)) +
geom_point(size=pointsize, aes(color=NAFLD, shape=NAFLD)) +
labs(x="Ethanol (mM)", y="Fungi", title = "Post prandial") +
scale_color_manual(values = colvec_myco) +
#scale_y_continuous(limits = c(1.5, 2.8)) +
scale_y_log10(limits=c(0.00000001, 0.006)) +
stat_cor() +
NULL; pa4
p.fungi <- ggarrange(plotlist = list(pa1, pa2, pa3, pa4),
ncol = 2, nrow = 2,
labels = c("A","B","C","D"), common.legend = T, legend = "bottom");p.fungi
ggsave(plot = p.fungi, filename = "./Figures/Figure.S9.Fungi.pdf", height = 6, width = 9)
ggsave(plot = p.fungi, filename = "./Figures/Figure.S9.Fungi.png", height = 6, width = 9)S10 - Intervention Microbiome
ps.intervention.t0 <- prune_samples(ps.intervention@sam_data$Time_Point==0, ps.intervention)
ps.intervention.t0 <- prune_taxa(taxa_sums(ps.intervention.t0)>0, ps.intervention.t0)
ps.intervention.deseq <- ps.intervention.t0
ps.intervention.deseq@sam_data$Group <- as.character(ps.intervention.deseq@sam_data$Group)
diagdds = phyloseq_to_deseq2(ps.intervention.deseq, ~ Group)
diagdds = DESeq(diagdds, test="Wald", fitType="parametric")
resultsNames(diagdds)## [1] "Intercept" "Group_NASH_vs_Healthy"res = results(diagdds,
cooksCutoff = FALSE)
title <- gsub(".*\\): ","",res@elementMetadata$description[2])
res$padj <- p.adjust(res$pvalue, "BH")
res <- cbind(res, ps.intervention.deseq@tax_table[rownames(res),]@.Data)
res$mean_abundance <- (taxa_sums(ps.intervention)/nsamples(ps.intervention))[rownames(res)]
y <- datatable(data.frame(res),
filter = 'top',
options = list(pageLength = 10, scipen=3),
rownames = T) %>%
formatSignif(columns = c('baseMean', 'log2FoldChange', 'lfcSE','pvalue', 'padj', 'stat'), digits = 3)
phyla_abundance <- aggregate(res$mean_abundance, by=list(res$Phylum), FUN=sum)
topx <- phyla_abundance[tail(order(phyla_abundance$x), n=5),1]
topx <- topx[topx != "UNMAPPED"]
df.res.sub <- data.frame(res[res$Phylum %in% topx,])
Order_abundance <- aggregate(df.res.sub$mean_abundance, by=list(df.res.sub$Order), FUN=sum)
topx <- Order_abundance[tail(order(Order_abundance$x), n=9),1]
df.res.sub$Order[!df.res.sub$Order %in% topx] <- "Minor orders"
df.res.sub$Order <- factor(df.res.sub$Order, levels=c(rev(topx),"Minor orders"))p.taxa <- ggplot(df.res.sub, aes(x=log2FoldChange, y=-log(padj), color=Order, fill=Order, size=log(mean_abundance))) +
theme_bw() +
geom_point(alpha=1,shape = 21, colour = "black") +
facet_wrap(~Phylum) +
geom_vline(xintercept = -1, alpha=0.5) +
geom_vline(xintercept = 1, alpha=0.5) +
geom_hline(yintercept = -log(0.05), alpha=0.5) +
guides(color=guide_legend(ncol=2)) +
#labs(title = paste0("Differential abundant taxa NASH vs Healthy")) +
scale_size(guide="none") +
guides(color=guide_legend(ncol=1)) +
labs(y="p.adjusted (-log10)") +
scale_color_manual(values = c(brewer.pal(n=9,name = "Spectral"),"#FFFFFF")) +
scale_fill_manual(values = c(brewer.pal(n=9,name = "Spectral"),"#FFFFFF")) +
guides(fill = guide_legend(override.aes = list(size = 5))) +
NULL; plot(p.taxa)
Group="NAFLD"
Subject_ID="Subject_ID" # paired analysis
Time_Point="Time_Point" # time serie analysis
measures = c("Observed","Shannon")
df <- estimate_richness(ps.intervention.t0, measures = c("Shannon", "Observed"))
df <- cbind(df, ps.intervention.t0@sam_data)
vars <- c(Group, Subject_ID, Time_Point)
df.ad <- cbind(estimate_richness(ps.intervention.t0, measures = measures),ps.intervention.t0@sam_data[,c(vars)])
df.ad.long <- melt(df.ad, id.vars = vars)
p.alpha <- hospital_names <- c(
`Observed` = "Observed",
`Shannon` = "Shannon",
`FPD` = "Phylogenetic Diversity",
`InvSimpson` = "Inverse Simpson"
)
p.alpha <- ggplot(df.ad.long, aes(x=NAFLD, y=value, group=NAFLD, color=NAFLD, shape=NAFLD)) +
geom_boxplot(outlier.colour = NA) +
ggbeeswarm::geom_beeswarm(size=2, cex = 1,priority = "random") +
facet_wrap(~variable, scale="free_y", labeller = as_labeller(hospital_names)) +
scale_color_manual(values=colvec[c(1,3)]) +
labs(x="Time Point",y="") +
theme_bw() +
#geom_line(aes(group=Relatives), color="black") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
labs(x=NULL) +
stat_compare_means() +
scale_shape_manual(values = shape_vec[c(1,3)]) +
NULL; plot(p.alpha)
ps.intervention.t0.trans <- transform_sample_counts(ps.intervention.t0, function(x) x/sum(x))
d <- phyloseq::distance(ps.intervention.t0.trans, "bray")
adonis <- vegan::adonis(d~ps.intervention.t0.trans@sam_data$NAFLD)
ord <- ordinate(ps.intervention.t0.trans, "PCoA", "bray")
p.beta <- plot_ordination(ps.intervention.t0.trans, ord, color="NAFLD", shape="NAFLD") +
theme_classic2() +
scale_color_manual(values=colvec) +
geom_point(size=2) +
annotate("text",x = 0.05, y=-0.2, label=paste0("Adonis p.value:",adonis$aov.tab$`Pr(>F)`[1])) +
#ggrepel::geom_label_repel(aes(label=Subject_ID)) +
scale_shape_manual(values = shape_vec[c(1,3)]) +
NULL; plot(p.beta)
adonis <- vegan::adonis(d~ps.intervention.t0.trans@sam_data$Delta)
ps.intervention.t0.trans@sam_data$'Ethanol increase (mM/h)' <- ps.intervention.t0.trans@sam_data$Delta
p.beta2 <- plot_ordination(ps.intervention.t0.trans, ord, color="Delta", shape="NAFLD") +
theme_classic2() +
#scale_color_manual(values=colvec) +
scale_color_viridis_c(option="B", end=0.8, name="Ethanol \n increase \n (mM/h)") +
#guides(color=guide_legend(title="'Ethanol increase (mM/h)'")) +
#scale_coloe=_discrete(name = "New Legend Title")
geom_point(size=2) +
annotate("text",x = 0.05, y=-0.2, label=paste0("Adonis p.value:",adonis$aov.tab$`Pr(>F)`[1])) +
theme(legend.position = "right") +
scale_shape_manual(values = shape_vec[c(1,3)]) +
NULL; plot(p.beta2)
p.ab <- ggarrange(p.alpha, p.beta, common.legend = T, legend = "right", labels = c("A","B")); p.ab
p.ab2 <- ggarrange(p.beta2, p.taxa, common.legend = F, labels = c("C","D")); p.ab2
p.abt <- ggarrange(p.ab, p.ab2, common.legend = F, ncol = 1); p.abt
ggsave(plot = p.abt, "./Figures/Figure.S10.ETHANASH.metagenome.pdf", height = 9, width = 9)
ggsave(plot = p.abt, "./Figures/Figure.S10.ETHANASH.metagenome.png", height = 9, width = 9)ST.2 - Correlations and differential abundance Table
correlations_wrapper <- function(ps, var){
ps.cor <- prune_samples(!is.na(get(var, ps@sam_data)), ps)
ps.cor <- prune_taxa(taxa_sums(ps.cor)!=0, ps.cor)
corlist <- list()
for(taxa in rownames(ps.cor@tax_table)){
corlist[[taxa]] <- cor.test(t(ps.cor@otu_table[taxa,]@.Data), get(var, ps.cor@sam_data), method = "spearman")
}
df.cor <- data.frame(ps.cor@tax_table@.Data)
df.cor$p.value <- unlist(lapply(corlist, function(x) x$p.value))[rownames(df.cor)]
df.cor$estimate <- unlist(lapply(corlist, function(x) unname(x$estimate)))[rownames(df.cor)]
df.cor$Mean_Abundance <- (taxa_sums(ps)/sum(ps.cor@otu_table))[rownames(df.cor)]
df.cor$padj <- p.adjust(df.cor$p.value, method = "fdr", n=ntaxa(ps))
return(df.cor)
}
df.cor.fasted <- correlations_wrapper(ps.BARIA.compositional, "Fasted")
df.cor.postp <- correlations_wrapper(ps.BARIA.compositional, "Post.prandial")
df.cor.portal <- correlations_wrapper(ps.BARIA.compositional, "Portal")
ps.BARIA@sam_data$NAFLD2 <- ps.BARIA@sam_data$NAFLD==names(colvec)[1]
deseq.models <- list()
deseq.models.nafld <- phyloseq_to_deseq2(ps.BARIA, ~ NAFLD2)
deseq.models.nafld <- DESeq(deseq.models.nafld, test="Wald", fitType="parametric")
resultsNames(deseq.models.nafld)## [1] "Intercept" "NAFLD2TRUE"res = results(deseq.models.nafld, cooksCutoff = FALSE, name = "NAFLD2TRUE")
df.deseq <- data.frame(res)
df.deseq$Family <- as.vector(ps.BARIA.compositional@tax_table[rownames(df.deseq),5])
df.deseq$Phylum <- as.vector(ps.BARIA.compositional@tax_table[rownames(df.deseq),2])
baria.res.list <- list("Correlations.fasted"=df.cor.fasted,
"Correlations.postprandial"=df.cor.postp,
"Correlations.portal"=df.cor.portal,
"DESeq_NAFLD"=df.deseq)
rio::export(baria.res.list, file = "./Tables/Supplementary.Table.BARIA.Metagenomics.xlsx", rowNames=T)S11 - Correlations vulcano/t(MA)
rho_p_postp = min(df.cor.postp[df.cor.postp$padj<0.05,]$estimate)
rho_abundance_plot <- function(df, select, rho_p){
p.cortop <- ggplot(df, aes(x=estimate, y=Mean_Abundance*100)) +
geom_point(aes(color=Family), size=4) +
scale_y_log10(limits=c(0.01, 4)) +
scale_color_manual(values = colvec.fam) +
#ggrepel::geom_text_repel(data=df.cor[!df.cor$Species %in% select, ], aes(label=Species)) +
#ggrepel::geom_text_repel(data=df.cor[df.cor$Species %in% select, ], aes(label=Species), col="red") +
ggrepel::geom_text_repel(
data=df[df$Species %in% select, ],
aes(
label=Species
),
#col="red",
force = 0.4,
nudge_x = 0.1,
direction = "y",
hjust = 0,
#segment.size = 0.5,
segment.curvature = -0,
box.padding = 0.0,
segment.alpha=0.2,
size=3.5
) +
labs(x="Spearman's Rho", y="Mean relative abundance (%)") +
NULL
if (!is.null(rho_p)){p.cortop = p.cortop + geom_vline(xintercept = rho_p, col="blue", alpha=0.4)}
return (p.cortop)
}
cor_vulcano_plot <- function(df, select, rho_p){
p.vulc <- ggplot(df, aes(x=estimate, y=-log(padj))) +
geom_point(aes(color=Family), size=3) +
#scale_y_log10(limits=c(0.01, 4)) +
scale_color_manual(values = colvec.fam) +
#ggrepel::geom_text_repel(data=df.cor[!df.cor$Species %in% select, ], aes(label=Species)) +
#ggrepel::geom_text_repel(data=df.cor[df.cor$Species %in% select, ], aes(label=Species), col="red") +
ggrepel::geom_text_repel(
data=df[df$Species %in% select, ],
aes(
label=Species
),
#col="red",
force = 0.4,
nudge_x = 0.1,
direction = "y",
hjust = 0,
#segment.size = 0.5,
segment.curvature = -0,
box.padding = 0.0,
segment.alpha=0.2,
size=3.5
) +
labs(x="Spearman's Rho", y="-log(p.adjusted)") +
NULL
if (!is.null(rho_p)){p.vulc = p.vulc + geom_vline(xintercept = rho_p, col="blue", alpha=0.4)}
return (p.vulc)
}
p.cor.A <- rho_abundance_plot(df.cor.postp, select=NULL, rho_p=rho_p_postp) +
scale_x_continuous(limits = c(-0.3,0.6)); plot(p.cor.A)
p.cor.B <- rho_abundance_plot(df.cor.portal, select=NULL, rho_p=NULL) +
scale_x_continuous(limits = c(-0.3,0.6)); plot(p.cor.B)
p.cor.C <- rho_abundance_plot(df.cor.fasted, select=NULL, rho_p=NULL) +
scale_x_continuous(limits = c(-0.3,0.6)); plot(p.cor.C)
p.vul.A <- cor_vulcano_plot(df.cor.postp[df.cor.postp$Family %in% top13.fams & df.cor.postp$Mean_Abundance>0.0001,], select=NULL, rho_p=NULL) + geom_hline(yintercept = -log(0.05), col="blue") +
scale_y_continuous(limits = c(0,6)) + scale_x_continuous(limits = c(-0.3,0.6)); plot(p.vul.A)
p.vul.B <- cor_vulcano_plot(df.cor.portal[df.cor.portal$Family %in% top13.fams & df.cor.portal$Mean_Abundance>0.0001,], select=NULL, rho_p=NULL) + geom_hline(yintercept = -log(0.05), col="blue") +
scale_y_continuous(limits = c(0,6)) + scale_x_continuous(limits = c(-0.3,0.6)); plot(p.vul.B)
p.vul.C <- cor_vulcano_plot(df.cor.fasted[df.cor.fasted$Family %in% top13.fams & df.cor.fasted$Mean_Abundance>0.0001,], select=NULL, rho_p=NULL) + geom_hline(yintercept = -log(0.05), col="blue") +
scale_y_continuous(limits = c(0,6)) + scale_x_continuous(limits = c(-0.3,0.6)); plot(p.vul.C)
p.correlation.pannels <- ggarrange(plotlist = list(p.cor.B,p.cor.C,p.cor.A,p.vul.B,p.vul.C,p.vul.A),
nrow = 2, ncol=3,
common.legend = T, legend = "right", labels = c("A","B","C","D","E","F")
); plot(p.correlation.pannels)
ggsave(plot = p.correlation.pannels, filename = "./Figures/Figure.S11.Ethanol.13fam.correlations.pdf", device = "pdf", width = 14, height = 7)
ggsave(plot = p.correlation.pannels, filename = "./Figures/Figure.S11.Ethanol.13fam.correlations.png", device = "png", width = 14, height = 7)Figure 2
Baria_metagenomics <- rio::import_list("./Tables/Supplementary.File.3.BARIA.Metagenomics.xlsx")
df.deseq.sub <- Baria_metagenomics$DESeq_NAFLD[Baria_metagenomics$DESeq_NAFLD$Family %in% top13.fams,]
p.deseq <- ggplot(df.deseq.sub, aes(x=-log2FoldChange, y=-log(padj), color=Family)) +
geom_point(size=2) +
facet_wrap(~Phylum, scales="free_x") +
scale_color_manual(values=colvec.fam) +
labs(y="-log(p.adjusted)") +
geom_hline(yintercept = -log(0.01)) +
#theme(legend.spacing.x = unit(0.20, 'cm')) +
theme(legend.key.size = unit(0.2, "cm")) +
NULL; p.deseq
select_pp=c("Lactobacillus fermentum",
"Lactobacillus casei",
"Lactobacillus reuteri",
"Lactobacillus salivarius",
"Lactobacillus johnsonii",
"Streptococcus thermophilus")
p.cortop <- rho_abundance_plot(df.cor.postp, select=select_pp, rho_p=rho_p_postp) +
scale_x_continuous(limits = c(-0.3,1.1)) +
NULL; plot(p.cortop)
df <- ps.BARIA.fam@sam_data
df$LAB <- t(ps.BARIA.fam@otu_table[ps.BARIA.fam@tax_table[,5] %in% "Lactobacillaceae",]@.Data)
p.cor2 <- ggplot(df, aes(x=LAB, y=Post.prandial, color=NAFLD, shape=NAFLD)) +
geom_point(size=1.5) +
scale_x_log10() +
scale_y_continuous(limits = c(0,230)) +
scale_color_manual(values=colvec[1:3]) +
labs(x=substitute(paste(italic(Lactobacillaceae)," abundance")), y = "Post prandial ethanol (mM)") +
stat_cor(aes(group=NA), method = "s", label.x = -5, label.y = 230, show.legend = F) +
#facet_wrap(~NAFLD) +
geom_smooth(method = "lm", aes(group=NA), show.legend = F) +
NULL; p.cor2
p.cor1 <- ggplot(df, aes(x=LAB, y=Fasted, color=NAFLD, shape=NAFLD)) +
geom_point(size=1.5) +
scale_x_log10() +
scale_y_continuous(limits = c(0,230)) +
scale_color_manual(values=colvec[1:3]) +
labs(x=substitute(paste(italic(Lactobacillaceae)," abundance")), y = "Fasted ethanol (mM)") +
stat_cor(aes(group=NA), method = "s", label.x = -5, label.y = 230, show.legend = F) +
#facet_wrap(~NAFLD) +
geom_smooth(method = "lm", aes(group=NA), show.legend = F) +
NULL; p.cor1
p.spec <- ggarrange(plotlist = list(p.deseq, p.cortop), labels = c("A","B"),
common.legend = T,
legend = "right"); plot(p.spec)
p.cor <- ggarrange(p.cor1, p.cor2, labels=c("C","D"),
common.legend = T, legend = "right")
p.metagenome <- ggarrange(p.spec, p.cor, ncol=1); p.metagenome 
ggsave(plot = p.metagenome, filename = "./Figures/Figure.2.metagenome.panels.png", height = 6, width = 9, device = "png")
ggsave(plot = p.metagenome, filename = "./Figures/Figure.2.metagenome.panels.pdf", height = 6, width = 9, device = "pdf")test medication
ps.BARIA <- readRDS("./Data/BARIA.metagenomics.RDS")
meds <- rio::import("Data/medicatie.xlsx")
meds[is.na(meds)] <- 0
meds <- do.call(cbind,lapply(meds, as.numeric)) %>% as.data.frame()
rownames(meds) <- paste0("SB_",meds$ID)
sample_data(ps.BARIA) <- sample_data(cbind(ps.BARIA@sam_data, meds[ps.BARIA@sam_data$Subject_ID,]))
ps.BARIA.compositional <- transform_sample_counts(ps.BARIA, function(x) x/sum(x))
ps.BARIA.fam <- microbiome::aggregate_taxa(ps.BARIA.compositional, "Family")
df <- ps.BARIA.fam@sam_data
df$LAB <- t(unclass((ps.BARIA.fam@otu_table[ps.BARIA.fam@tax_table[,5] %in% "Lactobacillaceae",]@.Data)))[,1]
df2 <- df[!is.na(df$Post.prandial),]
summary(lm(Post.prandial~ PPI+Metformine+Statines+LAB, data = data.frame(df)))##
## Call:
## lm(formula = Post.prandial ~ PPI + Metformine + Statines + LAB,
## data = data.frame(df))
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.172 -34.115 -1.039 20.646 143.076
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.890 5.196 14.605 < 2e-16 ***
## PPI 20.551 10.806 1.902 0.060002 .
## Metformine 18.437 15.538 1.187 0.238132
## Statines -17.086 16.215 -1.054 0.294502
## LAB 534.034 143.600 3.719 0.000326 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 44.53 on 103 degrees of freedom
## (14 observations deleted due to missingness)
## Multiple R-squared: 0.2036, Adjusted R-squared: 0.1726
## F-statistic: 6.581 on 4 and 103 DF, p-value: 9.32e-05summary(lm(Portal~ PPI+Metformine+Statines+LAB, data = data.frame(df)))##
## Call:
## lm(formula = Portal ~ PPI + Metformine + Statines + LAB, data = data.frame(df))
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.686 -13.349 -2.446 12.787 47.106
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.250 4.907 2.700 0.012 *
## PPI 8.101 8.425 0.962 0.345
## Metformine 3.397 12.006 0.283 0.779
## Statines -14.438 12.881 -1.121 0.273
## LAB 3024.708 2711.146 1.116 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.63 on 26 degrees of freedom
## (91 observations deleted due to missingness)
## Multiple R-squared: 0.1456, Adjusted R-squared: 0.01418
## F-statistic: 1.108 on 4 and 26 DF, p-value: 0.374summary(lm(Fasted~ PPI+Metformine+Statines+LAB, data = data.frame(df)))##
## Call:
## lm(formula = Fasted ~ PPI + Metformine + Statines + LAB, data = data.frame(df))
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.10 -18.22 -4.13 19.30 63.37
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.632 2.678 17.785 <2e-16 ***
## PPI 12.092 5.530 2.187 0.0308 *
## Metformine 7.940 7.938 1.000 0.3192
## Statines -11.875 8.468 -1.402 0.1635
## LAB 61.496 77.663 0.792 0.4301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.29 on 117 degrees of freedom
## Multiple R-squared: 0.06279, Adjusted R-squared: 0.03075
## F-statistic: 1.96 on 4 and 117 DF, p-value: 0.1052ggplot(df2, aes(x=factor(PPI), y=Post.prandial)) +
gghalves::geom_half_violin() + gghalves::geom_half_point() +
stat_compare_means()
ggplot(df2, aes(x=factor(PPI), y=LAB)) +
gghalves::geom_half_violin() + gghalves::geom_half_point() +
scale_y_log10() +
stat_compare_means()
df2$NAFLD2 <- as.numeric((factor(df2$NAFLD)))
summary(aov(NAFLD2~Post.prandial + PPI, data.frame(df2)))## Df Sum Sq Mean Sq F value Pr(>F)
## Post.prandial 1 7.29 7.292 8.591 0.00415 **
## PPI 1 0.50 0.498 0.587 0.44543
## Residuals 105 89.13 0.849
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df$NAFLD3 <- df$NAFLD=="No Steatosis"
chisq.test(df$PPI, df$NAFLD3)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: df$PPI and df$NAFLD3
## X-squared = 1.506, df = 1, p-value = 0.2198test function
These analysis are to test whether functional profiles or specific functions of the microbiome relate to plasma ethanol levels or NAFLD.
Correlate ethanol
correlations_wrapper <- function(ps, var, method="spearman", taxa=NULL){
ps.cor <- ps
ps.cor <- prune_samples(!is.na(get(var, ps@sam_data)), ps.cor)
ps.cor <- prune_taxa(taxa_sums(ps.cor)!=0, ps.cor)
if (is.null(taxa)){taxa=taxa_names(ps.cor)} else {ps.cor <- prune_taxa(taxa, ps.cor)}
corlist <- list()
for(taxa in taxa){
corlist[[taxa]] <- cor.test(ps.cor@otu_table[,taxa]@.Data, get(var, ps.cor@sam_data), method = method)
}
df.cor <- data.frame(taxa=taxa_names(ps.cor), row.names=taxa_names(ps.cor))
df.cor$p.value <- unlist(lapply(corlist, function(x) x$p.value))[rownames(df.cor)]
df.cor$estimate <- unlist(lapply(corlist, function(x) unname(x$estimate)))[rownames(df.cor)]
df.cor$Mean_Abundance <- (taxa_sums(ps)/sum(ps.cor@otu_table))[rownames(df.cor)]
df.cor$padj <- p.adjust(df.cor$p.value, method = "fdr", n=ntaxa(ps))
return(df.cor)
}
corres.pp <- correlations_wrapper(ps = ps.BARIA.ko, var = "Post_Prandial", method = "spearman", taxa=NULL)
corres.fa <- correlations_wrapper(ps = ps.BARIA.ko, var = "Fasted")
corres.po <- correlations_wrapper(ps = ps.BARIA.ko, var = "Portal")KO correlation tables
top
head(corres.fa %>% arrange(padj))## taxa p.value estimate Mean_Abundance padj
## K03969 K03969 2.056853e-05 -0.3781772 9.818170e-06 0.1491835
## K02110 K02110 4.103177e-04 -0.3174863 4.197645e-05 0.4548749
## K01715 K01715 4.390079e-04 -0.3159607 6.296877e-05 0.4548749
## K01286 K01286 2.329158e-04 -0.3299620 1.966092e-05 0.4548749
## K00248 K00248 2.706677e-04 -0.3267042 1.118136e-04 0.4548749
## K06381 K06381 3.937042e-04 -0.3184153 2.008878e-04 0.4548749head(corres.pp %>% arrange(padj))## taxa p.value estimate Mean_Abundance padj
## K05823 K05823 3.822188e-05 0.3887333 9.331590e-06 0.07964749
## K02859 K02859 5.090899e-05 0.3829396 7.027598e-07 0.07964749
## K02244 K02244 1.397202e-04 0.3616174 1.054280e-05 0.07964749
## K11704 K11704 1.652919e-04 0.3579188 1.264205e-06 0.07964749
## K08987 K08987 7.978974e-05 0.3736320 2.453135e-06 0.07964749
## K03471 K03471 1.239573e-04 0.3642245 1.634012e-05 0.07964749head(corres.po %>% arrange(padj))## taxa p.value estimate Mean_Abundance padj
## K01752 K01752 7.416231e-05 -0.65931045 1.267872e-03 0.5378993
## K01361 K01361 7.511056e-01 0.06042378 3.603966e-04 1.0000000
## K01362 K01362 8.149968e-01 0.04459320 1.901621e-04 1.0000000
## K05849 K05849 5.075132e-01 -0.12586003 2.873581e-07 1.0000000
## K05841 K05841 6.309775e-01 -0.09139954 4.741408e-07 1.0000000
## K12762 K12762 1.986458e-01 -0.24145357 8.031658e-07 1.0000000None of the kegg orthologs significantly correlate with the plasma ethanol levels.
Fasted
DT::datatable(corres.fa)Post Prandial
DT::datatable(corres.fa)Portal
DT::datatable(corres.fa)DESeq NAFLD
deseq.ko <- phyloseq_to_deseq2(ps.BARIA.ko, ~ NAFLD)
deseq.ko <- DESeq(deseq.ko, test="Wald", fitType="parametric")resultsNames(deseq.ko)## [1] "Intercept" "NAFLD_NAFL_vs_No.Steatosis"
## [3] "NAFLD_NASH_vs_No.Steatosis"resultsNames(deseq.ko)## [1] "Intercept" "NAFLD_NAFL_vs_No.Steatosis"
## [3] "NAFLD_NASH_vs_No.Steatosis"res_NAFL = results(deseq.ko, cooksCutoff = FALSE, name = "NAFLD_NAFL_vs_No.Steatosis")
res_NASH = results(deseq.ko, cooksCutoff = FALSE, name = "NAFLD_NASH_vs_No.Steatosis")res_NAFL %>% data.frame() %>% arrange(padj) %>% head## baseMean log2FoldChange lfcSE stat pvalue padj
## K14194 37.800407 -2.197480 0.3969061 -5.536522 3.085363e-08 0.0001711451
## K16199 4.034747 2.085127 0.4222269 4.938403 7.876479e-07 0.0017094031
## K11777 52.503635 2.422286 0.4936325 4.907063 9.245014e-07 0.0017094031
## K03030 24.826365 -3.309159 0.6861511 -4.822785 1.415674e-06 0.0019631852
## K16013 470.632288 1.487739 0.3264136 4.557833 5.168400e-06 0.0057338227
## K10681 4.260890 2.572521 0.6300713 4.082904 4.447637e-05 0.0411184016res_NASH %>% data.frame() %>% arrange(padj) %>% head## baseMean log2FoldChange lfcSE stat pvalue padj
## K01157 14.98020 -20.325274 1.9597553 -10.371332 3.348236e-25 1.683828e-21
## K09136 506.89465 -2.502615 0.5966891 -4.194169 2.738737e-05 5.115310e-02
## K14194 37.80041 -2.900380 0.6956031 -4.169590 3.051487e-05 5.115310e-02
## K07348 16.60025 -4.545975 1.2599257 -3.608129 3.084131e-04 9.972008e-02
## K04765 77.47199 -1.662725 0.4545739 -3.657767 2.544223e-04 9.972008e-02
## K01184 63.32333 -4.959110 1.3728826 -3.612188 3.036243e-04 9.972008e-02ko.sighits <- res_NAFL %>% data.frame() %>% subset(padj<0.05)While there are a few significant hits these seem to be false positive hits with low prevalence. They are not annotated with functions likely to be involed in NAFL development.
Export Kegg ortholog - ethanol NAFLD correlations
cor.res.list <- list(corres.fa,corres.pp, corres.po, res_NAFL, res_NASH)
names(cor.res.list) <- c("Correlations_KO_Fasted","Correlations_KO_Postprandial","Correlations_KO_Portal","DESeq_KO_NAFL","DESeq_KO_NASH")
rio::export(x = cor.res.list, file = "Tables/Supplementary.Table.BARIA.Metagenomics.KO.xlsx")Alcohol specifc KO
Due to multiple testing issues; were testing 7253 KEGG orthologs; we might be missing something. Also checking these specifically
alcoholDH1111 <- c("K00001","K00121","K04072","K11440","K13951","K13952","K13953","K13954","K13980","K18857")
acetaldehydeDH12110 <- c("K00132","K04072","K04073","K18366")
KO_select <- c(alcoholDH1111[alcoholDH1111 %in% taxa_names(ps.BARIA.ko)],acetaldehydeDH12110[acetaldehydeDH12110 %in% taxa_names(ps.BARIA.ko)])res_NAFL[KO_select,]## log2 fold change (MLE): NAFLD NAFL vs No.Steatosis
## Wald test p-value: NAFLD NAFL vs No.Steatosis
## DataFrame with 10 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## K00001 1.79766e+03 -0.0132220 0.1254503 -0.105397 0.9160611 0.998325
## K00121 6.28905e+01 0.1594518 0.3363747 0.474030 0.6354783 0.998325
## K04072 1.55072e+04 0.0900515 0.0788318 1.142325 0.2533191 0.998325
## K11440 3.96785e+00 -0.6448604 0.5341738 -1.207211 0.2273510 0.998325
## K13951 1.57895e-01 -0.0930079 3.0903779 -0.030096 0.9759905 NA
## K13953 2.59395e+02 0.5550805 0.2706014 2.051285 0.0402392 0.855199
## K13954 6.50818e+02 0.3243162 0.1543689 2.100917 0.0356483 0.800571
## K00132 1.00197e+02 0.0671477 0.2102646 0.319348 0.7494623 0.998325
## K04072 1.55072e+04 0.0900515 0.0788318 1.142325 0.2533191 0.998325
## K04073 1.17462e+01 -0.6056677 0.4925668 -1.229615 0.2188412 0.998325res_NASH[KO_select,]## log2 fold change (MLE): NAFLD NASH vs No.Steatosis
## Wald test p-value: NAFLD NASH vs No.Steatosis
## DataFrame with 10 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## K00001 1.79766e+03 0.2000848 0.217077 0.9217223 0.356673 0.751389
## K00121 6.28905e+01 -0.4267302 0.583276 -0.7316096 0.464407 0.779543
## K04072 1.55072e+04 0.1433680 0.136425 1.0508941 0.293307 0.723060
## K11440 3.96785e+00 -1.4765800 0.953902 -1.5479362 0.121638 0.602082
## K13951 1.57895e-01 -0.5187948 5.348146 -0.0970046 0.922723 NA
## K13953 2.59395e+02 -0.0930737 0.468552 -0.1986409 0.842544 0.949396
## K13954 6.50818e+02 -0.1365042 0.267291 -0.5106944 0.609565 0.844645
## K00132 1.00197e+02 0.4688185 0.363237 1.2906676 0.196819 0.661633
## K04072 1.55072e+04 0.1433680 0.136425 1.0508941 0.293307 0.723060
## K04073 1.17462e+01 -0.9648938 0.857684 -1.1249988 0.260590 0.703926corres.fa[KO_select,]## taxa p.value estimate Mean_Abundance padj
## K00001 K00001 0.44635344 0.07016456 3.870717e-05 1.0000000
## K00121 K00121 0.70909530 0.03440665 1.345116e-06 1.0000000
## K04072 K04072 0.43728140 0.07156921 3.380731e-04 1.0000000
## K11440 K11440 0.09530871 -0.15297367 8.406974e-08 0.9677093
## K13951 K13951 0.20975879 0.11531955 3.384163e-09 1.0000000
## K13953 K13953 0.78920676 -0.02465916 5.482879e-06 1.0000000
## K13954 K13954 0.53659480 0.05696312 1.392263e-05 1.0000000
## K00132 K00132 0.33921385 0.08799654 2.184923e-06 1.0000000
## K04072.1 K04072 0.43728140 0.07156921 3.380731e-04 1.0000000
## K04073 K04073 0.22043033 0.11268460 2.530998e-07 1.0000000corres.pp[KO_select,]## taxa p.value estimate Mean_Abundance padj
## K00001 K00001 0.36215611 0.08939237 4.362229e-05 0.8459640
## K00121 K00121 0.20740017 0.12344926 1.515922e-06 0.7632032
## K04072 K04072 0.09171381 0.16463250 3.810024e-04 0.6685430
## K11440 K11440 0.05844146 -0.18441024 9.474511e-08 0.5704925
## K13951 K13951 0.35501383 0.09072527 3.813892e-09 0.8420468
## K13953 K13953 0.04431755 0.19575973 6.179107e-06 0.4840892
## K13954 K13954 0.88972024 0.01362865 1.569055e-05 1.0000000
## K00132 K00132 0.49670465 0.06673435 2.462369e-06 0.9025050
## K04072.1 K04072 0.09171381 0.16463250 3.810024e-04 0.6685430
## K04073 K04073 0.29126855 0.10345528 2.852390e-07 0.8234219corres.po[KO_select,]## taxa p.value estimate Mean_Abundance padj
## K00001 K00001 0.8976371 0.02452626 1.561195e-04 1
## K00121 K00121 0.1269675 0.28493416 5.425321e-06 1
## K04072 K04072 0.3066332 -0.19308855 1.363567e-03 1
## K11440 K11440 0.4975153 -0.12881510 3.390825e-07 1
## K13951 K13951 0.1137102 0.29485026 1.364951e-08 1
## K13953 K13953 0.3921125 -0.16209628 2.211436e-05 1
## K13954 K13954 0.2524551 -0.21563210 5.615480e-05 1
## K00132 K00132 0.8372678 0.03914795 8.812554e-06 1
## K04072.1 K04072 0.3066332 -0.19308855 1.363567e-03 1
## K04073 K04073 0.7400754 0.06319427 1.020840e-06 1Some associations have p.values just below 0.05 but with the number of test, these are regarded as likely false positives.
sessionInfo
sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_U.S. Virgin Islands.1252
## [2] LC_CTYPE=English_U.S. Virgin Islands.1252
## [3] LC_MONETARY=English_U.S. Virgin Islands.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_U.S. Virgin Islands.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggpmisc_0.3.9 gghalves_0.1.1
## [3] phyloseq_1.34.0 enrichR_3.0
## [5] biomaRt_2.46.3 rio_0.5.26
## [7] reshape2_1.4.4 RColorBrewer_1.1-2
## [9] DT_0.18 DESeq2_1.30.1
## [11] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [13] MatrixGenerics_1.2.1 matrixStats_0.58.0
## [15] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [17] IRanges_2.24.1 S4Vectors_0.28.1
## [19] BiocGenerics_0.36.1 vegan_2.5-7
## [21] lattice_0.20-41 permute_0.9-5
## [23] lmerTest_3.1-3 lme4_1.1-26
## [25] Matrix_1.3-2 forcats_0.5.1
## [27] stringr_1.4.0 dplyr_1.0.4
## [29] purrr_0.3.4 readr_1.4.0
## [31] tidyr_1.1.3 tibble_3.0.4
## [33] tidyverse_1.3.1 ggpubr_0.4.0
## [35] ggsci_2.9 ggplot2_3.3.5
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 BiocFileCache_1.14.0
## [4] igraph_1.2.6 plyr_1.8.6 splines_4.0.5
## [7] crosstalk_1.1.1 BiocParallel_1.24.1 digest_0.6.27
## [10] foreach_1.5.1 htmltools_0.5.1.1 fansi_0.4.2
## [13] magrittr_2.0.1 memoise_2.0.0 cluster_2.1.1
## [16] openxlsx_4.2.3 Biostrings_2.58.0 annotate_1.68.0
## [19] modelr_0.1.8 askpass_1.1 prettyunits_1.1.1
## [22] colorspace_2.0-0 ggrepel_0.9.0 blob_1.2.1
## [25] rvest_1.0.0 rappdirs_0.3.3 haven_2.3.1
## [28] xfun_0.29 crayon_1.4.1 RCurl_1.98-1.3
## [31] jsonlite_1.7.2 genefilter_1.72.1 ape_5.5
## [34] iterators_1.0.13 survival_3.2-10 glue_1.4.2
## [37] gtable_0.3.0 zlibbioc_1.36.0 XVector_0.30.0
## [40] DelayedArray_0.16.3 car_3.0-10 Rhdf5lib_1.12.1
## [43] abind_1.4-5 scales_1.1.1 DBI_1.1.1
## [46] rstatix_0.7.0 Rcpp_1.0.7 viridisLite_0.4.0
## [49] xtable_1.8-4 progress_1.2.2 foreign_0.8-81
## [52] bit_4.0.4 htmlwidgets_1.5.3 httr_1.4.2
## [55] ellipsis_0.3.1 farver_2.1.0 pkgconfig_2.0.3
## [58] XML_3.99-0.6 sass_0.4.0 dbplyr_2.1.1
## [61] locfit_1.5-9.4 utf8_1.2.1 labeling_0.4.2
## [64] tidyselect_1.1.1 rlang_0.4.10 AnnotationDbi_1.52.0
## [67] munsell_0.5.0 cellranger_1.1.0 tools_4.0.5
## [70] cachem_1.0.4 cli_2.5.0 generics_0.1.0
## [73] RSQLite_2.2.7 ade4_1.7-16 broom_0.7.9
## [76] biomformat_1.18.0 evaluate_0.14 fastmap_1.1.0
## [79] yaml_2.2.1 knitr_1.33 bit64_4.0.5
## [82] fs_1.5.0 zip_2.1.1 nlme_3.1-152
## [85] xml2_1.3.2 compiler_4.0.5 rstudioapi_0.13
## [88] beeswarm_0.3.1 curl_4.3 ggsignif_0.6.1
## [91] reprex_2.0.0 statmod_1.4.36 geneplotter_1.68.0
## [94] bslib_0.2.5.1 stringi_1.5.3 highr_0.9
## [97] nloptr_1.2.2.2 microbiome_1.12.0 multtest_2.46.0
## [100] vctrs_0.3.6 rhdf5filters_1.2.0 pillar_1.6.0
## [103] lifecycle_1.0.0 jquerylib_0.1.4 cowplot_1.1.1
## [106] data.table_1.14.0 bitops_1.0-7 R6_2.5.0
## [109] gridExtra_2.3 vipor_0.4.5 codetools_0.2-18
## [112] boot_1.3-27 MASS_7.3-53.1 assertthat_0.2.1
## [115] rhdf5_2.34.0 rjson_0.2.20 openssl_1.4.4
## [118] withr_2.4.2 GenomeInfoDbData_1.2.4 mgcv_1.8-34
## [121] hms_1.0.0 grid_4.0.5 minqa_1.2.4
## [124] rmarkdown_2.11 carData_3.0-4 Rtsne_0.15
## [127] numDeriv_2016.8-1.1 lubridate_1.7.10 ggbeeswarm_0.6.0