This analysis pipe noted https://github.com/minsiksudo/BTE3207_Advanced_Biostatistics.
Roading the data
from previous lecture.. we generated a phyloseq object. it can be saved to hard disk using the below command.
saveRDS(ps, "phyloseq_example.rds")
And it can be loaded to a new R session like below!
getwd()## [1] "/Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics"library(phyloseq)
library(tidyverse)
library(microbiome)
library(vegan)
library(ggplot2)
#source("https://raw.githubusercontent.com/joey711/phyloseq/master/inst/scripts/installer.R",
# local = TRUE)
phyloseq <- readRDS("/Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/phyloseq_example.rds")
phyloseq_rel <- transform_sample_counts(phyloseq, function(x){x/sum(x)})Recall
From the barplot…
# Create a df with all `Top 10 Genus` and flag the top 20
plot_bar(phyloseq_rel, fill = "Genus") +
ylab("Total read counts") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90),
legend.position = "none") +
#scale_fill_brewer(type = "qual", palette = 3) +
facet_wrap(~When, scales = "free") +
ggtitle("Microbiome composition by sampling event time")
cowplot::get_legend(
plot_bar(phyloseq_rel, fill = "Genus") +
ylab("Total read counts") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90)) +
labs(fill = "ASVs and their Genus") +
#scale_fill_brewer(type = "qual", palette = 3) +
facet_wrap(~When, scales = "free") +
ggtitle("Microbiome composition by sampling event time")
) %>%
ggpubr::as_ggplot()
Can you see the difference?
Differential abundance analysis using your model
Using linear models 1
#Loading OTU table
otu_table <- phyloseq_rel %>% otu_table %>% t
#tax table for different sample type
taxa_data <- otu_table[1] %>% t %>% data.frame()
#Making a merged dataframe having sample data and CLR transformed output
lme_data <- merge(taxa_data, sample_data(phyloseq_rel), by = 0) %>% column_to_rownames("Row.names")
ggplot(lme_data, aes(x = Day, y = ASV1)) +
geom_point()
It looks like ther is some difference (in ther relative abundacnes)
Can we test them?
lm(ASV1 ~ When, data = lme_data) %>%
summary()##
## Call:
## lm(formula = ASV1 ~ When, data = lme_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.061144 -0.029114 -0.005525 0.020382 0.085583
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.12328 0.01258 9.799 2.08e-08 ***
## WhenLate -0.02523 0.01734 -1.455 0.164
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03774 on 17 degrees of freedom
## Multiple R-squared: 0.1107, Adjusted R-squared: 0.05842
## F-statistic: 2.117 on 1 and 17 DF, p-value: 0.1639THe difference is not significant
Using linear models 2
How does it look like with the second ASV
#Loading OTU table
otu_table <- phyloseq_rel %>% otu_table %>% t
#tax table for different sample type
taxa_data <- otu_table[2] %>% t %>% data.frame()
#Making a merged dataframe having sample data and CLR transformed output
lme_data <- merge(taxa_data, sample_data(phyloseq_rel), by = 0) %>% column_to_rownames("Row.names")
ggplot(lme_data, aes(x = Day, y = ASV2)) +
geom_point()
How about statistics?
lm(ASV2 ~ When, data = lme_data) %>%
summary()##
## Call:
## lm(formula = ASV2 ~ When, data = lme_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.032858 -0.013527 -0.004422 0.011586 0.042471
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.086563 0.006825 12.682 4.3e-10 ***
## WhenLate -0.007591 0.009408 -0.807 0.431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02048 on 17 degrees of freedom
## Multiple R-squared: 0.03688, Adjusted R-squared: -0.01977
## F-statistic: 0.651 on 1 and 17 DF, p-value: 0.4309No difference…
Can we run all the statistical test at once?
Introduction to for loop
for(i in 1:10) {
print(i)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10x = 0
for(i in 1:10) {
x <- x + i
print(x + i)
}## [1] 2
## [1] 5
## [1] 9
## [1] 14
## [1] 20
## [1] 27
## [1] 35
## [1] 44
## [1] 54
## [1] 65For-loop for analysis result
THis is the example code for running lm test for the 1st ASV.
#with all samples
lm_taxa_all <- data.frame()
for(i in 1:1) {
#Creating a data frame that includes CLR transformed data of i-th bug.
#making differnt otu tables for each sample type
otu_table <- phyloseq_rel %>% otu_table
#tax table for different sample type
taxa_data <- otu_table[,i] %>% data.frame()
#Making a merged dataframe having sample data and CLR transformed output
lm_data <- merge(taxa_data, sample_data(phyloseq_rel), by = 0) %>% column_to_rownames("Row.names")
#generating a character of formula.
#Here, as the taxa are already CLR transformed, I did not make model at log-scale.
lmformula <- paste(names(taxa_data), "~", "When")
#BAL stratified analysis
all_result <- lm(formula = lmformula,
data = lm_data) %>%
summary %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(feature = names(taxa_data)) %>%
rownames_to_column("value") %>%
.[2,] %>%
remove_rownames()
#row binding all the associations of i-th taxa to one data frame
lm_taxa_all <- rbind(lm_taxa_all,
all_result) %>%
remove_rownames()
}
head(lm_taxa_all)## value Estimate Std. Error t value Pr(>|t|) feature
## 1 WhenLate -0.02523023 0.01734153 -1.454903 0.1639126 ASV1For-loop for DA analysis
By changing the number of for loop, we can run multiple test and aggregate the result into one dataframe.
#with all samples
lm_taxa_all <- data.frame()
all <- phyloseq_rel %>% taxa_sums() %>% length()
for(i in 1:all) {
#Creating a data frame that includes CLR transformed data of i-th bug.
#making differnt otu tables for each sample type
otu_table <- phyloseq_rel %>% otu_table
#tax table for different sample type
taxa_data <- otu_table[,i] %>% data.frame()
#Making a merged dataframe having sample data and CLR transformed output
lm_data <- merge(taxa_data, sample_data(phyloseq_rel), by = 0) %>% column_to_rownames("Row.names")
#generating a character of formula.
#Here, as the taxa are already CLR transformed, I did not make model at log-scale.
lmformula <- paste(names(taxa_data), "~", "When")
#BAL stratified analysis
all_result <- lm(formula = lmformula,
data = lm_data) %>%
summary %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(feature = names(taxa_data)) %>%
rownames_to_column("value") %>%
.[2,] %>%
remove_rownames()
#row binding all the associations of i-th taxa to one data frame
lm_taxa_all <- rbind(lm_taxa_all,
all_result) %>%
remove_rownames()
}
head(lm_taxa_all)## value Estimate Std. Error t value Pr(>|t|) feature
## 1 WhenLate -0.025230234 0.017341527 -1.4549027 1.639126e-01 ASV1
## 2 WhenLate -0.007591208 0.009408264 -0.8068660 4.308881e-01 ASV2
## 3 WhenLate -0.002902008 0.006347035 -0.4572227 6.532998e-01 ASV3
## 4 WhenLate 0.049923591 0.007639995 6.5345057 5.098529e-06 ASV4
## 5 WhenLate -0.021115907 0.013164545 -1.6039982 1.271301e-01 ASV5
## 6 WhenLate 0.067059346 0.009262459 7.2399073 1.382214e-06 ASV6Multple testing correction
Since we tested multiple tests at once, the p-value needs to be adjusted
lm_taxa_all <- lm_taxa_all %>%
subset(., !is.nan(.$`Pr(>|t|)`)) %>%
mutate(qval = qvalue::qvalue(`Pr(>|t|)`, lambda = 0)$qvalue,
p_bh = p.adjust(p = `Pr(>|t|)`, method = "BH"))
#Benjamini-Hochberg correction
head(lm_taxa_all)## value Estimate Std. Error t value Pr(>|t|) feature
## 1 WhenLate -0.025230234 0.017341527 -1.4549027 1.639126e-01 ASV1
## 2 WhenLate -0.007591208 0.009408264 -0.8068660 4.308881e-01 ASV2
## 3 WhenLate -0.002902008 0.006347035 -0.4572227 6.532998e-01 ASV3
## 4 WhenLate 0.049923591 0.007639995 6.5345057 5.098529e-06 ASV4
## 5 WhenLate -0.021115907 0.013164545 -1.6039982 1.271301e-01 ASV5
## 6 WhenLate 0.067059346 0.009262459 7.2399073 1.382214e-06 ASV6
## qval p_bh
## 1 0.3247614964 0.3247614964
## 2 0.4991709066 0.4991709066
## 3 0.6924978165 0.6924978165
## 4 0.0002161776 0.0002161776
## 5 0.2848574361 0.2848574361
## 6 0.0001465147 0.0001465147To visualize the output.
volcano plot
ggplot(lm_taxa_all, aes(x = Estimate, y = -log10(p_bh))) +
geom_point() +
geom_hline(yintercept = -log10(0.05), col = "red") +
annotate("text", x = 0.05, y = 1.5, label = "p-value = 0.05", col = "red")
It looks weird, as we only tested the relative abundance!
For-loop for DA analysis (transformed data)
We can make CLR transfomred data, to consider the caveat of relative abundance
phyloseq_clr <- microbiome::transform(phyloseq, transform = "clr")
#with all samples
lm_taxa_all <- data.frame()
all <- phyloseq_clr %>% taxa_sums() %>% length()
for(i in 1:all) {
#Creating a data frame that includes CLR transformed data of i-th bug.
#making differnt otu tables for each sample type
otu_table <- phyloseq_clr %>% otu_table
#tax table for different sample type
taxa_data <- otu_table[,i] %>% data.frame()
#Making a merged dataframe having sample data and CLR transformed output
lm_data <- merge(taxa_data, sample_data(phyloseq_clr), by = 0) %>% column_to_rownames("Row.names")
#generating a character of formula.
#Here, as the taxa are already CLR transformed, I did not make model at log-scale.
lmformula <- paste(names(taxa_data), "~", "When")
#BAL stratified analysis
all_result <- lm(formula = lmformula,
data = lm_data) %>%
summary %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(feature = names(taxa_data)) %>%
rownames_to_column("value") %>%
.[2,] %>%
remove_rownames()
#row binding all the associations of i-th taxa to one data frame
lm_taxa_all <- rbind(lm_taxa_all,
all_result) %>%
remove_rownames()
}
lm_taxa_all <- lm_taxa_all %>%
subset(., !is.nan(.$`Pr(>|t|)`)) %>%
mutate(qval = qvalue::qvalue(`Pr(>|t|)`, lambda = 0)$qvalue,
p_bh = p.adjust(p = `Pr(>|t|)`, method = "BH"))
#Benjamini-Hochberg correction
head(lm_taxa_all)## value Estimate Std. Error t value Pr(>|t|) feature qval
## 1 WhenLate 0.10211879 0.2343786 0.4357000 6.685362e-01 ASV1 0.8281112120
## 2 WhenLate 0.17875798 0.2171687 0.8231294 4.218316e-01 ASV2 0.6220427897
## 3 WhenLate 0.23437974 0.1929377 1.2147951 2.410485e-01 ASV3 0.4047456597
## 4 WhenLate 1.12402925 0.2423343 4.6383407 2.350677e-04 ASV4 0.0017630077
## 5 WhenLate 0.05008611 0.3177319 0.1576364 8.766005e-01 ASV5 0.9273201971
## 6 WhenLate 3.61718038 0.7028486 5.1464577 8.075004e-05 ASV6 0.0009562504
## p_bh
## 1 0.8281112120
## 2 0.6220427897
## 3 0.4047456597
## 4 0.0017630077
## 5 0.9273201971
## 6 0.0009562504volcano plot
ggplot(lm_taxa_all, aes(x = Estimate, y = -log10(p_bh))) +
geom_point() +
geom_hline(yintercept = -log10(0.05), col = "red") +
annotate("text", x = 0.05, y = 1.5, label = "p-value = 0.05", col = "red")
Volcano plot looks ok.
Using MaAslin
MaAsLin2 is comprehensive R package for efficiently determining multivariable association between phenotypes, environments, exposures, covariates and microbial meta’omic features. MaAsLin2 relies on general linear models to accommodate most modern epidemiological study designs, including cross-sectional and longitudinal, and offers a variety of data exploration, normalization, and transformation methods.
https://huttenhower.sph.harvard.edu/maaslin/
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Maaslin2")
library(Maaslin2)Running maaslin
capture.output(
maaslin_output<-
Maaslin2(input_data = otu_table(phyloseq_rel) %>%
data.frame(),
input_metadata = phyloseq_rel %>%
sample_data %>%
data.frame(check.names = F),
output = paste0(getwd(), "/maaslin_output") ,
fixed_effects = c("When"),
transform = "LOG", #default
normalization = "TSS", #Choose your own normalization method
reference = c("When,Early")))## [1] "[1] \"Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/maaslin2.log\""
## [2] "2023-11-15 15:37:38.831979 INFO::Writing function arguments to log file"
## [3] "2023-11-15 15:37:38.848363 INFO::Verifying options selected are valid"
## [4] "2023-11-15 15:37:38.884709 INFO::Determining format of input files"
## [5] "2023-11-15 15:37:38.885629 INFO::Input format is data samples as rows and metadata samples as rows"
## [6] "2023-11-15 15:37:38.907258 INFO::Formula for fixed effects: expr ~ When"
## [7] "2023-11-15 15:37:38.908122 INFO::Filter data based on min abundance and min prevalence"
## [8] "2023-11-15 15:37:38.908656 INFO::Total samples in data: 19"
## [9] "2023-11-15 15:37:38.909172 INFO::Min samples required with min abundance for a feature not to be filtered: 1.900000"
## [10] "2023-11-15 15:37:38.911586 INFO::Total filtered features: 48"
## [11] "2023-11-15 15:37:38.912209 INFO::Filtered feature names from abundance and prevalence filtering: ASV34, ASV42, ASV45, ASV51, ASV60, ASV70, ASV83, ASV91, ASV98, ASV102, ASV103, ASV107, ASV113, ASV139, ASV144, ASV154, ASV159, ASV169, ASV182, ASV186, ASV188, ASV192, ASV194, ASV195, ASV196, ASV198, ASV202, ASV204, ASV205, ASV207, ASV208, ASV209, ASV210, ASV211, ASV212, ASV213, ASV214, ASV215, ASV216, ASV217, ASV218, ASV219, ASV220, ASV221, ASV222, ASV223, ASV224, ASV225"
## [12] "2023-11-15 15:37:38.915025 INFO::Total filtered features with variance filtering: 0"
## [13] "2023-11-15 15:37:38.915694 INFO::Filtered feature names from variance filtering:"
## [14] "2023-11-15 15:37:38.916234 INFO::Running selected normalization method: TSS"
## [15] "2023-11-15 15:37:38.919081 INFO::Applying z-score to standardize continuous metadata"
## [16] "2023-11-15 15:37:38.924656 INFO::Running selected transform method: LOG"
## [17] "2023-11-15 15:37:38.927161 INFO::Running selected analysis method: LM"
## [18] "2023-11-15 15:37:38.93208 INFO::Fitting model to feature number 1, ASV1"
## [19] "2023-11-15 15:37:38.935748 INFO::Fitting model to feature number 2, ASV2"
## [20] "2023-11-15 15:37:38.937568 INFO::Fitting model to feature number 3, ASV3"
## [21] "2023-11-15 15:37:38.939429 INFO::Fitting model to feature number 4, ASV4"
## [22] "2023-11-15 15:37:38.941312 INFO::Fitting model to feature number 5, ASV5"
## [23] "2023-11-15 15:37:38.943175 INFO::Fitting model to feature number 6, ASV6"
## [24] "2023-11-15 15:37:38.945043 INFO::Fitting model to feature number 7, ASV7"
## [25] "2023-11-15 15:37:38.946921 INFO::Fitting model to feature number 8, ASV8"
## [26] "2023-11-15 15:37:38.949298 INFO::Fitting model to feature number 9, ASV9"
## [27] "2023-11-15 15:37:38.951275 INFO::Fitting model to feature number 10, ASV10"
## [28] "2023-11-15 15:37:38.953255 INFO::Fitting model to feature number 11, ASV11"
## [29] "2023-11-15 15:37:38.955386 INFO::Fitting model to feature number 12, ASV12"
## [30] "2023-11-15 15:37:38.957617 INFO::Fitting model to feature number 13, ASV13"
## [31] "2023-11-15 15:37:38.959716 INFO::Fitting model to feature number 14, ASV14"
## [32] "2023-11-15 15:37:38.961851 INFO::Fitting model to feature number 15, ASV15"
## [33] "2023-11-15 15:37:38.964076 INFO::Fitting model to feature number 16, ASV16"
## [34] "2023-11-15 15:37:38.966185 INFO::Fitting model to feature number 17, ASV17"
## [35] "2023-11-15 15:37:38.968248 INFO::Fitting model to feature number 18, ASV18"
## [36] "2023-11-15 15:37:38.970722 INFO::Fitting model to feature number 19, ASV19"
## [37] "2023-11-15 15:37:38.992317 INFO::Fitting model to feature number 20, ASV20"
## [38] "2023-11-15 15:37:38.994732 INFO::Fitting model to feature number 21, ASV21"
## [39] "2023-11-15 15:37:39.014759 INFO::Fitting model to feature number 22, ASV22"
## [40] "2023-11-15 15:37:39.018486 INFO::Fitting model to feature number 23, ASV23"
## [41] "2023-11-15 15:37:39.022334 INFO::Fitting model to feature number 24, ASV24"
## [42] "2023-11-15 15:37:39.026264 INFO::Fitting model to feature number 25, ASV25"
## [43] "2023-11-15 15:37:39.028624 INFO::Fitting model to feature number 26, ASV26"
## [44] "2023-11-15 15:37:39.030923 INFO::Fitting model to feature number 27, ASV27"
## [45] "2023-11-15 15:37:39.033305 INFO::Fitting model to feature number 28, ASV28"
## [46] "2023-11-15 15:37:39.035571 INFO::Fitting model to feature number 29, ASV29"
## [47] "2023-11-15 15:37:39.037887 INFO::Fitting model to feature number 30, ASV30"
## [48] "2023-11-15 15:37:39.040174 INFO::Fitting model to feature number 31, ASV31"
## [49] "2023-11-15 15:37:39.042432 INFO::Fitting model to feature number 32, ASV32"
## [50] "2023-11-15 15:37:39.044615 INFO::Fitting model to feature number 33, ASV33"
## [51] "2023-11-15 15:37:39.046826 INFO::Fitting model to feature number 34, ASV35"
## [52] "2023-11-15 15:37:39.049097 INFO::Fitting model to feature number 35, ASV36"
## [53] "2023-11-15 15:37:39.051319 INFO::Fitting model to feature number 36, ASV37"
## [54] "2023-11-15 15:37:39.053583 INFO::Fitting model to feature number 37, ASV38"
## [55] "2023-11-15 15:37:39.055972 INFO::Fitting model to feature number 38, ASV39"
## [56] "2023-11-15 15:37:39.058402 INFO::Fitting model to feature number 39, ASV40"
## [57] "2023-11-15 15:37:39.060747 INFO::Fitting model to feature number 40, ASV41"
## [58] "2023-11-15 15:37:39.063038 INFO::Fitting model to feature number 41, ASV43"
## [59] "2023-11-15 15:37:39.065307 INFO::Fitting model to feature number 42, ASV44"
## [60] "2023-11-15 15:37:39.067544 INFO::Fitting model to feature number 43, ASV46"
## [61] "2023-11-15 15:37:39.069765 INFO::Fitting model to feature number 44, ASV47"
## [62] "2023-11-15 15:37:39.071942 INFO::Fitting model to feature number 45, ASV48"
## [63] "2023-11-15 15:37:39.074331 INFO::Fitting model to feature number 46, ASV49"
## [64] "2023-11-15 15:37:39.076662 INFO::Fitting model to feature number 47, ASV50"
## [65] "2023-11-15 15:37:39.078624 INFO::Fitting model to feature number 48, ASV52"
## [66] "2023-11-15 15:37:39.080945 INFO::Fitting model to feature number 49, ASV53"
## [67] "2023-11-15 15:37:39.083139 INFO::Fitting model to feature number 50, ASV54"
## [68] "2023-11-15 15:37:39.085345 INFO::Fitting model to feature number 51, ASV55"
## [69] "2023-11-15 15:37:39.087499 INFO::Fitting model to feature number 52, ASV56"
## [70] "2023-11-15 15:37:39.089536 INFO::Fitting model to feature number 53, ASV57"
## [71] "2023-11-15 15:37:39.091774 INFO::Fitting model to feature number 54, ASV58"
## [72] "2023-11-15 15:37:39.094287 INFO::Fitting model to feature number 55, ASV59"
## [73] "2023-11-15 15:37:39.096634 INFO::Fitting model to feature number 56, ASV61"
## [74] "2023-11-15 15:37:39.098851 INFO::Fitting model to feature number 57, ASV62"
## [75] "2023-11-15 15:37:39.10085 INFO::Fitting model to feature number 58, ASV63"
## [76] "2023-11-15 15:37:39.102944 INFO::Fitting model to feature number 59, ASV64"
## [77] "2023-11-15 15:37:39.105613 INFO::Fitting model to feature number 60, ASV65"
## [78] "2023-11-15 15:37:39.108295 INFO::Fitting model to feature number 61, ASV66"
## [79] "2023-11-15 15:37:39.123815 INFO::Fitting model to feature number 62, ASV67"
## [80] "2023-11-15 15:37:39.126476 INFO::Fitting model to feature number 63, ASV68"
## [81] "2023-11-15 15:37:39.148329 INFO::Fitting model to feature number 64, ASV69"
## [82] "2023-11-15 15:37:39.15079 INFO::Fitting model to feature number 65, ASV71"
## [83] "2023-11-15 15:37:39.153092 INFO::Fitting model to feature number 66, ASV72"
## [84] "2023-11-15 15:37:39.155369 INFO::Fitting model to feature number 67, ASV73"
## [85] "2023-11-15 15:37:39.157624 INFO::Fitting model to feature number 68, ASV74"
## [86] "2023-11-15 15:37:39.159852 INFO::Fitting model to feature number 69, ASV75"
## [87] "2023-11-15 15:37:39.162042 INFO::Fitting model to feature number 70, ASV76"
## [88] "2023-11-15 15:37:39.164258 INFO::Fitting model to feature number 71, ASV77"
## [89] "2023-11-15 15:37:39.166489 INFO::Fitting model to feature number 72, ASV78"
## [90] "2023-11-15 15:37:39.168667 INFO::Fitting model to feature number 73, ASV79"
## [91] "2023-11-15 15:37:39.170749 INFO::Fitting model to feature number 74, ASV80"
## [92] "2023-11-15 15:37:39.172877 INFO::Fitting model to feature number 75, ASV81"
## [93] "2023-11-15 15:37:39.175255 INFO::Fitting model to feature number 76, ASV82"
## [94] "2023-11-15 15:37:39.177566 INFO::Fitting model to feature number 77, ASV84"
## [95] "2023-11-15 15:37:39.179877 INFO::Fitting model to feature number 78, ASV85"
## [96] "2023-11-15 15:37:39.182434 INFO::Fitting model to feature number 79, ASV86"
## [97] "2023-11-15 15:37:39.184526 INFO::Fitting model to feature number 80, ASV87"
## [98] "2023-11-15 15:37:39.186474 INFO::Fitting model to feature number 81, ASV88"
## [99] "2023-11-15 15:37:39.189395 INFO::Fitting model to feature number 82, ASV89"
## [100] "2023-11-15 15:37:39.191884 INFO::Fitting model to feature number 83, ASV90"
## [101] "2023-11-15 15:37:39.194874 INFO::Fitting model to feature number 84, ASV92"
## [102] "2023-11-15 15:37:39.197556 INFO::Fitting model to feature number 85, ASV93"
## [103] "2023-11-15 15:37:39.199802 INFO::Fitting model to feature number 86, ASV94"
## [104] "2023-11-15 15:37:39.202035 INFO::Fitting model to feature number 87, ASV95"
## [105] "2023-11-15 15:37:39.204382 INFO::Fitting model to feature number 88, ASV96"
## [106] "2023-11-15 15:37:39.206885 INFO::Fitting model to feature number 89, ASV97"
## [107] "2023-11-15 15:37:39.209156 INFO::Fitting model to feature number 90, ASV99"
## [108] "2023-11-15 15:37:39.211362 INFO::Fitting model to feature number 91, ASV100"
## [109] "2023-11-15 15:37:39.213694 INFO::Fitting model to feature number 92, ASV101"
## [110] "2023-11-15 15:37:39.215914 INFO::Fitting model to feature number 93, ASV104"
## [111] "2023-11-15 15:37:39.218141 INFO::Fitting model to feature number 94, ASV105"
## [112] "2023-11-15 15:37:39.220837 INFO::Fitting model to feature number 95, ASV106"
## [113] "2023-11-15 15:37:39.224133 INFO::Fitting model to feature number 96, ASV108"
## [114] "2023-11-15 15:37:39.226523 INFO::Fitting model to feature number 97, ASV109"
## [115] "2023-11-15 15:37:39.228819 INFO::Fitting model to feature number 98, ASV110"
## [116] "2023-11-15 15:37:39.231255 INFO::Fitting model to feature number 99, ASV111"
## [117] "2023-11-15 15:37:39.252094 INFO::Fitting model to feature number 100, ASV112"
## [118] "2023-11-15 15:37:39.25435 INFO::Fitting model to feature number 101, ASV114"
## [119] "2023-11-15 15:37:39.257073 INFO::Fitting model to feature number 102, ASV115"
## [120] "2023-11-15 15:37:39.259044 INFO::Fitting model to feature number 103, ASV116"
## [121] "2023-11-15 15:37:39.261439 INFO::Fitting model to feature number 104, ASV117"
## [122] "2023-11-15 15:37:39.263362 INFO::Fitting model to feature number 105, ASV118"
## [123] "2023-11-15 15:37:39.26527 INFO::Fitting model to feature number 106, ASV119"
## [124] "2023-11-15 15:37:39.26716 INFO::Fitting model to feature number 107, ASV120"
## [125] "2023-11-15 15:37:39.269045 INFO::Fitting model to feature number 108, ASV121"
## [126] "2023-11-15 15:37:39.271143 INFO::Fitting model to feature number 109, ASV122"
## [127] "2023-11-15 15:37:39.273346 INFO::Fitting model to feature number 110, ASV123"
## [128] "2023-11-15 15:37:39.275343 INFO::Fitting model to feature number 111, ASV124"
## [129] "2023-11-15 15:37:39.277298 INFO::Fitting model to feature number 112, ASV125"
## [130] "2023-11-15 15:37:39.279248 INFO::Fitting model to feature number 113, ASV126"
## [131] "2023-11-15 15:37:39.281199 INFO::Fitting model to feature number 114, ASV127"
## [132] "2023-11-15 15:37:39.283135 INFO::Fitting model to feature number 115, ASV128"
## [133] "2023-11-15 15:37:39.285277 INFO::Fitting model to feature number 116, ASV129"
## [134] "2023-11-15 15:37:39.287209 INFO::Fitting model to feature number 117, ASV130"
## [135] "2023-11-15 15:37:39.289656 INFO::Fitting model to feature number 118, ASV131"
## [136] "2023-11-15 15:37:39.291561 INFO::Fitting model to feature number 119, ASV132"
## [137] "2023-11-15 15:37:39.293502 INFO::Fitting model to feature number 120, ASV133"
## [138] "2023-11-15 15:37:39.295438 INFO::Fitting model to feature number 121, ASV134"
## [139] "2023-11-15 15:37:39.297406 INFO::Fitting model to feature number 122, ASV135"
## [140] "2023-11-15 15:37:39.299399 INFO::Fitting model to feature number 123, ASV136"
## [141] "2023-11-15 15:37:39.301669 INFO::Fitting model to feature number 124, ASV137"
## [142] "2023-11-15 15:37:39.303648 INFO::Fitting model to feature number 125, ASV138"
## [143] "2023-11-15 15:37:39.305603 INFO::Fitting model to feature number 126, ASV140"
## [144] "2023-11-15 15:37:39.307542 INFO::Fitting model to feature number 127, ASV141"
## [145] "2023-11-15 15:37:39.309473 INFO::Fitting model to feature number 128, ASV142"
## [146] "2023-11-15 15:37:39.311414 INFO::Fitting model to feature number 129, ASV143"
## [147] "2023-11-15 15:37:39.313416 INFO::Fitting model to feature number 130, ASV145"
## [148] "2023-11-15 15:37:39.31535 INFO::Fitting model to feature number 131, ASV146"
## [149] "2023-11-15 15:37:39.317284 INFO::Fitting model to feature number 132, ASV147"
## [150] "2023-11-15 15:37:39.319213 INFO::Fitting model to feature number 133, ASV148"
## [151] "2023-11-15 15:37:39.321695 INFO::Fitting model to feature number 134, ASV149"
## [152] "2023-11-15 15:37:39.325387 INFO::Fitting model to feature number 135, ASV150"
## [153] "2023-11-15 15:37:39.327669 INFO::Fitting model to feature number 136, ASV151"
## [154] "2023-11-15 15:37:39.347105 INFO::Fitting model to feature number 137, ASV152"
## [155] "2023-11-15 15:37:39.349903 INFO::Fitting model to feature number 138, ASV153"
## [156] "2023-11-15 15:37:39.352171 INFO::Fitting model to feature number 139, ASV155"
## [157] "2023-11-15 15:37:39.354345 INFO::Fitting model to feature number 140, ASV156"
## [158] "2023-11-15 15:37:39.356386 INFO::Fitting model to feature number 141, ASV157"
## [159] "2023-11-15 15:37:39.358519 INFO::Fitting model to feature number 142, ASV158"
## [160] "2023-11-15 15:37:39.360317 INFO::Fitting model to feature number 143, ASV160"
## [161] "2023-11-15 15:37:39.362108 INFO::Fitting model to feature number 144, ASV161"
## [162] "2023-11-15 15:37:39.363952 INFO::Fitting model to feature number 145, ASV162"
## [163] "2023-11-15 15:37:39.365756 INFO::Fitting model to feature number 146, ASV163"
## [164] "2023-11-15 15:37:39.367624 INFO::Fitting model to feature number 147, ASV164"
## [165] "2023-11-15 15:37:39.369586 INFO::Fitting model to feature number 148, ASV165"
## [166] "2023-11-15 15:37:39.371401 INFO::Fitting model to feature number 149, ASV166"
## [167] "2023-11-15 15:37:39.373206 INFO::Fitting model to feature number 150, ASV167"
## [168] "2023-11-15 15:37:39.375078 INFO::Fitting model to feature number 151, ASV168"
## [169] "2023-11-15 15:37:39.376869 INFO::Fitting model to feature number 152, ASV170"
## [170] "2023-11-15 15:37:39.378662 INFO::Fitting model to feature number 153, ASV171"
## [171] "2023-11-15 15:37:39.38058 INFO::Fitting model to feature number 154, ASV172"
## [172] "2023-11-15 15:37:39.382466 INFO::Fitting model to feature number 155, ASV173"
## [173] "2023-11-15 15:37:39.384253 INFO::Fitting model to feature number 156, ASV174"
## [174] "2023-11-15 15:37:39.386094 INFO::Fitting model to feature number 157, ASV175"
## [175] "2023-11-15 15:37:39.387888 INFO::Fitting model to feature number 158, ASV176"
## [176] "2023-11-15 15:37:39.390557 INFO::Fitting model to feature number 159, ASV177"
## [177] "2023-11-15 15:37:39.392451 INFO::Fitting model to feature number 160, ASV178"
## [178] "2023-11-15 15:37:39.394437 INFO::Fitting model to feature number 161, ASV179"
## [179] "2023-11-15 15:37:39.396469 INFO::Fitting model to feature number 162, ASV180"
## [180] "2023-11-15 15:37:39.398321 INFO::Fitting model to feature number 163, ASV181"
## [181] "2023-11-15 15:37:39.400116 INFO::Fitting model to feature number 164, ASV183"
## [182] "2023-11-15 15:37:39.401907 INFO::Fitting model to feature number 165, ASV184"
## [183] "2023-11-15 15:37:39.403691 INFO::Fitting model to feature number 166, ASV185"
## [184] "2023-11-15 15:37:39.406046 INFO::Fitting model to feature number 167, ASV187"
## [185] "2023-11-15 15:37:39.40812 INFO::Fitting model to feature number 168, ASV189"
## [186] "2023-11-15 15:37:39.409939 INFO::Fitting model to feature number 169, ASV190"
## [187] "2023-11-15 15:37:39.411739 INFO::Fitting model to feature number 170, ASV191"
## [188] "2023-11-15 15:37:39.413622 INFO::Fitting model to feature number 171, ASV193"
## [189] "2023-11-15 15:37:39.415631 INFO::Fitting model to feature number 172, ASV197"
## [190] "2023-11-15 15:37:39.417466 INFO::Fitting model to feature number 173, ASV199"
## [191] "2023-11-15 15:37:39.41929 INFO::Fitting model to feature number 174, ASV200"
## [192] "2023-11-15 15:37:39.421456 INFO::Fitting model to feature number 175, ASV201"
## [193] "2023-11-15 15:37:39.425652 INFO::Fitting model to feature number 176, ASV203"
## [194] "2023-11-15 15:37:39.428895 INFO::Fitting model to feature number 177, ASV206"
## [195] "2023-11-15 15:37:39.448862 INFO::Counting total values for each feature"
## [196] "2023-11-15 15:37:39.456794 INFO::Writing filtered data to file /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/features/filtered_data.tsv"
## [197] "2023-11-15 15:37:39.479832 INFO::Writing filtered, normalized data to file /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/features/filtered_data_norm.tsv"
## [198] "2023-11-15 15:37:39.502827 INFO::Writing filtered, normalized, transformed data to file /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/features/filtered_data_norm_transformed.tsv"
## [199] "2023-11-15 15:37:39.531336 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/fits/residuals.rds"
## [200] "2023-11-15 15:37:39.533454 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/fits/residuals.rds"
## [201] "2023-11-15 15:37:39.535765 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/fits/fitted.rds"
## [202] "2023-11-15 15:37:39.537362 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/fits/fitted.rds"
## [203] "2023-11-15 15:37:39.557493 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/all_results.tsv"
## [204] "2023-11-15 15:37:39.560099 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/significant_results.tsv"
## [205] "2023-11-15 15:37:39.562769 INFO::Writing heatmap of significant results to file: /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output/heatmap.pdf"
## [206] "[1] \"There is not enough metadata in the associations to create a heatmap plot. Please review the associations in text output file.\""
## [207] "2023-11-15 15:37:39.565925 INFO::Writing association plots (one for each significant association) to output folder: /Users/minsikkim/Dropbox (Personal)/Inha/5_Lectures/Advanced metagenomics/scripts/IBS7048_Advanced_metagenomics/maaslin_output"
## [208] "2023-11-15 15:37:39.569565 INFO::Plotting associations from most to least significant, grouped by metadata"
## [209] "2023-11-15 15:37:39.570347 INFO::Plotting data for metadata number 1, When"
## [210] "2023-11-15 15:37:39.591008 INFO::Creating boxplot for categorical data, When vs ASV41"
## [211] "2023-11-15 15:37:39.794174 INFO::Creating boxplot for categorical data, When vs ASV36"
## [212] "2023-11-15 15:37:39.904872 INFO::Creating boxplot for categorical data, When vs ASV50"
## [213] "2023-11-15 15:37:40.056341 INFO::Creating boxplot for categorical data, When vs ASV61"
## [214] "2023-11-15 15:37:40.190964 INFO::Creating boxplot for categorical data, When vs ASV44"
## [215] "2023-11-15 15:37:40.340704 INFO::Creating boxplot for categorical data, When vs ASV25"
## [216] "2023-11-15 15:37:40.47179 INFO::Creating boxplot for categorical data, When vs ASV49"
## [217] "2023-11-15 15:37:40.606972 INFO::Creating boxplot for categorical data, When vs ASV35"
## [218] "2023-11-15 15:37:40.729461 INFO::Creating boxplot for categorical data, When vs ASV12"
## [219] "2023-11-15 15:37:40.882979 INFO::Creating boxplot for categorical data, When vs ASV6"
## [220] "2023-11-15 15:37:41.006994 INFO::Creating boxplot for categorical data, When vs ASV54"
## [221] "2023-11-15 15:37:41.140253 INFO::Creating boxplot for categorical data, When vs ASV22"
## [222] "2023-11-15 15:37:41.265097 INFO::Creating boxplot for categorical data, When vs ASV26"
## [223] "2023-11-15 15:37:41.395684 INFO::Creating boxplot for categorical data, When vs ASV109"
## [224] "2023-11-15 15:37:41.525673 INFO::Creating boxplot for categorical data, When vs ASV7"
## [225] "2023-11-15 15:37:41.679827 INFO::Creating boxplot for categorical data, When vs ASV57"
## [226] "2023-11-15 15:37:41.80813 INFO::Creating boxplot for categorical data, When vs ASV16"
## [227] "2023-11-15 15:37:41.927941 INFO::Creating boxplot for categorical data, When vs ASV108"
## [228] "2023-11-15 15:37:42.046657 INFO::Creating boxplot for categorical data, When vs ASV106"
## [229] "2023-11-15 15:37:42.199469 INFO::Creating boxplot for categorical data, When vs ASV19"
## [230] "2023-11-15 15:37:42.321867 INFO::Creating boxplot for categorical data, When vs ASV15"
## [231] "2023-11-15 15:37:42.450373 INFO::Creating boxplot for categorical data, When vs ASV4"
## [232] "2023-11-15 15:37:42.570806 INFO::Creating boxplot for categorical data, When vs ASV11"
## [233] "2023-11-15 15:37:42.73767 INFO::Creating boxplot for categorical data, When vs ASV59"
## [234] "2023-11-15 15:37:42.857416 INFO::Creating boxplot for categorical data, When vs ASV52"
## [235] "2023-11-15 15:37:43.00013 INFO::Creating boxplot for categorical data, When vs ASV67"
## [236] "2023-11-15 15:37:43.170808 INFO::Creating boxplot for categorical data, When vs ASV62"
## [237] "2023-11-15 15:37:43.292815 INFO::Creating boxplot for categorical data, When vs ASV69"
## [238] "2023-11-15 15:37:43.421452 INFO::Creating boxplot for categorical data, When vs ASV117"
## [239] "2023-11-15 15:37:43.544716 INFO::Creating boxplot for categorical data, When vs ASV97"
## [240] "2023-11-15 15:37:43.715448 INFO::Creating boxplot for categorical data, When vs ASV29"
## [241] "2023-11-15 15:37:43.838724 INFO::Creating boxplot for categorical data, When vs ASV104"
## [242] "2023-11-15 15:37:43.992208 INFO::Creating boxplot for categorical data, When vs ASV76"
## [243] "2023-11-15 15:37:44.115074 INFO::Creating boxplot for categorical data, When vs ASV142"
## [244] "2023-11-15 15:37:44.277281 INFO::Creating boxplot for categorical data, When vs ASV125"
## [245] "2023-11-15 15:37:44.398846 INFO::Creating boxplot for categorical data, When vs ASV87"
## [246] "2023-11-15 15:37:44.551135 INFO::Creating boxplot for categorical data, When vs ASV100"
## [247] "2023-11-15 15:37:44.694018 INFO::Creating boxplot for categorical data, When vs ASV14"
## [248] "2023-11-15 15:37:44.849033 INFO::Creating boxplot for categorical data, When vs ASV71"
## [249] "2023-11-15 15:37:44.96733 INFO::Creating boxplot for categorical data, When vs ASV77"
## [250] "2023-11-15 15:37:45.10056 INFO::Creating boxplot for categorical data, When vs ASV32"
## [251] "2023-11-15 15:37:45.223059 INFO::Creating boxplot for categorical data, When vs ASV80"
## [252] "2023-11-15 15:37:45.356542 INFO::Creating boxplot for categorical data, When vs ASV110"
## [253] "2023-11-15 15:37:45.510903 INFO::Creating boxplot for categorical data, When vs ASV13"
## [254] "2023-11-15 15:37:45.642071 INFO::Creating boxplot for categorical data, When vs ASV84"
## [255] "2023-11-15 15:37:45.757654 INFO::Creating boxplot for categorical data, When vs ASV114"
## [256] "2023-11-15 15:37:45.914874 INFO::Creating boxplot for categorical data, When vs ASV130"
## [257] "2023-11-15 15:37:46.030211 INFO::Creating boxplot for categorical data, When vs ASV94"
## [258] "2023-11-15 15:37:46.175879 INFO::Creating boxplot for categorical data, When vs ASV122"
## [259] "2023-11-15 15:37:46.296273 INFO::Creating boxplot for categorical data, When vs ASV138"
## [260] "2023-11-15 15:37:46.430572 INFO::Creating boxplot for categorical data, When vs ASV23"
## [261] "2023-11-15 15:37:46.551883 INFO::Creating boxplot for categorical data, When vs ASV147"
## [262] "2023-11-15 15:37:46.718461 INFO::Creating boxplot for categorical data, When vs ASV66"
## [263] "2023-11-15 15:37:46.861427 INFO::Creating boxplot for categorical data, When vs ASV115"
## [264] "2023-11-15 15:37:46.997245 INFO::Creating boxplot for categorical data, When vs ASV8"
## [265] "2023-11-15 15:37:47.115875 INFO::Creating boxplot for categorical data, When vs ASV164"
## [266] "2023-11-15 15:37:47.28678 INFO::Creating boxplot for categorical data, When vs ASV30"
## [267] "2023-11-15 15:37:47.412871 INFO::Creating boxplot for categorical data, When vs ASV46"
## [268] "2023-11-15 15:37:47.549891 INFO::Creating boxplot for categorical data, When vs ASV126"
## [269] "2023-11-15 15:37:47.695014 INFO::Creating boxplot for categorical data, When vs ASV27"
## [270] "2023-11-15 15:37:47.839722 INFO::Creating boxplot for categorical data, When vs ASV90"
## [271] "2023-11-15 15:37:47.959634 INFO::Creating boxplot for categorical data, When vs ASV137"
## [272] "2023-11-15 15:37:48.095924 INFO::Creating boxplot for categorical data, When vs ASV141"
## [273] "2023-11-15 15:37:48.238626 INFO::Creating boxplot for categorical data, When vs ASV167"
## [274] "2023-11-15 15:37:48.402845 INFO::Creating boxplot for categorical data, When vs ASV111"
## [275] "2023-11-15 15:37:48.522426 INFO::Creating boxplot for categorical data, When vs ASV28"
## [276] "2023-11-15 15:37:48.667893 INFO::Creating boxplot for categorical data, When vs ASV96"
## [277] "2023-11-15 15:37:48.791836 INFO::Creating boxplot for categorical data, When vs ASV85"
## [278] "2023-11-15 15:37:48.943695 INFO::Creating boxplot for categorical data, When vs ASV156"
## [279] "2023-11-15 15:37:49.071795 INFO::Creating boxplot for categorical data, When vs ASV17"
## [280] "2023-11-15 15:37:49.212301 INFO::Creating boxplot for categorical data, When vs ASV116"
## [281] "2023-11-15 15:37:49.347789 INFO::Creating boxplot for categorical data, When vs ASV140"
## [282] "2023-11-15 15:37:49.468278 INFO::Creating boxplot for categorical data, When vs ASV10"
## [283] "2023-11-15 15:37:49.60082 INFO::Creating boxplot for categorical data, When vs ASV105"
## [284] "2023-11-15 15:37:49.720548 INFO::Creating boxplot for categorical data, When vs ASV95"
## [285] "2023-11-15 15:37:50.193691 INFO::Creating boxplot for categorical data, When vs ASV162"
## [286] "2023-11-15 15:37:50.311123 INFO::Creating boxplot for categorical data, When vs ASV9"
## [287] "2023-11-15 15:37:50.449521 INFO::Creating boxplot for categorical data, When vs ASV118"
## [288] "2023-11-15 15:37:50.616939 INFO::Creating boxplot for categorical data, When vs ASV18"
## [289] "2023-11-15 15:37:50.785543 INFO::Creating boxplot for categorical data, When vs ASV121"
## [290] "2023-11-15 15:37:50.937384 INFO::Creating boxplot for categorical data, When vs ASV163"
## [291] "2023-11-15 15:37:51.046009 INFO::Creating boxplot for categorical data, When vs ASV155"
## [292] "2023-11-15 15:37:51.155707 INFO::Creating boxplot for categorical data, When vs ASV148"
## [293] "2023-11-15 15:37:51.313547 INFO::Creating boxplot for categorical data, When vs ASV65"
## [294] "2023-11-15 15:37:51.434255 INFO::Creating boxplot for categorical data, When vs ASV172"
## [295] "2023-11-15 15:37:51.557774 INFO::Creating boxplot for categorical data, When vs ASV128"
## [296] "2023-11-15 15:37:51.671529 INFO::Creating boxplot for categorical data, When vs ASV132"
## [297] "2023-11-15 15:37:51.780464 INFO::Creating boxplot for categorical data, When vs ASV176"
## [298] "2023-11-15 15:37:51.888496 INFO::Creating boxplot for categorical data, When vs ASV53"
## [299] "2023-11-15 15:37:52.030214 INFO::Creating boxplot for categorical data, When vs ASV92"
## [300] "2023-11-15 15:37:52.153231 INFO::Creating boxplot for categorical data, When vs ASV37"
## [301] "2023-11-15 15:37:52.284619 INFO::Creating boxplot for categorical data, When vs ASV146"
## [302] "2023-11-15 15:37:52.409481 INFO::Creating boxplot for categorical data, When vs ASV143"
## [303] "2023-11-15 15:37:52.52513 INFO::Creating boxplot for categorical data, When vs ASV206"
## [304] "2023-11-15 15:37:52.641094 INFO::Creating boxplot for categorical data, When vs ASV33"
## [305] "2023-11-15 15:37:52.778305 INFO::Creating boxplot for categorical data, When vs ASV170"
## [306] "2023-11-15 15:37:52.889769 INFO::Creating boxplot for categorical data, When vs ASV179"
## [307] "2023-11-15 15:37:53.004444 INFO::Creating boxplot for categorical data, When vs ASV199"
## [308] "2023-11-15 15:37:53.117166 INFO::Creating boxplot for categorical data, When vs ASV72"
## [309] "2023-11-15 15:37:53.236631 INFO::Creating boxplot for categorical data, When vs ASV93"
## [310] "2023-11-15 15:37:53.37566 INFO::Creating boxplot for categorical data, When vs ASV193"
## [311] "2023-11-15 15:37:53.503754 INFO::Creating boxplot for categorical data, When vs ASV203"maaslin_all <- read.csv(paste0(getwd(), "/maaslin_output/all_results.tsv"), sep = "\t")
head(maaslin_all)## feature metadata value coef stderr N N.not.0 pval
## 1 ASV41 When Late -2.306703 0.2261749 19 9 1.159387e-08
## 2 ASV36 When Late -3.361515 0.3657296 19 12 5.251227e-08
## 3 ASV50 When Late -3.342318 0.3553378 19 11 3.769650e-08
## 4 ASV61 When Late 2.826934 0.3068431 19 10 5.077002e-08
## 5 ASV44 When Late -3.144538 0.3718139 19 12 1.698523e-07
## 6 ASV25 When Late -3.885594 0.4924141 19 14 4.398006e-07
## qval
## 1 2.052115e-06
## 2 2.323668e-06
## 3 2.323668e-06
## 4 2.323668e-06
## 5 6.012771e-06
## 6 1.112067e-05Volcano plot of MaAslin output
ggplot(maaslin_all, aes(x = coef, y = -log10(qval))) +
geom_point() +
geom_hline(yintercept = -log10(0.05), col = "red") +
annotate("text", x = 0.1, y = 1.5, label = "q-value = 0.1", col = "red")
install.packages("devtools")##
## The downloaded binary packages are in
## /var/folders/qj/y615pkf53ms_hcxg1f9zjkyc0000gp/T//Rtmp1P1C5w/downloaded_packageslibrary(devtools)
devtools::install_github("slowkow/ggrepel")
library(ggrepel)
ggplot(maaslin_all, aes(x = coef, y = -log10(qval))) +
geom_point() +
geom_hline(yintercept = -log10(0.05), col = "red") +
annotate("text", x = 0.1, y = 1.5, label = "q-value = 0.1", col = "red") +
geom_text_repel(aes(label = feature)) +
theme_classic() 
Using SIAMCAT
Due to limited number of samples / metadat in our dataset, we are going to use examples et from SIAMCAT pacakge.
**https://siamcat.embl.de/articles/SIAMCAT_vignette.html*
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SIAMCAT")
#or,
#require("devtools")
#devtools::install_github(repo = 'zellerlab/siamcat')
library(SIAMCAT)Loading example data set
data("feat_crc_zeller", package="SIAMCAT")
data("meta_crc_zeller", package="SIAMCAT")
siamcat_feature <- feat.crc.zeller
siamcat_meta <- meta.crc.zellerPlease note that SIAMCAT is supposed to work with relative abundances. Other types of data (e.g. counts) will also work, but not all functions of the package will result in meaningful outputs.
Checking dimensions of SIAMCAT
dim(siamcat_feature)## [1] 1754 141Generating labels for SIAMCAT analysis
siamcat_label <- create.label(meta=siamcat_meta,
label='Group', case='CRC')
head(siamcat_label)## $label
## CCIS27304052ST-3-0 CCIS15794887ST-4-0 CCIS74726977ST-3-0 CCIS16561622ST-4-0
## -1 -1 -1 -1
## CCIS79210440ST-3-0 CCIS82507866ST-3-0 CCIS98482370ST-3-0 CCIS00281083ST-3-0
## -1 -1 -1 -1
## CCIS32452666ST-4-0 CCIS94417875ST-3-0 CCIS51667829ST-4-0 CCIS27927933ST-4-0
## -1 -1 -1 -1
## CCIS45571137ST-3-0 CCIS81139242ST-4-0 CCIS21126322ST-4-0 CCIS17669415ST-4-0
## -1 -1 -1 1
## CCIS44676181ST-4-0 CCIS29210128ST-4-0 CCIS48725289ST-4-0 CCIS64785924ST-20-0
## -1 -1 1 -1
## CCIS03473770ST-4-0 CCIS10706551ST-3-0 CCIS64773582ST-4-0 CCIS06260551ST-3-0
## -1 -1 -1 1
## CCIS40244499ST-3-0 CCIS85214191ST-3-0 CCIS08668806ST-3-0 CCIS50369211ST-20-0
## 1 1 -1 -1
## CCIS87167916ST-4-0 CCIS94496512ST-4-0 CCIS58234805ST-4-0 CCIS91228662ST-4-0
## 1 -1 1 1
## CCIS71301801ST-4-0 CCIS24898163ST-4-0 CCIS53355328ST-4-0 CCIS41222843ST-4-0
## 1 -1 1 1
## CCIS81887263ST-4-0 CCIS29688262ST-20-0 CCIS87116798ST-4-0 CCIS87605453ST-4-0
## 1 -1 1 1
## CCIS83870198ST-4-0 CCIS61287323ST-4-0 CCIS11362406ST-4-0 CCIS14449628ST-4-0
## 1 1 -1 1
## CCIS21278152ST-4-0 CCIS25399172ST-4-0 CCIS34055159ST-4-0 CCIS76845094ST-20-0
## 1 -1 1 1
## CCIS41548810ST-4-0 CCIS98512455ST-4-0 CCIS83445808ST-4-0 CCIS37250421ST-4-0
## 1 1 -1 -1
## CCIS35092938ST-4-0 CCIS50148151ST-4-0 CCIS62605362ST-3-0 CCIS87252800ST-4-0
## -1 1 1 1
## CCIS22958137ST-20-0 CCIS48507077ST-4-0 CCIS32105356ST-4-0 CCIS81710917ST-20-0
## 1 -1 -1 -1
## CCIS16326685ST-4-0 CCIS53043478ST-4-0 CCIS98832363ST-4-0 CCIS02124300ST-4-0
## -1 1 1 -1
## CCIS44743950ST-4-0 CCIS10793554ST-4-0 CCIS59132091ST-4-0 CCIS96387239ST-4-0
## -1 -1 1 -1
## CCIS82146115ST-4-0 CCIS31434951ST-20-0 CCIS33816588ST-4-0 CCIS12656533ST-4-0
## 1 -1 1 1
## CCIS41806458ST-4-0 CCIS53557295ST-4-0 CCIS47284573ST-4-0 CCIS76563044ST-4-0
## -1 -1 1 -1
## CCIS93040568ST-20-0 CCIS45793747ST-4-0 CCIS95097901ST-4-0 CCIS80834637ST-4-0
## -1 -1 -1 -1
## CCIS63910149ST-4-0 CCIS84543192ST-4-0 CCIS38765456ST-20-0 CCIS82944710ST-20-0
## -1 1 1 -1
## CCIS55230578ST-4-0 CCIS24254057ST-4-0 CCIS78318719ST-4-0 CCIS44757994ST-4-0
## 1 1 1 1
## CCIS65479369ST-4-0 CCIS62794166ST-4-0 CCIS77252613ST-4-0 CCIS11015875ST-4-0
## 1 1 1 1
## CCIS55531770ST-4-0 CCIS70398272ST-4-0 CCIS35100175ST-4-0 CCIS22416007ST-4-0
## -1 -1 1 1
## CCIS41288781ST-4-0 CCIS59903910ST-4-0 CCIS12370844ST-4-0 CCIS72607085ST-4-0
## 1 -1 1 1
## CCIS46047672ST-4-0 CCIS22906510ST-20-0 CCIS52370277ST-4-0 CCIS02379307ST-4-0
## 1 -1 1 1
## CCIS48579360ST-4-0 CCIS46467422ST-4-0 CCIS54027808ST-4-0 CCIS15704761ST-4-0
## -1 1 -1 1
## CCIS11558985ST-4-0 CCIS13047523ST-4-0 CCIS94603952ST-4-0 CCIS02856720ST-4-0
## -1 -1 -1 -1
## CCIS07648107ST-4-0 CCIS34604008ST-4-0 CCIS63468405ST-4-0 CCIS16383318ST-4-0
## -1 -1 -1 -1
## CCIS36699628ST-4-0 CCIS95409808ST-4-0 CCIS74239020ST-4-0 CCIS05314658ST-4-0
## -1 -1 -1 -1
## CCIS71578391ST-4-0 CCIS50003399ST-4-0 CCIS36797902ST-4-0 CCIS11354283ST-4-0
## -1 -1 -1 -1
## CCIS50471204ST-4-0 CCIS20795251ST-4-0 CCIS07277498ST-4-0 CCIS90164298ST-4-0
## -1 -1 -1 -1
## CCIS23164343ST-4-0 CCIS03857607ST-4-0 CCIS63448910ST-4-0 CCIS88007743ST-4-0
## -1 -1 -1 -1
## CCIS51595129ST-4-0 CCIS09568613ST-4-0 CCIS07539127ST-4-0 CCIS41692898ST-4-0
## -1 -1 -1 -1
## CCIS44093303ST-4-0 CCIS00146684ST-4-0 CCIS78100604ST-4-0 CCIS48174381ST-4-0
## -1 -1 -1 -1
## CCIS88317640ST-4-0
## -1
##
## $info
## CTR CRC
## -1 1
##
## $type
## [1] "BINARY"Creating a SIAMCAT object
sc.obj <- siamcat(feat=siamcat_feature,
label=siamcat_label,
meta=siamcat_meta)
sc.obj## siamcat-class object
## label() Label object: 88 CTR and 53 CRC samples
##
## contains phyloseq-class experiment-level object @phyloseq:
## phyloseq@otu_table() OTU Table: [ 1754 taxa and 141 samples ]
## phyloseq@sam_data() Sample Data: [ 141 samples by 6 sample variables ]Filtering low abundant taxta
sc.obj <- filter.features(sc.obj,
filter.method = 'abundance',
cutoff = 0.001)Testing associations
Associations between microbial species and the label can be tested with the check.associations function. The function computes for each species the significance using a non-parametric Wilcoxon test and different effect sizes for the association (e.g. AUC or fold change).
The function produces a pdf file as output, since the plot is optimized for a landscape DIN-A4 layout, but can also used to plot on an active graphic device, e.g. in RStudio. The resulting plot then looks like that:
sc.obj <- check.associations(sc.obj, log.n0 = 1e-06, alpha = 0.05)
association.plot(sc.obj, sort.by = 'fc',
panels = c('fc', 'prevalence', 'auroc'))## Are you sure that you want to continue? (Yes/no/cancel)
Confounder testing
check.confounders(sc.obj, fn.plot = 'confounder_plots.pdf',
meta.in = NULL, feature.type = 'filtered')Data normalization
sc.obj <- normalize.features(sc.obj, norm.method = "log.unit",
norm.param = list(log.n0 = 1e-06, n.p = 2,norm.margin = 1))Cross-validation preparing
sc.obj <- create.data.split(sc.obj, num.folds = 5, num.resample = 2)Training model using lasso (least absolute shrinkage and selection operator; also Lasso or LASSO), a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the resulting statistical model.
Lasso is a supervised regularization method used in machine learning.
method = specifies the type of model to be trained, may
be one of these:
c('lasso', 'enet', 'ridge', 'lasso_ll', 'ridge_ll', 'randomForest')
#Model training
sc.obj.lasso <- train.model(sc.obj, method = "lasso")Model evaluation
This function calculates the Area Under the Receiver Operating Characteristic (ROC) Curve (AU-ROC) )and the Precision Recall (PR) Curve for each resampled cross-validation run.
#Making predictions
sc.obj.lasso <- make.predictions(sc.obj.lasso)
sc.obj.lasso <- evaluate.predictions(sc.obj.lasso)
model.evaluation.plot(sc.obj.lasso)
Making interpretation plot
model.interpretation.plot(sc.obj.lasso, fn.plot = 'interpretation.pdf',
consens.thres = 0.5, limits = c(-3, 3), heatmap.type = 'zscore')Bibliography
## Computing. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>. We have invested a lot of time and effort in creating R, please cite it when using it for data analysis. See also 'citation("pkgname")' for citing R packages.
## Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
## J, reikoch, Beasley W, O'Connor B, Warnes GR, Quinn M, Kamvar ZN (2023). yaml: Methods to Convert R Data to YAML and Back_. R package version 2.3.7, <https://CRAN.R-project.org/package=yaml>. ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see 'help("citation")'.
## Generation in R. R package version 1.45, <https://yihui.org/knitr/>. Yihui Xie (2015) Dynamic Documents with R and knitr. 2nd edition. Chapman and Hall/CRC. ISBN 978-1498716963 Yihui Xie (2014) knitr: A Comprehensive Tool for Reproducible Research in R. In Victoria Stodden, Friedrich Leisch and Roger D. Peng, editors, Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595
## Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2023). rmarkdown: Dynamic Documents for R. R package version 2.25, <https://github.com/rstudio/rmarkdown>. Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, <https://bookdown.org/yihui/rmarkdown>. Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook_. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, <https://bookdown.org/yihui/rmarkdown-cookbook>.
## 'rmarkdown' Documents. R package version 1.0.4, <https://CRAN.R-project.org/package=rmdformats>.