Nothing
core.OTU <- function(data, meta, meta.factor="", percent=1){
if ( class(data) != "list" ) {
stop("please provide otu tables as list. See ?RAM.input.formatting for details")
}
.valid.data(data=data, is.OTU=TRUE)
list.core <-list()
labels <- names(data)
for ( i in 1:length(data) ) {
label <- names(data)[i]
otu <- data[[i]]
if ( is.null(otu) ) { break }
# valid otu & meta
valid.OTU(otu)
#.valid.meta(otu1=otu, meta=meta)
samp <- colnames(otu)[1:(ncol(otu)-1)]
meta <- meta[match(samp, rownames(meta)), ]
if ( !identical(rownames(meta), samp) ) {
stop("Error: metadata and otu table do not have same samples")
}
# check meta category
.valid.factor(meta, meta.factor)
fac.levels <- levels(factor(meta[[meta.factor]]))
otu <- LCA.OTU(otu, strip.format=FALSE, drop=TRUE)
otu$LCA <- paste(otu$LCA, rownames(otu), sep="_")
###check for zero entries?
otu.fac.list=list()
otu.fac.id = list()
otu.fac.rank = list()
for (j in 1:length(fac.levels)) {
# each fac levels
meta.fac <- meta[grep(fac.levels[j], meta[[meta.factor]]), ]
# presence/absence
otu.pa <- decostand(otu[, -ncol(otu)], "pa")
# select only otus present in more than percent of samples
# in each category
col <- which(colnames(otu.pa) %in% rownames(meta.fac))
if ( length(col) == 1 ) {
vec <- otu.pa[, col]
sel <- which(vec == 1 )
} else {
sel <- which(rowSums(otu.pa[, col, drop=FALSE]) >= round(percent*nrow(meta.fac)))
}
otu.fac <- otu[, c(col, ncol(otu))]
if ( length(sel) != 0) {
otu.fac.ids <- as.character(rownames(otu.fac)[sel])
otu.fac.id[[j]] <- as.character(rownames(otu.fac)[sel])
otu.fac.ranks <- unique(as.character(factor(otu.fac[["LCA"]][sel])))
otu.fac.rank[[j]] <- unique(as.character(factor(otu.fac[["LCA"]][sel])))
otu.fac.ids.num <- length(otu.fac.ids)
} else {
otu.fac.ids <- ""
otu.fac.id[[j]] <- ""
otu.fac.ranks <- ""
otu.fac.rank[[j]] <- ""
otu.fac.ids.num <- 0
}
# summarize
description <- paste(otu.fac.ids.num, "_otus_found_in_",
percent*100, "%_",
fac.levels[j], "_samples;", sep="")
otu.fac.list[[j]] <- list(summary=description, otuID=otu.fac.ids,
LCA=otu.fac.ranks)
names(otu.fac.list)[[j]] <- fac.levels[j]
}
#return(otu.fac.list)
# find otus across all categories
otu.fac.list.num <- length(otu.fac.list)
otu.fac.id[[otu.fac.list.num+1]] <- as.character(factor(
Reduce(intersect, otu.fac.id)))
otu.fac.rank[[otu.fac.list.num+1]] <- as.character(factor(
Reduce(intersect, otu.fac.rank)))
#return(list(otu.fac.id, otu.fac.list.num))
# otu ids in all samples
# percent of sequences in core otus across all categories
if ( length(otu.fac.id[[otu.fac.list.num+1]]) == 0 ||
unique(unique(otu.fac.id[[otu.fac.list.num+1]]) == "") ) {
otu.fac.id.all.num <- 0
otu.fac.rank.all.num <- 0
otu.fac.perc <- 0
} else {
otu.fac.id.all.num <- length(unique(otu.fac.id[[otu.fac.list.num+1]]))
otu.fac.rank.all.num <- length(unique(otu.fac.rank[[otu.fac.list.num+1]]))
otu.sel <- otu[rownames(otu) %in% otu.fac.id[[otu.fac.list.num+1]], ]
otu.fac.perc <- 100*sum(otu.sel[, -ncol(otu.sel)])/sum(otu[, -dim(otu)[2]])
}
description <- paste(otu.fac.id.all.num, "_otus_present_in_core_OTU_of_each_",
meta.factor, sep="")
otu.fac.list[[otu.fac.list.num+1]] <- list(summary=description,
otuID=otu.fac.id[[otu.fac.list.num+1]],
LCA=otu.fac.rank[[otu.fac.list.num+1]])
names(otu.fac.list)[[otu.fac.list.num+1]] <- meta.factor
list.core[[label]] <- otu.fac.list
}
return(list.core)
}
core.OTU.rank<-function(data, rank="g", drop.unclassified=TRUE,
meta, meta.factor="", percent=1){
if ( class(data) != "list" ) {
stop("please provide otu tables as list. See ?RAM.input.formatting for details")
}
.valid.data(data=data, is.OTU=TRUE)
list.core <-list()
labels <- names(data)
for ( i in 1:length(data) ) {
label <- names(data)[i]
otu <- data[[i]]
if ( is.null(otu) ) { break }
# valid otu & meta
valid.OTU(otu)
#.valid.meta(otu1=otu, meta=meta)
samp <- colnames(otu)[1:(ncol(otu)-1)]
meta <- meta[match(samp, rownames(meta)), ]
if ( !identical(rownames(meta), samp) ) {
stop("Error: metadata and otu table do not have same samples")
}
# check meta category
.valid.factor(meta, meta.factor)
fac.levels <- levels(factor(meta[[meta.factor]]))
# get rank name & pattern
rank.name <- .get.rank.name(rank, plural=TRUE)
rank_name <- .get.rank.name(rank)
rank_pat <- .get.rank.pat(rank)
# tax.split
otu.all <- tax.split(otu1=otu, rank=rank)
otu <- otu.all[otu.all[[rank_name]] !="", ]
rm <- paste(.blacklist(), collapse="|")
if ( drop.unclassified ) {
rm.row <- which(grepl(rm, otu[[rank_name]], ignore.case=TRUE))
if ( length(rm.row) == 0 ) {
rm.row <- NULL
} else {
rm.row <- rm.row
}
} else {
rm.row <- NULL
}
if ( !is.null(rm.row) ) {
otu <- otu[-rm.row, ]
} else {
otu <- otu
}
# presence/absence
otu.pa <- decostand(otu[, -ncol(otu)], "pa")
###check for zero entries?
otu.fac.list=list()
otu.fac.id = list()
otu.fac.rank = list()
for (j in 1:length(fac.levels)) {
# each fac levels
meta.fac <- meta[grep(fac.levels[j], meta[[meta.factor]]), ]
# select only otus present in more than percent of samples
# in each category
col <- which(colnames(otu.pa) %in% rownames(meta.fac))
if ( length(col) == 1 ) {
vec <- otu.pa[, col]
sel <- which(vec == 1 )
} else {
sel <- which(rowSums(otu.pa[, col]) >= round(percent*nrow(meta.fac)))
}
otu.fac <- otu[, c(col, ncol(otu))]
otu.all.fac <- otu.all[, c(col, ncol(otu.all))]
if ( length(sel) != 0) {
otu.fac.ids <- as.character(rownames(otu.fac)[sel])
otu.fac.id[[j]] <- as.character(rownames(otu.fac)[sel])
otu.fac.ranks <- unique(as.character(factor(otu.fac[[rank_name]][sel])))
otu.fac.rank[[j]] <- unique(as.character(factor(otu.fac[[rank_name]][sel])))
otu.fac.ids.num <- length(otu.fac.ids)
otu.fac.ranks.num <- length(otu.fac.ranks)
otu.fac.ids.perc <- 100*sum(otu.fac[c(otu.fac.ids),
-ncol(otu.fac)])/sum(otu.all.fac[,
-ncol(otu.all.fac)])
} else {
otu.fac.ids <- ""
otu.fac.id[[j]] <- ""
otu.fac.ranks <- ""
otu.fac.rank[[j]] <- ""
otu.fac.ids.num <- 0
otu.fac.ranks.num <- 0
otu.fac.ids.perc <- 0
}
# summarize
description <- paste(otu.fac.ids.num, "_otus_found_in_", percent*100, "%_",
fac.levels[j], "_samples; ","assigned_to_", otu.fac.ranks.num, "_",
rank.name, "; ", otu.fac.ids.perc, "%_of_total_sequences_in_",
fac.levels[j], sep="")
otu.fac.list[[j]]<- list(summary=description, otuID=otu.fac.ids, taxa=otu.fac.ranks)
names(otu.fac.list)[[j]] <- fac.levels[j]
}
# find otus across all categories
otu.fac.list.num <- length(otu.fac.list)
otu.fac.id[[otu.fac.list.num+1]] <- as.character(factor(Reduce(intersect, otu.fac.id)))
otu.fac.rank[[otu.fac.list.num+1]] <- as.character(factor(Reduce(intersect, otu.fac.rank)))
# otu ids in all samples
# percent of sequences in core otus across all categories
if ( length(otu.fac.id[[otu.fac.list.num+1]]) == 0 ||
unique(unique(otu.fac.id[[otu.fac.list.num+1]]) == "") ) {
otu.fac.id.all.num <- 0
otu.fac.rank.all.num <- 0
otu.fac.perc <- 0
} else {
otu.fac.id.all.num <- length(unique(otu.fac.id[[otu.fac.list.num+1]]))
otu.fac.rank.all.num <- length(unique(otu.fac.rank[[otu.fac.list.num+1]]))
otu.sel <- otu[rownames(otu) %in% otu.fac.id[[otu.fac.list.num+1]], ]
otu.fac.perc <- 100*sum(otu.sel[, -ncol(otu.sel)])/sum(otu.all[, -dim(otu.all)[2]])
}
description <- paste(otu.fac.id.all.num, "_otus_present_in_core_OTU_of_each_",
meta.factor, "; ", "assigned_to_", otu.fac.rank.all.num, "_",
rank.name, "; ", otu.fac.perc, "%_of_total_sequences",sep="")
otu.fac.list[[otu.fac.list.num+1]] <- list(summary=description,
otuID=otu.fac.id[[otu.fac.list.num+1]],
taxa=otu.fac.rank[[otu.fac.list.num+1]])
names(otu.fac.list)[[otu.fac.list.num+1]] <- meta.factor
list.core[[label]] <- otu.fac.list
}
return(list.core)
}
core.Taxa<-function(data, is.OTU=FALSE, rank="g",
drop.unclassified=TRUE,
meta, meta.factor="", percent=1){
if ( class(data) != "list" ) {
stop("please provide otu tables as list. See ?RAM.input.formatting for details")
}
.valid.data(data=data, is.OTU=is.OTU)
# valide rank
if ( is.OTU ) {
rank.name <- .get.rank.name(rank, plural=TRUE)
rank_name <- .get.rank.name(rank)
rank_pat <- .get.rank.pat(rank)
} else {
warning("data are not otu tables, will ignore the rank provided")
rank_name <- "taxa"
rank.name <- "taxa"
rank_pat <- ""
}
list.core <-list()
labels <- names(data)
for ( i in 1:length(data) ) {
label <- names(data)[i]
elem <- data[[i]]
if ( is.null(elem) ) { break }
if ( is.OTU ) {
valid.OTU(elem)
.valid.meta(otu1=elem, meta=meta)
# get rank name & pattern
rank.name <- .get.rank.name(rank, plural=TRUE)
rank_name <- .get.rank.name(rank)
rank_pat <- .get.rank.pat(rank)
# tax.abund
tax.all <- tax.abund(elem, rank=rank, drop.unclassified=FALSE,
count=TRUE)
} else {
elem <- elem[match(rownames(meta), rownames(elem)), ]
if (identical(rownames(elem), rownames(meta))) {
tax.all <- elem
} else {
stop("Error: metadata and data don't have same subjects")
}
}
# check meta category
if( !length(meta.factor)==1 ) {
stop("Error: please provide one factor variable in metadata")
} else {
.valid.factor(meta, meta.factor)
fac.levels <- levels(factor(meta[[meta.factor]]))
}
if ( drop.unclassified ) {
# remove unclassified columns
# this selects all columns NOT containing in the blacklist
remove.pat <- paste(.blacklist(), collapse="|")
tax <- tax.all[ , !grepl(remove.pat, names(tax.all),
ignore.case=TRUE), drop=FALSE]
} else {
tax <- tax.all
}
###check for zero entries?
tax.fac.list=list()
tax.fac.id = list()
for ( j in 1:length(fac.levels) ) {
# each fac levels
meta.fac <- meta[grep(fac.levels[j], meta[[meta.factor]]), ]
tax.fac <- tax[match(rownames(meta.fac), rownames(tax)),]
tax.fac <- tax.fac[rowSums(tax.fac)>0, colSums(tax.fac)>0]
# presence or absence
tax.fac.pa <- vegan::decostand(tax.fac, "pa")
tax.all.fac <- tax.all[match(rownames(meta.fac), rownames(tax.all)),]
tax.all.fac <- tax.all.fac[rowSums(tax.all.fac)>0, colSums(tax.all.fac)>0]
# presence or absence
tax.all.fac.pa <- vegan::decostand(tax.all.fac, "pa")
# select taxa present in more than percent of samples in each category
sel <- which(colSums(tax.fac.pa) >= round(percent*nrow(tax.fac.pa)))
if (length(sel) != 0) {
tax.fac.ids <- as.character(colnames(tax.fac)[sel])
tax.fac.id[[j]] <- as.character(colnames(tax.fac)[sel])
tax.fac.ids.num <- length(tax.fac.ids)
tax.fac.ids.perc <- 100*sum(tax.fac[, c(tax.fac.ids)])/sum(tax.all.fac)
} else {
tax.fac.ids <- ""
tax.fac.id[[j]] <- ""
tax.fac.ids.num <- 0
tax.fac.ids.perc <- 0
}
# summarize
description <- paste(tax.fac.ids.num, "_", rank.name, "_found_in_",
percent*100, "%_", fac.levels[j], "_samples; ",
tax.fac.ids.perc, "%_of_total_sequences_in_",
fac.levels[j], sep="")
tax.fac.list[[j]]<- list(summary=description, taxa=tax.fac.ids)
names(tax.fac.list)[[j]] <- fac.levels[j]
}
# find taxa across all categories
tax.fac.list.num <- length(tax.fac.list)
tax.fac.id[[tax.fac.list.num+1]] <- as.character(factor(
Reduce(intersect, tax.fac.id)))
# taxa in all samples
# percent of sequences in core taxa across all categories
if( unique(tax.fac.id[[tax.fac.list.num+1]]) == "" ||
length(unique(tax.fac.id[[tax.fac.list.num+1]])) == 0 ) {
tax.fac.id.all.num <- 0
tax.fac.perc <- 0
tax.fac.id.all <- ""
} else {
tax.fac.id.all.num <- length(tax.fac.id[[tax.fac.list.num+1]])
tax.fac.id.all <- tax.fac.id[[tax.fac.list.num+1]]
tax.sel <- tax[, which(colnames(tax) %in% tax.fac.id.all)]
tax.fac.perc <- 100*sum(tax.sel)/sum(tax.all)
}
description <- paste(tax.fac.id.all.num, "_", rank.name, "_found_in_",
percent*100, "%_of_samples_at_each_", meta.factor,
"; ", tax.fac.perc, "%_of_total_sequences",
sep="")
tax.fac.list[[tax.fac.list.num+1]] <- list(summary=description, taxa=tax.fac.id.all)
names(tax.fac.list)[[tax.fac.list.num+1]] <- meta.factor
list.core[[label]] <- tax.fac.list
}
return(list.core)
}
group.OTU <- function(otu, rank="g", otuIDs="", meta,
meta.factor="", boxplot=TRUE,
main="", file=NULL, ext=NULL,
height=8, width=16) {
if(!length(meta.factor)==1) {
stop("Error: please provide one factor variable in metadata")
} else {
fac.levels <- levels(factor(meta[[meta.factor]]))
}
# get the groups
if(!length(meta.factor)==1) {
stop ("Error: please provide one factor")
}
if(length(otuIDs)==0) {
stop ("Error: please provide list of otuIDs to plot")
}
# validate inputs
valid.OTU(otu)
.valid.meta(otu1=otu, meta=meta)
# tax.split
if ( is.null(rank) ) {
otu <- LCA.OTU(otu, strip.format=FALSE, drop=TRUE)
rank_name <- "LCA"
} else {
otu <- tax.split(otu1=otu, rank=rank)
rank.name <- .get.rank.name(rank, plural=TRUE)
rank_name <- .get.rank.name(rank)
rank_pat <- .get.rank.pat(rank)
}
otu.p <- cbind(decostand(otu[, -ncol(otu)], "total",
MARGIN=2), otu[[rank_name]])
colnames(otu.p)[ncol(otu.p)] <- rank_name
otu.p <- otu.p[otu.p[[rank_name]] !="", ]
if( length(which(rownames(otu.p) %in% otuIDs))==0) {
stop ("Error: no selected otuIDs present in the otu table")
}
otuIDs <- unique(otuIDs)
otu.p.sel <- otu.p[which(rownames(otu.p) %in% c(otuIDs)), ]
otu.sel <- otu[which(rownames(otu) %in% c(otuIDs)), ]
otus.ex<-vector()
for (i in otuIDs) {
if(!(i %in% rownames(otu.p))) {
otus.ex<-unique(c(otus.ex, i))
}
}
if (length(otus.ex) != 0) {
print(paste("The following selected otus in your list were not in the otu table: ", paste(otus.ex, collapse=", "), sep=""))
}
# combine otu with taxonomy
rownames(otu.p.sel) <- paste(paste0("OTU_ID|",rownames(otu.p.sel)),
otu.p.sel[[rank_name]], sep=": ")
rownames(otu.sel) <- paste(paste0("OTU_ID|",rownames(otu.sel)),
otu.sel[[rank_name]], sep=": ")
# transpose OTU table
otu.p.sel.t <- as.data.frame(t(otu.p.sel[, -ncol(otu.p.sel)]))
otu.p.sel.tax <- otu.p.sel[[rank_name]]
otu.sel.t <- as.data.frame(t(otu.sel[, -ncol(otu.sel)]))
otu.sel.tax <- otu.sel[[rank_name]]
if(identical(rownames(otu.p.sel.t), rownames(meta))) {
otu.p.sel.fac <- cbind(otu.p.sel.t, meta[[meta.factor]])
names(otu.p.sel.fac)[ncol(otu.p.sel.fac)] <- meta.factor
otu.sel.fac <- cbind(otu.sel.t, meta[[meta.factor]])
names(otu.sel.fac)[ncol(otu.sel.fac)] <- meta.factor
} else {
stop("samples not identical in otu and metadata")
}
# melt by sample
#if (!require("reshape2")) {
# stop("package 'reshape2' is required to use this function")
# }
otu.sel.fac.m <- reshape2::melt(cbind(otu.sel.fac, Sample=rownames(otu.sel.fac)),
tax=otu.sel.fac$tax, variable.name="OTU",
value.name="Count")
names(otu.sel.fac.m)[ncol(otu.sel.fac.m)-1] <- "OTU"
names(otu.sel.fac.m)[ncol(otu.sel.fac.m)] <- "Count"
# boxplot or barplot
if(boxplot) {
# boxplot
otu.p.sel.fac.m <- reshape2::melt(cbind(otu.p.sel.fac,
Sample=rownames(otu.p.sel.fac)),
variable.name="OTU", value.name="RA")
names(otu.p.sel.fac.m)[ncol(otu.p.sel.fac.m)-1] <- "OTU"
names(otu.p.sel.fac.m)[ncol(otu.p.sel.fac.m)] <- "RA"
otu.p.sel.fac.m <- cbind(otu.p.sel.fac.m, do.call(rbind,
strsplit(x=as.character(otu.p.sel.fac.m$OTU),split=":")))
names(otu.p.sel.fac.m)[ncol(otu.p.sel.fac.m)-1] <- "otu"
names(otu.p.sel.fac.m)[ncol(otu.p.sel.fac.m)] <- "tax"
otu.p.sel.fac.m <- otu.p.sel.fac.m[order(otu.p.sel.fac.m$tax,
otu.p.sel.fac.m$otu), ]
} else {
# barplot
otu.sel.fac.agg = stats::aggregate(otu.sel.t,
by=list(meta[[meta.factor]]), FUN=sum)
otu.sel.fac.agg <- cbind(otu.sel.fac.agg[,1],
vegan::decostand(otu.sel.fac.agg[, -1],
MARGIN=2, "total"))
names(otu.sel.fac.agg)[1] <- meta.factor
# label counts as "RA" to be consistent with previous df
otu.sel.fac.agg.m <- reshape2::melt(otu.sel.fac.agg,
is.vars=c(meta.factor), value.name="RA",
variable.name="OTU")
names(otu.sel.fac.agg.m)[ncol(otu.sel.fac.agg.m)-1] <- "OTU"
names(otu.sel.fac.agg.m)[ncol(otu.sel.fac.agg.m)] <- "RA"
otu.sel.fac.agg.m <- cbind(otu.sel.fac.agg.m,
do.call(rbind,strsplit(x=as.character(otu.sel.fac.agg.m$OTU),
split=":")))
names(otu.sel.fac.agg.m)[ncol(otu.sel.fac.agg.m)-1] <- "otu"
names(otu.sel.fac.agg.m)[ncol(otu.sel.fac.agg.m)] <- "tax"
otu.sel.fac.agg.m <- otu.sel.fac.agg.m[order(otu.sel.fac.agg.m$tax,
otu.sel.fac.agg.m$otu), ]
otu.p.sel.fac.m <- otu.sel.fac.agg.m
#return(otu.p.sel.fac.m)
}
#calculate total count of each Taxon
total= sapply(levels(otu.sel.fac.m$OTU),
function(x){sum(otu.sel.fac.m$Count[otu.sel.fac.m$OTU==x])},
USE.NAMES=F)
total.count <- data.frame(OTU=levels(otu.sel.fac.m$OTU),total)
# split otu&tax
total.count<-cbind(total.count,
do.call(rbind,strsplit(x=as.character(total.count$OTU),
split=":")))
names(total.count)[ncol(total.count)-1] <- "otu"
names(total.count)[ncol(total.count)] <- "tax"
#total.count <- total.count[order(total.count$tax,
#total.count$OTU, total.count$total, decreasing=TRUE), ]
#return(list(otu.sel.fac.m, otu.p.sel.fac.m, total.count))
# merge total.count and otu.p.sel.fac.m
otu.p.sel.fac.m.m <- merge(otu.p.sel.fac.m, total.count[, c("OTU", "total")], by="OTU", all=TRUE)
# reorder levels based on order of appearance in total.count,
#otu.p.sel.fac.m$OTU <- factor(as.character(otu.p.sel.fac.m$OTU),
# levels=total.count$OTU, ordered=TRUE)
# we need to use aes_string to pass CRAN check; see
# http://goo.gl/JxgZ9u
# the distribution of each taxon in each category of meta.factor
if (!requireNamespace("scales")) {
stop("package 'scales' is required to use this function")
}
if (!requireNamespace("RColorBrewer")) {
stop("package 'RColorBrewer' is required to use this function")
}
if (!requireNamespace("grid")) {
stop("package 'grid' is required to use this function")
}
ylab <- "Relative Abundance (%)"
# sort by meta.fac and RA
aa<-otu.p.sel.fac.m.m
#return(aa)
#aa<-aa[order(aa[, meta.factor], aa[, "RA"]),]
#aa<-aa[order(aa[, meta.factor], aa[, "RA"], decreasing=TRUE),]
aa<-aa[order(aa[, meta.factor], aa[, "RA"], decreasing=FALSE),]
# sort each meta.factor level by RA
#for ( i in levels(aa[, meta.factor])) {
# bb <-aa[aa[, meta.factor]==i, ]
#bb <- bb[order(bb[, "RA"], decreasing=TRUE),]
# bb <- bb[order(bb[, "RA"]),]
# aa[aa[, meta.factor]==i, ] <- bb
#}
#return(aa)
# re-order meta.factor levels
#aa[, meta.factor] <- factor(as.character(aa[, meta.factor] ),
# levels= unique(sort(aa[, meta.factor],
# decreasing=FALSE)), ordered=TRUE)
#Turn your 'OTU' column into a character vector
#aa$OTU <- as.character(aa$OTU)
#Then turn it back into an ordered factor
#aa$OTU <- factor(aa$OTU, levels=unique(aa$OTU), ordered=TRUE)
otu.p.sel.fac.m <-aa
otu.p.sel.fac.m$OTU <- factor(as.character(otu.p.sel.fac.m$OTU),
levels= unique(otu.p.sel.fac.m$OTU), ordered=TRUE)
#return(otu.p.sel.fac.m)
if(boxplot) {
p <- ggplot(otu.p.sel.fac.m, aes_string(x="OTU", y="RA",
fill=meta.factor)) +
geom_boxplot()
} else {
p <- ggplot(otu.p.sel.fac.m, aes_string(x="OTU", y="RA",
fill=meta.factor)) +
geom_bar(colour="white",stat="identity",
position="fill")
}
cols.needed <- length(levels(factor(meta[[meta.factor]])))
if (cols.needed <= 12 ) {
p <- p + scale_fill_brewer(type="div", palette="Set3")
} else {
p <- p + scale_fill_manual(values=RAM.pal(cols.needed))
}
p <- p + scale_x_discrete(labels = paste(otu.p.sel.fac.m$otu, "\n",
otu.p.sel.fac.m$tax,":",otu.p.sel.fac.m$total, sep="")) +
theme(legend.key.size=unit(6,"mm"),
legend.text = element_text(size=10),
axis.text.y=element_text(size=7),
legend.position="right") +
xlab("OTU") +
ylab(ylab) +
ggtitle(main) +
coord_flip()
p <- p + scale_y_continuous(labels=scales::percent_format())
.valid.plot.settings(file, ext)
save <- !is.null(file)
if (save) {
.ggsave.helper(file, ext, width, height, plot=p)
} else {
p
}
}
# venn for classified taxa
group.venn <- function(vectors, cat.cex=1.5, cex=1,
cat.pos=NULL, cat.dist=NULL,
label=TRUE, lab.cex=1,
lab.col= "black", fill=NULL,
file=NULL, ext=NULL, width=8, height=8) {
save <- !is.null(file)
if (save) { .get.dev(file, ext, height=height, width=width) }
if ( !requireNamespace("VennDiagram") ) {
stop("package 'VennDiagram' is required for this function")
}
if ( !requireNamespace("RColorBrewer") ) {
stop("package 'RColorBrewer' is required for this function")
}
if (!requireNamespace("grid")) {
stop("package 'grid' is required to use this function")
}
# Generate plot
# number of vectors to plot
len <- length(vectors)
if ( is.null(fill) ) {
if ( len == 2 ) {
fill = c("lightpink", "lightblue")
} else {
fill = RColorBrewer::brewer.pal(len, "Pastel1")
}
} else {
if ( length(fill) == len ) {
fill = fill
} else if ( length(fill) > len ) {
warning(paste("more colors being provided than required, will ignore ", length(fill)-len, " colors", sep=""))
fill = fill[1:len]
} else {
warning("not enough colors being provided, will use default")
if ( len == 2 ) {
fill = c("lightpink", "lightblue")
} else {
fill = RColorBrewer::brewer.pal(len, "Pastel1")
}
}
}
if ( len > 2 && label ) {
warning("currently only support 2 groups to have actual item labels; will only use numbers")
} else if ( len > 5 || len < 2 ) {
stop("please provide 2 to 5 vectors")
}
alpha = rep(0.5, len)
if ( !is.null(cat.pos) && !is.null(cat.dist) ) {
v <- VennDiagram::venn.diagram(vectors, fill = fill, alpha = alpha,
cat.dist=cat.dist, cat.pos=cat.pos,
cat.fontface = "bold", cat.cex = cat.cex, cex=cex,
filename=NULL)
} else if ( !is.null(cat.pos) && is.null(cat.dist) ) {
v <- VennDiagram::venn.diagram(vectors, fill = fill, alpha = alpha,
cat.pos=cat.pos, cat.fontface = "bold",
cat.cex = cat.cex, cex=cex,
filename=NULL)
} else if ( is.null(cat.pos) && !is.null(cat.dist) ) {
v <- VennDiagram::venn.diagram(vectors, fill = fill, alpha = alpha,
cat.fontface = "bold", cat.dist= cat.dist,
cat.cex = cat.cex, cex=cex,
filename=NULL)
} else {
v <- VennDiagram::venn.diagram(vectors, fill = fill, alpha = alpha,
cat.fontface = "bold", cat.cex = cat.cex, cex=cex,
filename=NULL)
}
if ( len > 2 && len <=5 ) {
grid::grid.newpage()
grid::grid.draw(v)
}
if ( len==2 ) {
if ( !label ) {
grid::grid.newpage()
grid::grid.draw(v)
} else {
name <- lapply(v, names)
lapply(v, function(i) i$label)
# plot with labels
# first find out whether the vectors got rearranged
v.labels <- lapply(v, function(i) i$label)
v.lab <- vector()
for ( i in 1:length(v.labels) ) {
if ( length(v.labels[[i]] %in% names(vectors)) !=0 &&
isTRUE(v.labels[[i]] %in% names(vectors)) ) {
v.lab <- c(v.lab, v.labels[[i]])
}
}
v1 <- vectors[[v.lab[1]]]
v2 <- vectors[[v.lab[2]]]
v[[5]]$label <- paste(c(v[[5]]$label,setdiff(v1, v2)), collapse="\n")
v[[5]]$gp$cex <- lab.cex
v[[5]]$gp$col <- lab.col
# in baa only
v[[6]]$label <- paste(c(v[[6]]$label,setdiff(v2, v1)) , collapse="\n")
v[[6]]$gp$cex <- lab.cex
v[[6]]$gp$col <- lab.col
# intesection
v[[7]]$label <- paste(c(v[[7]]$label,intersect(v1, v2)), collapse="\n")
v[[7]]$gp$cex <- lab.cex
# plot with labels
v[[7]]$gp$col <- lab.col
# plot with labels
grid::grid.newpage()
grid::grid.draw(v)
}
}
if (save) { dev.off() }
invisible()
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.