.Rhistory

tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
tibble::add_column(., timepoint = stringr::str_sub(.$SampleID, 12, 13)) %>%
dplyr::group_by(timepoint, OTU) %>%
dplyr::select(-SampleID)%>%
dplyr::summarise(Abundance=mean(Abundance)) %>%
tidyr::spread(timepoint, Abundance) %>%
tibble::column_to_rownames(var = 'OTU') %>%
data.frame()
gg_topi_meta_samples <- phyloseq(otu_table(meta_samples, taxa_are_rows = TRUE), tax_table(gg_topi), phy_tree(gg_topi))
otus <- data.frame(otu_table(gg_topi.genus.relab), check.names = FALSE) %>%
tibble::rownames_to_column('OTU') %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
left_join(data.frame(tax_table(gg_topi.genus.relab)[,-7], check.names = FALSE) %>% tibble::rownames_to_column('OTU')) %>%
tibble::add_column(., timepoint = stringr::str_sub(.$SampleID, 12, 13)) %>%
tibble::add_column(., mouseID = stringr::str_sub(.$SampleID, 8, 10))
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU ~ timepoint, value.var = 'Abundance')
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU ~ timepoint + mouseID, value.var = 'Abundance')
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU + ...~ timepoint + mouseID, value.var = 'Abundance')
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU ~ timepoint + mouseID + ... , value.var = 'Abundance')
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU ~ timepoint + mouseID , value.var = 'Abundance')
otus <- data.frame(otu_table(gg_topi.genus.relab), check.names = FALSE) %>%
tibble::rownames_to_column('OTU') %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
tibble::add_column(., timepoint = stringr::str_sub(.$SampleID, 12, 13)) %>%
tibble::add_column(., mouseID = stringr::str_sub(.$SampleID, 8, 10))
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU ~ timepoint + mouseID , value.var = 'Abundance')
otus
otus %>%
# filter(mouseID==865) %>%
# head %>%
reshape2::dcast(OTU ~ timepoint + mouseID , value.var = 'Abundance') %>%
left_join(data.frame(tax_table(gg_topi.genus.relab)[,-7], check.names = FALSE) %>% tibble::rownames_to_column('OTU'))
otus %>%
filter(mouseID==865) %>%
reshape2::dcast(OTU ~ timepoint + mouseID , value.var = 'Abundance') %>%
left_join(data.frame(tax_table(gg_topi.genus.relab)[,-7], check.names = FALSE) %>% tibble::rownames_to_column('OTU'))
otus %>%
filter(mouseID==865) %>%
reshape2::dcast(OTU ~ timepoint, value.var = 'Abundance') %>%
left_join(data.frame(tax_table(gg_topi.genus.relab)[,-7], check.names = FALSE) %>% tibble::rownames_to_column('OTU'))
libray
library(xlsx)
install.packages('xlsx')
library(xlsx)
library(xlsx)
library(xlsx)
library(xlsx)
library(xlsx)
Sys.getenv(c("JAVA_HOME", "LD_LIBRARY_PATH"))
install.packages('rJava')
library(xlsx)
remove.packages(c('rJava', 'xlsx'))
remove.packages(c('rJava', 'xlsx'))
library(xlsx)
gg_topi <- import_biom(BIOMfilename = '/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/otus/otu_table_mc10_w_tax_rdp.biom',
treefilename = '/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/otus/rep_set.tre',
refseqfilename = '/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/otus/rep_set.fna',
refseqFunction = parse_taxonomy_default
)
knitr::opts_chunk$set(message = FALSE, warning = FALSE, paged.print=FALSE, echo=FALSE)
knitr::opts_chunk$set(fig.width=12, fig.height=8)
library(phyloseq)
library(magrittr)
library(tidyverse)
library(reshape2)
library(ggsignif)
library(DESeq2)
gg_topi <- import_biom(BIOMfilename = '/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/otus/otu_table_mc10_w_tax_rdp.biom',
treefilename = '/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/otus/rep_set.tre',
refseqfilename = '/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/otus/rep_set.fna',
refseqFunction = parse_taxonomy_default
)
mappingfile <- read.delim('/shares/CIBIO-Storage/CM/mir/projects/fbeghini_16s/20180731_run186ms_segata/grandi/mappingfile', row.names = 1, check.names = FALSE)
mappingfile$Description %<>% str_replace_all(.,'_','-')
mappingfile$condizione <- factor(mappingfile$condizione, labels = c('Pentatome + Bifido', 'Pentatome + PBS', 'Empty + Bifido', 'Empty + PBS'))
sample_names(gg_topi) <- mappingfile[sample_names(gg_topi), 'Description']
rownames(mappingfile) <- mappingfile$Description
sample_data(gg_topi) <- mappingfile
sample_data(gg_topi)$id.topo <- factor(sample_data(gg_topi)$id.topo)
sample_data(gg_topi)$pentatome <- factor(stringr::str_split(sample_data(gg_topi)$condizione, ' \\+ ', simplify = TRUE)[,1])
sample_data(gg_topi)$bifido <- factor(stringr::str_split(sample_data(gg_topi)$condizione, ' \\+ ', simplify = TRUE)[,2])
sample_data(gg_topi)$cond_all <- factor(with(sample_data(gg_topi), paste(timepoint,real_vaccine, real_bifido, sep = '_')))
sample_data(gg_topi)$cond_vac <- factor(with(sample_data(gg_topi), paste(merge_timepoint,real_vaccine, sep = '_')))
sample_data(gg_topi)$cond_bif <- factor(with(sample_data(gg_topi), paste(merge_timepoint,real_bifido, sep = '_')))
sample_data(gg_topi)$vacc_bif <- strtoi(paste(ifelse(sample_data(gg_topi)$real_vaccine=='YES','1','0'),ifelse(sample_data(gg_topi)$real_bifido=='YES','1','0'), sep = ''), 2)
colnames(tax_table(gg_topi)) <- c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species")
tax_table_filtered <- tax_table(gg_topi) %>%
as.data.frame %>%
mutate_all( .funs = function(x) stringr:::str_replace_all(x,"D_.__",""))
rownames(tax_table_filtered) <- rownames(tax_table(gg_topi))
tax_table(gg_topi) <-  as.matrix(tax_table_filtered)
unass_otus <- tax_table(gg_topi) %>%
data.frame %>%
tibble::rownames_to_column('OTU') %>%
filter(Domain == 'Unassigned') %>%
select(OTU)
unknown_otus <- tax_table(gg_topi) %>%
data.frame %>%
tibble::rownames_to_column('OTU') %>%
filter(is.na(Genus)) %>%
select(OTU)
unknown_family <- tax_table(gg_topi) %>%
data.frame %>%
tibble::rownames_to_column('OTU') %>%
filter(is.na(Family)) %>%
select(OTU)
tax_table(gg_topi)[unass_otus$OTU, ] <- rep.int('Unassigned', 7)
tax_table(gg_topi)[unknown_otus$OTU, 'Genus'] <- 'Unknown'
tax_table(gg_topi)[unknown_family$OTU, 'Family'] <- 'Unknown'
gg_topi.relab <- transform_sample_counts(gg_topi, function(x) x/sum(x))
gg_topi.genus <- tax_glom(gg_topi,taxrank = "Genus")
gg_topi.genus.relab <- transform_sample_counts(gg_topi.genus, function(x) x/sum(x))
gg_topi.family <- tax_glom(gg_topi,taxrank = "Genus")
gg_topi.family.relab <- transform_sample_counts(gg_topi.family, function(x) x/sum(x))
ggtopi_subs <-subset_samples(gg_topi, timepoint %in% c('T3','T7'))
ggtopi_subs.genus <- tax_glom(ggtopi_subs,taxrank = "Genus")
ggtopi_subs.genus.relab <- transform_sample_counts(ggtopi_subs.genus, function(x) x/sum(x))
mmtx <- model.matrix(object = ~timepoint + timepoint:real_vaccine + timepoint:real_bifido + timepoint:real_vaccine:real_bifido, data = data.frame(ggtopi_subs@sam_data))
idx <- which(apply(mmtx, 2, function(x) all(x==0)))
mmtx <- mmtx[,-idx]
pds_37 <- phyloseq_to_deseq2(ggtopi_subs, design = ~ timepoint + timepoint:real_vaccine + timepoint:real_bifido + timepoint:real_vaccine:real_bifido, ignoreRank = TRUE)
dds_37 <- DESeq2::DESeq(pds_37, parallel = TRUE, quiet = TRUE, full = mmtx)
pds_genus <- phyloseq_to_deseq2(ggtopi_subs.genus, design =  ~ timepoint + timepoint:real_vaccine + timepoint:real_bifido + timepoint:real_vaccine:real_bifido, ignoreRank = TRUE)
dds_genus <- DESeq2::DESeq(pds_genus, parallel = TRUE, quiet = TRUE, full = mmtx)
res37 <- DESeq2::results(dds_37, name = "timepointT7", alpha = 0.05)
res37.filt <- res37[!is.na(res37$padj),]
res37.filt <- res37.filt[res37.filt$padj<0.05,]
res37.filt
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs)[rownames(res37.filt),])) %>%
ggplot(aes(log2FoldChange, interaction(Family, Genus))) +
geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
geom_point(size=3, aes(color = baseMean)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y",switch = "y",as.table = T)+
theme_bw() +
scale_y_discrete(position = "right") +
theme(axis.text = element_text(hjust = 0, vjust=0.5, size = 11),
strip.text = element_text(size = 9),
axis.title = element_text(size = 9),
axis.title.y = element_blank(),
strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 180)
)
t37_samples_idx <- grep('T7|T3', sample_names(ggtopi_subs), perl = TRUE)
otu_ab_t37 <- otu_table(ggtopi_subs)[rownames(res37.filt), t37_samples_idx] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13))
otu_ab_t37
otu_ab_t37 %>% bind_cols(data.frame(tax_table(ggtopi_filt)[otu_ab_t37$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
otu_ab_t37 %>% bind_cols(data.frame(tax_table(ggtopi_subs)[otu_ab_t37$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(gg_topi.relab)[rownames(res37.filt),])) %>%
ggplot(aes(log2FoldChange, interaction(Family, Genus))) +
geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
geom_point(size=3, aes(color = baseMean)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y",switch = "y",as.table = T)+
theme_bw() +
scale_y_discrete(position = "right") +
theme(axis.text = element_text(hjust = 0, vjust=0.5, size = 11),
strip.text = element_text(size = 9),
axis.title = element_text(size = 9),
axis.title.y = element_blank(),
strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 180)
)
otu_ab_t37 <- otu_table(gg_topi.relab)[rownames(res37.filt), t37_samples_idx] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13))
otu_ab_t37 %>% bind_cols(data.frame(tax_table(gg_topi.relab)[otu_ab_t37$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
t37_samples_idx <- grep('T7|T3', sample_names(ggtopi_subs), perl = TRUE)
otu_ab_t37 <- otu_table(gg_topi.relab)[rownames(res37.filt), t37_samples_idx] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13))
otu_ab_t37
otu_ab_t37 %>% bind_cols(data.frame(tax_table(gg_topi.relab)[otu_ab_t37$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
otu_ab_t37
otu_ab_t37 %>% bind_cols(data.frame(tax_table(ggtopi_subs)[otu_ab_t37$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
t37_samples_idx <- grep('T7|T3', sample_names(ggtopi_subs), perl = TRUE)
ggtopi_subs.relab <- transform_sample_counts(ggtopi_subs, function(x) x/sum(x))
otu_ab_t37 <- otu_table(ggtopi_subs.relab)[rownames(res37.filt), t37_samples_idx] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13))
otu_ab_t37 %>% bind_cols(data.frame(tax_table(ggtopi_subs.relab)[otu_ab_t37$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs)[rownames(res37.filt),]))
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs)[rownames(res37.filt),])) %>% filter(Genus %in% c('Escherichia'))
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs)[rownames(res37.filt),])) %>% filter(Genus %in% c('Escherichia-Shigella'))
res37 <- DESeq2::results(dds_genus, name = "timepointT7", alpha = 0.05)
res37.filt <- res37[!is.na(res37$padj),]
res37.filt <- res37.filt[res37.filt$padj<0.05,]
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs.genus)[rownames(res37.filt),])) %>%
ggplot(aes(log2FoldChange, interaction(Family, Genus))) +
geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
geom_point(size=3, aes(color = baseMean)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y",switch = "y",as.table = T)+
theme_bw() +
scale_y_discrete(position = "right") +
theme(axis.text = element_text(hjust = 0, vjust=0.5, size = 11),
strip.text = element_text(size = 9),
axis.title = element_text(size = 9),
axis.title.y = element_blank(),
strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 180)
)
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs.genus)[rownames(res37.filt),-7])) %>%
ggplot(aes(log2FoldChange, interaction(Family, Genus))) +
geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
geom_point(size=3, aes(color = baseMean)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y",switch = "y",as.table = T)+
theme_bw() +
scale_y_discrete(position = "right") +
theme(axis.text = element_text(hjust = 0, vjust=0.5, size = 11),
strip.text = element_text(size = 9),
axis.title = element_text(size = 9),
axis.title.y = element_blank(),
strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 180)
)
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13))
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
gg_topi_no247 <- subset_taxa(gg_topi, !Family %in% c("Bacteroidales S24-7 group"))
ggtopi_subs <-subset_samples(gg_topi_no247, timepoint %in% c('T3','T7'))
ggtopi_subs.genus <- tax_glom(ggtopi_subs,taxrank = "Genus")
ggtopi_subs.genus.relab <- transform_sample_counts(ggtopi_subs.genus, function(x) x/sum(x))
mmtx <- model.matrix(object = ~timepoint + timepoint:real_vaccine + timepoint:real_bifido + timepoint:real_vaccine:real_bifido, data = data.frame(ggtopi_subs@sam_data))
idx <- which(apply(mmtx, 2, function(x) all(x==0)))
mmtx <- mmtx[,-idx]
pds_37 <- phyloseq_to_deseq2(ggtopi_subs, design = ~ timepoint + timepoint:real_vaccine + timepoint:real_bifido + timepoint:real_vaccine:real_bifido, ignoreRank = TRUE)
dds_37 <- DESeq2::DESeq(pds_37, parallel = TRUE, quiet = TRUE, full = mmtx)
pds_genus <- phyloseq_to_deseq2(ggtopi_subs.genus, design =  ~ timepoint + timepoint:real_vaccine + timepoint:real_bifido + timepoint:real_vaccine:real_bifido, ignoreRank = TRUE)
dds_genus <- DESeq2::DESeq(pds_genus, parallel = TRUE, quiet = TRUE, full = mmtx)
res37 <- DESeq2::results(dds_genus, name = "timepointT7", alpha = 0.05)
res37.filt <- res37[!is.na(res37$padj),]
res37.filt <- res37.filt[res37.filt$padj<0.05,]
res37.filt
res37.filt %>% as.data.frame %>%
cbind(as.matrix(tax_table(ggtopi_subs.genus)[rownames(res37.filt),-7])) %>%
ggplot(aes(log2FoldChange, interaction(Family, Genus))) +
geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
geom_point(size=3, aes(color = baseMean)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y",switch = "y",as.table = T)+
theme_bw() +
scale_y_discrete(position = "right") +
theme(axis.text = element_text(hjust = 0, vjust=0.5, size = 11),
strip.text = element_text(size = 9),
axis.title = element_text(size = 9),
axis.title.y = element_blank(),
strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 180)
)
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13))
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
genus_ab_t
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')]))
order(genus_ab_t[['Genus']], genus_ab_t$logFC, decreasing=TRUE)
order(genus_ab_t[[Genus]], genus_ab_t$logFC, decreasing=TRUE)
order(genus_ab_t[[Genus]], genus_ab_t$Abundance, decreasing=TRUE)
order(genus_ab_t[['Genus']], genus_ab_t$Abundance, decreasing=TRUE)
order(genus_ab_t$Abundance, decreasing=TRUE)
order(genus_ab_t[['Genus']], genus_ab_t$Abundance, decreasing=TRUE)
genus_ab_t
genus_ab_t %<>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')]))
order(genus_ab_t[['Genus']], genus_ab_t$Abundance, decreasing=TRUE)
ord <- order(genus_ab_t[['Genus']], genus_ab_t$Abundance, decreasing=TRUE)
genus_ab_t$Genus <- factor(genus_ab_t$Genus, levels = unique(genus_ab_t$Genus[ord]))
genus_ab_t$Genus <- factor(genus_ab_t$Genus, levels = unique(genus_ab_t$Genus[ord]))
genus_ab_t$Genus
ggplot(genus_ab_t) +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0)) +
theme_bw()
ggplot(genus_ab_t) +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme_bw() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
genus_ab_t
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=timepoint)) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme_bw() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
genus_ab_t
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')]))
stringr::str_sub(.$SampleID,12, 13)
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 7))
genus_ab_t
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6))
genus_ab_t
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=interaction(timepoint, condition)) +
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')])) %>%
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=interaction(timepoint, condition))) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme_bw() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
?mutate_if
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6)) %>%
mutate_if(timepoint == 'T3', condition = '')
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6)) %>%
mutate_if('timepoint' == 'T3', condition = '')
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6)) %>%
mutate_if(.$timepoint == 'T3', condition = '')
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6)) %>%
mutate_if(.$timepoint == 'T3', .$condition = '')
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6)) %>%
mutate_if(equals(timepoint, 'T3'), condition = '')
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition=stringr::str_sub(.$SampleID,4, 6)) %>%
mutate_if(equals(.$timepoint, 'T3'), condition = '')
genus_ab_t <- otu_table(ggtopi_subs.genus.relab)[rownames(res37.filt), ] %>%
data.frame %>%
tibble::rownames_to_column("OTU") %>%
tidyr::gather_("SampleID", "Abundance", setdiff(colnames(.), "OTU")) %>%
mutate(timepoint=stringr::str_sub(.$SampleID,12, 13)) %>%
mutate(condition= if_else(timepoint=='T3','',stringr::str_sub(.$SampleID,4, 6)))
genus_ab_t
genus_ab_t %>% bind_cols(data.frame(tax_table(ggtopi_subs.genus.relab)[genus_ab_t$OTU,c('Phylum','Family','Genus')])) %>%
ggplot() +
geom_boxplot(aes(interaction(Family, Genus), Abundance, fill=interaction(timepoint, condition))) +
facet_grid(Phylum ~ ., scales = "free", space = "free_y", as.table= T)+
coord_flip() +
theme_bw() +
theme(strip.background = element_rect(colour = "transparent", fill = "transparent"),
strip.text.y = element_text(angle = 0))
ggplot(x) +
geom_boxplot(aes(timepoint, value)) +
facet_grid(variable~condizione,scales = 'free_y')
knitr::opts_chunk$set(message = FALSE, warning = FALSE, paged.print=FALSE, echo=FALSE)
knitr::opts_chunk$set(fig.width=12, fig.height=8)
library(phyloseq)
library(magrittr)
library(tidyverse)
library(reshape2)
library(ggsignif)
library(DESeq2)
ggplot(x) +
geom_boxplot(aes(timepoint, value)) +
facet_grid(variable~condizione,scales = 'free_y')
x <- melt(alpha_div, id.vars = c('Description','condizione','timepoint'), measure.vars =  c('Observed', 'se.chao1', 'Simpson', 'Shannon'))
alpha_div <- estimate_richness(gg_topi)
alpha_div %<>% bind_cols(data.frame(sample_data(gg_topi)))
x <- melt(alpha_div, id.vars = c('Description','condizione','timepoint'), measure.vars =  c('Observed', 'se.chao1', 'Simpson', 'Shannon'))
x$value <- as.numeric(x$value)
ggplot(x) +
geom_boxplot(aes(timepoint, value)) +
facet_grid(variable~condizione,scales = 'free_y')
filt_alpha <- alpha_div %>%
filter(timepoint %in% c('T2','T3')) %>%
select('Description','condizione','timepoint','Observed', 'se.chao1', 'Simpson', 'Shannon') %>%
mutate(t=as.numeric(stringr::str_sub(timepoint,2,2)))
ttest_res <- data.frame()
for(measure in c('Observed', 'se.chao1', 'Simpson', 'Shannon')){
ttest_res %<>% bind_rows(data.frame(Measure=measure, pvalue=with(filt_alpha, t.test(as.formula(paste(measure, 't', sep = '~'))))$p.value))
}
filt_alpha %>%
melt(id.vars = c('Description','condizione','timepoint'), measure.vars =  c('Observed', 'se.chao1', 'Simpson', 'Shannon')) %>%
ggplot() +
geom_boxplot(aes(timepoint, value)) +
facet_grid(variable~.,scales = 'free_y')
knitr::kable(ttest_res)
source('R/loadQiimeData.R')
# Chunk 1: setup
knitr::opts_chunk$set(message = FALSE, warning = FALSE, paged.print=FALSE, echo=FALSE)
knitr::opts_chunk$set(fig.width=12, fig.height=8)
library(phyloseq)
library(magrittr)
library(tidyverse)
library(reshape2)
library(ggsignif)
library(DESeq2)
source('R/loadQiimeData.R')
gg_topi <- loadQiimeData()
gg_topi
devtools::create('ggtopi')