Commit 7b1157b1 authored by Sebastien Theil's avatar Sebastien Theil
Browse files

Merge branch 'master' of forgemia.inra.fr:umrf/ranomaly

parents 3fb77eb3 c439e564
Pipeline #44135 passed with stage
in 16 seconds
.Rproj.user
.Rhistory
README.Rmd
renv
......@@ -3,6 +3,7 @@
export(ASVenn_fun)
export(VENNFUN)
export(aggregate_fun)
export(aggregate_top_taxa)
export(assign_taxo_fun)
export(bars_fun)
export(check_tax_fun)
......@@ -61,9 +62,13 @@ importFrom(DECIPHER,AlignSeqs)
importFrom(DESeq2,DESeq)
importFrom(DESeq2,results)
importFrom(DESeq2,resultsNames)
importFrom(dplyr,`%>%`)
importFrom(dplyr,across)
importFrom(dplyr,arrange)
importFrom(dplyr,group_by)
importFrom(dplyr,group_map)
importFrom(dplyr,left_join)
importFrom(dplyr,mutate)
importFrom(ggpubr,ggtexttable)
importFrom(ggpubr,text_grob)
importFrom(ggpubr,ttheme)
......@@ -83,16 +88,21 @@ importFrom(metacoder,heat_tree)
importFrom(metacoder,parse_phyloseq)
importFrom(metacoder,taxon_names)
importFrom(metacoder,zero_low_counts)
importFrom(microbiome,aggregate_top_taxa)
importFrom(microbiome,aggregate_taxa)
importFrom(microbiome,core_members)
importFrom(microbiome,top_taxa)
importFrom(microbiome,transform)
importFrom(nlme,lme)
importFrom(phangorn,NJ)
importFrom(phangorn,dist.ml)
importFrom(phangorn,phyDat)
importFrom(phangorn,pml)
importFrom(plotly,config)
importFrom(plotly,ggplotly)
importFrom(plotly,layout)
importFrom(plotly,plot_ly)
importFrom(plotly,subplot)
importFrom(ranacapa,ggrare)
importFrom(reshape2,melt)
importFrom(tibble,rownames_to_column)
importFrom(venn,venn)
......@@ -10,6 +10,8 @@
#' @param lvls Vector comma separated list levels of factor to print in venn diagram (max. 5).
#' @param krona Krona of exclusive ASV or shared with informed level and others. Must be among levels of column1 argument.
#' @param shared shared [TRUE] or exclusive [FALSE] mode.
#' @param verbose Verbose level. (1: quiet, 2: print infos, 3: print infos + debug)
#' @param ggplotmode if TRUE plot the Venn diagram using ggplot
#'
#'
#' @return Returns list with venn diagram and table with shared features. Exports a venn diagram with corresponding tabulated file.
......@@ -22,10 +24,22 @@
ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
column1 = NULL, subset = "", lvls = NULL, krona = "",
shared = TRUE){
shared = TRUE, verbose = 2, ggplotmode = FALSE){
invisible(flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger"))
if(verbose == 3){
invisible(flog.threshold(DEBUG))
}
if(verbose == 2){
invisible(flog.threshold(INFO))
}
if(verbose == 1){
invisible(flog.threshold(ERROR))
}
if(!is.null(output)){
if(!dir.exists(output)){
dir.create(output, recursive=TRUE)
......@@ -68,8 +82,8 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
if(length(lvls)==0){
flog.error('You must provide levels...')
stop()
} else if(length(lvls)>5){
flog.error('Venn diagram is limited to 5 levels')
} else if(length(lvls)>7){
flog.error('Venn diagram is limited to 7 levels')
stop()
} else{
if(!all(lvls %in% na.omit(levels(as.factor(sample_data(data)[,column1]@.Data[[1]])) ))){
......@@ -79,19 +93,6 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
}
if(length(lvls)==0){
flog.error('You must provide levels...')
stop()
} else if(length(lvls)>5){
flog.error('Venn diagram is limited to 5 levels')
stop()
} else{
if(!all(lvls %in% na.omit(levels(as.factor(sample_data(data)[,column1]@.Data[[1]])) ))){
flog.error('Your levels are not present in metadata...')
stop()
}
}
#Nombre d'espèce par matrice
flog.info('Parsing factor ...')
......@@ -111,14 +112,14 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
for(i in 1:length(lvls)){
databak -> data
LOC=as.character(lvls[i])
print(LOC)
flog.info(LOC)
fun <- paste("data <- subset_samples(data, ",column1," %in% '",LOC,"')",sep="")
eval(parse(text=fun))
# print(rank)
sp_data <- prune_taxa(taxa_sums(data) > 0, data)
# cat(LOC,ntaxa(sp_data)," ", rank, " \n")
print(data)
# print(data)
ttable <- sp_data@tax_table@.Data
otable <- as.data.frame(otu_table(sp_data))
# print(nrow(ttable))
......@@ -213,7 +214,7 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
}
}
} else {
print("< 5 levels")
flog.info("< 5 levels")
if(!is.null(lvls)){
flog.info(glue('Selecting {lvls} ...'))
LVLs <- unlist(strsplit(lvls,","))
......@@ -238,6 +239,7 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
#' @param output Output path.
#' @param refseq1 Reference sequences.
#' @param alltax Taxonomy table.
#' @param ggplotmode if TRUE plot the Venn diagram using ggplot
#'
#'
#' @return Exports a venn diagram with corresponding tabulated file.
......@@ -247,20 +249,26 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
#' @export
VENNFUN <- function(TF = TF, mode = 1, TITRE = TITRE, output = "./", refseq1 = NULL, alltax=NULL){
VENNFUN <- function(TF = TF, mode = 1, TITRE = TITRE, output = "./", refseq1 = NULL, alltax=NULL, ggplotmode = FALSE){
if(mode==1){
# pdf(NULL)
venn::venn(TF, zcol = rainbow(7), ilcs = 2, sncs = 2, ggplot = FALSE) #, col=rainbow(7)
flog.info("VennDiagram")
venn::venn(TF, zcol = rainbow(7), ilcs = 1.5, sncs = 2, ggplot = ggplotmode) #, col=rainbow(7)
venn.plot <- recordPlot()
invisible(dev.off())
flog.info("Output figure")
if(!is.null(output)){
png(paste(output,'/',TITRE,'_venndiag.png',sep=''), width=20, height=20, units="cm", res=200)
replayPlot(venn.plot)
dev.off()
}
print("plotOK")
flog.info("plotOK")
ov <- calculate.overlap(TF)
......@@ -281,6 +289,7 @@ VENNFUN <- function(TF = TF, mode = 1, TITRE = TITRE, output = "./", refseq1 = N
names(TABf)[1] = names(TF)[j]
}
flog.info('Generating table ...')
if(!is.null(alltax)){
if(!is.null(refseq1)){
TABf <- cbind(TABf,alltax[as.character(TABf$TABf),2], refseq1[as.character(TABf$TABf),])
......@@ -291,6 +300,7 @@ VENNFUN <- function(TF = TF, mode = 1, TITRE = TITRE, output = "./", refseq1 = N
}
}
flog.info('Writing table ...')
if(!is.null(output)){
write.table(TABf, paste(output,"/",TITRE,"_venn_table.csv",sep=""), sep="\t", quote=FALSE, row.names=FALSE)
}
......@@ -299,6 +309,7 @@ VENNFUN <- function(TF = TF, mode = 1, TITRE = TITRE, output = "./", refseq1 = N
invisible(dev.off())
venn.plot <- venn::venn(TF, zcol = rainbow(7), ilcs = 2, sncs = 2, ggplot = FALSE) #, col=rainbow(7)
venn.plot <- recordPlot()
invisible(dev.off())
if(!is.null(output)){
......@@ -337,7 +348,9 @@ VENNFUN <- function(TF = TF, mode = 1, TITRE = TITRE, output = "./", refseq1 = N
}
LL=list()
LL$venn_plot = venn.plot
LL$TABf = TABf
LL$venn_plot <- venn.plot
LL$TABf <- TABf
flog.info('Finish ...')
return(LL)
}
......@@ -9,7 +9,7 @@
#' @param mgseq Path to metagenomeseq results folder
#' @param column1 Column name of factor to test (among sample_variables(data))
#' @param column2 Column name on which table were splitted
#' @param verbose Verbose level. (1: quiet, 3: verbal)
#' @param verbose Verbose level. (1: quiet, 2: print infos, 3: print infos + debug)
#' @param rank Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)
#' @param comp Comparison to test. Comma separated and comparisons are informed with a tilde (A~C,A~B,B~C). If empty, test all combination.
#' @param returnval Boolean for function to return values.
......@@ -34,16 +34,21 @@ aggregate_fun <- function(data = data, metacoder = NULL, deseq = NULL, mgseq = N
if(verbose == 3){
invisible(flog.threshold(DEBUG))
} else {
}
if(verbose == 2){
invisible(flog.threshold(INFO))
}
if(verbose == 1){
invisible(flog.threshold(ERROR))
}
if(!dir.exists(output)){
dir.create(output,recursive=T )
}
flog.info('Done.')
# data.glom <- tax_glom(data, taxrank=rank, h=TRUE)
flog.info('Glom Tax.')
data.glom <- tax_glom(data, taxrank=rank)
if(comp == ''){
fun <- paste('combinaisons <- combn(na.omit(unique(sample_data(data)$',column1,')),2)',sep='')
......@@ -79,13 +84,12 @@ aggregate_fun <- function(data = data, metacoder = NULL, deseq = NULL, mgseq = N
flog.warn('File does not exists.')
deseqT <- data.frame()
}
# print(head(deseqT))
flog.info('MetagenomeSeq.')
if(file.exists(paste(mgseq,"/signtab_",column1,"_",combinaisons[1,col],"_vs_",combinaisons[2,col],".csv",sep=""))){
mgseqT <- read.table(paste(mgseq,"/signtab_",column1,"_",combinaisons[1,col],"_vs_",combinaisons[2,col],".csv",sep=""), h=TRUE,sep="\t")
} else{
if(file.exists(paste(mgseq,"/signtab_",column1,"_",combinaisons[2,col],"_vs_",combinaisons[1,col],".csv",sep=""))){
mgseqT <- read.table(paste(mgseq,"/signtab_",column1,"_",combinaisons[1,col],"_vs_",combinaisons[2,col],".csv",sep=""), h=TRUE,sep="\t")
mgseqT <- read.table(paste(mgseq,"/signtab_",column1,"_",combinaisons[2,col],"_vs_",combinaisons[1,col],".csv",sep=""), h=TRUE,sep="\t")
}
else{
flog.warn('File does not exists.')
......@@ -172,7 +176,7 @@ aggregate_fun <- function(data = data, metacoder = NULL, deseq = NULL, mgseq = N
names(TABf)[ncol(TABf)] = names(TF)[j]
}
TABfbak <- TABf
TABfbak0 <- TABf
# add new columns, sumMethods, DeseqLFC, Mean Relative Abundance (TSS) condition 1 & 2
row.names(deseqT) = deseqT[,1]
......@@ -188,7 +192,7 @@ aggregate_fun <- function(data = data, metacoder = NULL, deseq = NULL, mgseq = N
# clr = function(x){log(x+1) - rowMeans(log(x+1))}
# otableNORM <- clr(otable)
normf = function(x){ x/sum(x) }
data.norm <- transform_sample_counts(data, normf)
data.norm <- transform_sample_counts(data.glom, normf)
otableNORM <- otu_table(data.norm)
Gtab <- cbind(as.data.frame(ssample), t(otableNORM))
......@@ -270,10 +274,9 @@ aggregate_fun <- function(data = data, metacoder = NULL, deseq = NULL, mgseq = N
# print(TABfinal)
if(length(TABfinal) > 0){
if(length(TABfinal$ListAllOtu) > 0){
flog.info('Building csv file...')
TAX <- cbind(seqid = rownames(ttax[as.character(TABfinal$ListAllOtu),]),as.data.frame(ttax[as.character(TABfinal$ListAllOtu),]))
sseq <- cbind(seqid = row.names(as.data.frame(sseq)), as.data.frame(sseq))
......@@ -324,7 +327,7 @@ aggregate_fun <- function(data = data, metacoder = NULL, deseq = NULL, mgseq = N
# save(data, file=paste(output,"/signif_phyloseq.rdata",sep=""))
#
}
}else{flog.info("No significant results in the three methods...")}
}
......@@ -31,6 +31,32 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){
}
#' aggregate_top_taxa from microbiome package
#'
#'
#' @param x phyloseq object
#' @param top Keep the top-n taxa, and merge the rest under the category 'Other'. Instead of top-n numeric this can also be a character vector listing the groups to combine.
#' @param level Summarization level (from ‘rank_names(pseq)’)
#'
#' @importFrom microbiome aggregate_taxa
#' @importFrom microbiome top_taxa
#'
#' @export
aggregate_top_taxa <- function (x, top, level){
x <- aggregate_taxa(x, level)
tops <- top_taxa(x, top)
tax <- tax_table(x)
inds <- which(!rownames(tax) %in% tops)
tax[inds, level] <- "Other"
tax_table(x) <- tax
tt <- tax_table(x)[, level]
tax_table(x) <- tax_table(tt)
aggregate_taxa(x, level)
}
#' Barplots plotly
#'
#'
......@@ -38,39 +64,44 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){
#' @param rank Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)
#' @param top Number of top taxa to plot
#' @param Ord1 Variable used to order sample (X axis) or split the barplot if split = TRUE
#' @param Fact1 Variable used to change X axis tick labels and color (when split = FALSE)
#' @param sample_labels If true, x axis labels are sample IDS, if false labels displayed are levels from Ord1 argument. Ignored if split = TRUE (FALSE)
#' @param split if TRUE make a facet_wrap like grouped by Ord1 (default FALSE)
#' @param relative Plot relative (TRUE, default) or raw abundance plot (FALSE)
#' @param autoorder Automatic ordering xaxis labels based on Ord1 factor levels with gtools::mixedorder function (TRUE).
#' @param ylab Y axis title ("Abundance")
#' @param outfile Output html file.
#'
#' @return Returns barplots in an interactive plotly community plot
#'
#' @importFrom plotly plot_ly subplot layout
#' @importFrom plotly plot_ly subplot layout ggplotly
#' @importFrom microbiome aggregate_top_taxa
#' @importFrom reshape2 melt
#' @importFrom gtools mixedsort
#' @importFrom dplyr group_map group_by across
#' @importFrom dplyr group_map group_by across `%>%` mutate arrange
#'
#'
#' @export
bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 = NULL, split = FALSE,
relative = TRUE,
outfile="plot_compo.html"){
bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, sample_labels = FALSE, split = FALSE, split_sid_order = FALSE,
relative = TRUE, autoorder = TRUE, ylab = "Abundance", outfile="plot_compo.html", verbose = TRUE){
# Fdata <- prune_samples(sample_names(r$data16S())[r$rowselect()], r$data16S())
# Fdata <- prune_taxa(taxa_sums(Fdata) > 0, Fdata)
# if( r$RankGlom() == "ASV"){
# Fdata <- prune_taxa(r$asvselect(), Fdata)
# }
if(verbose){
invisible(flog.threshold(INFO))
} else {
invisible(flog.threshold(ERROR))
}
if( all(Ord1 != sample_variables(data))){
stop(paste("Wrong value in Ord1, please use variables existing in the phyloseq object:", toString(sample_variables(data))))
}
flog.info('Preprocess...')
Fdata = data
# print("top")
psobj.top <- aggregate_top_taxa(Fdata, rank, top = top)
# print("get data")
sdata = as.data.frame(sample_data(psobj.top))
sdata <- as.data.frame(sample_data(psobj.top), stringsAsFactors = TRUE)
sdata$sample.id = sample_names(psobj.top)
otable = as.data.frame(otu_table(psobj.top))
row.names(otable) = tax_table(psobj.top)[,rank]
......@@ -78,6 +109,8 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 =
# print("melt data")
dat <- as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
flog.info(' Melting table...')
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
meltdat$variable <- factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
......@@ -87,23 +120,51 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 =
# print(levels(meltdat$sample.id))
# save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
# TODO: Message d'erreur si factor n'est pas dans les sample_data
fun = glue( "xform <- list(categoryorder = 'array',
categoryarray = unique(meltdat$sample.id[gtools::mixedorder(meltdat${Ord1})]),
title = 'Samples',
tickmode = 'array',
tickvals = 0:nrow(sdata),
ticktext = sdata[unique(meltdat$sample.id[gtools::mixedorder(meltdat${Ord1})]), '{Fact1}']@.Data[[1]],
tickangle = -90)")
eval(parse(text=fun))
if(autoorder){
flog.info(' Ordering samples...')
fun = glue( "labs = gtools::mixedorder(as.character(meltdat${Ord1}))" )
eval(parse(text=fun))
orderedIDS <- unique(meltdat$sample.id[gtools::mixedorder(as.character(meltdat[,Ord1]))])
orderedOrd1 <- meltdat[,Ord1][gtools::mixedorder(as.character(meltdat[,Ord1]))]
orderedOrd1 <- factor(orderedOrd1, levels = gtools::mixedsort(levels(orderedOrd1)))
}else{
labs = 1:nrow(meltdat)
orderedIDS <- unique(meltdat$sample.id)
orderedOrd1 <- meltdat[,Ord1]
}
if(sample_labels){
lab1 = "sample.id"
}else{
lab1 = Ord1
}
flog.info(' Set labels...')
xform <- list(categoryorder = 'array',
categoryarray = unique(meltdat$sample.id[labs]),
title = 'Samples',
tickmode = 'array',
tickvals = 0:nrow(sdata),
ticktext = sdata[as.character(unique(meltdat$sample.id[labs])), lab1]@.Data[[1]],
tickangle = -90)
# subplot to vizualize groups
# print(head(sdata))
df1 <- cbind.data.frame(x=sdata[unique(meltdat$sample.id[gtools::mixedorder(meltdat[,Ord1])]), "sample.id"]@.Data[[1]],
g=sdata[unique(meltdat$sample.id[gtools::mixedorder(meltdat[,Ord1])]), Fact1]@.Data[[1]],
flog.info(' Subplot...')
df1 <- cbind.data.frame(x=sdata[orderedIDS, "sample.id"]@.Data[[1]],
g=sdata[orderedIDS, Ord1]@.Data[[1]],
y=1)
fun = glue( "df1$g <- factor(df1$g, levels = as.character(unique(orderedOrd1)))")
eval(parse(text=fun))
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(unique(orderedOrd1)))") # keep sample ids order from metadata.
eval(parse(text=fun))
subp1 <- df1 %>% plot_ly(
type = 'bar',
x = ~x,
......@@ -116,14 +177,18 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 =
if(relative){
flog.info('Plotting relative...')
#relative abondance
plottitle = "Relative abundance"
otable=apply(otable,2, function(x){Tot=sum(x); x/Tot})
dat= as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
meltdat = reshape2::melt(dat, id.vars=1:ncol(sdata))
tt=levels(meltdat$variable)
meltdat$variable = factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
meltdat$variable <- factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(unique(orderedOrd1)))")
eval(parse(text=fun))
p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable
plotly::layout(title="Relative abundance", yaxis = list(title = 'Relative abundance'), xaxis = xform, barmode = 'stack')
......@@ -133,8 +198,8 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 =
plotly::layout(xaxis = xform)
}
}else{
flog.info('Plotting raw...')
#raw abundance
plottitle = "Raw abundance"
p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable
plotly::layout(title="Raw abundance", yaxis = list(title = 'Raw abundance'), xaxis = xform, barmode = 'stack')
......@@ -143,32 +208,40 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 =
}
}
# dir.create(outpath, recursive = TRUE)
# htmlwidgets::saveWidget(p1, glue::glue("{outpath}/{outfile}"))
if(!is.null(outfile)){
htmlwidgets::saveWidget(p1, outfile)
}
# facet_wrap output
# Splitted plot output
if(!split) {
if(!is.null(outfile)){
htmlwidgets::saveWidget(p1, outfile)
}
flog.info('Finish...')
return(p1)
} else {
p1 = meltdat %>% group_by(across({Ord1})) %>%
dplyr::group_map(~ plot_ly(data=., x = ~sample.id, y = ~value, type = 'bar',
name = ~variable,
color = ~variable, legendgroup = ~variable,
showlegend = (.y == levels(meltdat[, Ord1])[1])),
keep = TRUE) %>%
flog.info('Splitted plot...')
if(split_sid_order){
meltdat$sample.id = factor(meltdat$sample.id, levels = unique(meltdat$sample.id)) # keep sample ids order from metadata.
}
p1 = meltdat %>% dplyr::arrange(across({Ord1})) %>% dplyr::group_by(across({Ord1})) %>%
dplyr::mutate(across(where(is.character), as.factor)) %>%
dplyr::group_map(~ plot_ly(data=., x = ~sample.id, y = ~value, type = 'bar',
name = ~variable,
color = ~variable, legendgroup = ~variable,
showlegend = (.y == levels(meltdat[, Ord1])[1])),
keep = TRUE) %>%
plotly::subplot(nrows = 1, shareX = TRUE, shareY=TRUE, titleX = FALSE) %>%
plotly::layout(title=plottitle,
xaxis = list(title = glue("{Ord1} = {unique(meltdat[, Ord1])[1]}")),
yaxis = list(title = 'Relative abundance'),
plotly::layout(title="",
xaxis = list(title = glue("{Ord1} =\n{levels(meltdat[, Ord1])[1]}")),
yaxis = list(title = ylab),
barmode = 'stack')
for (i in 2:length(unique(meltdat[, Ord1]))) {
p1$x$layoutAttrs[[1]][[paste0("xaxis", i)]] = NULL
p1$x$layoutAttrs[[1]][[paste0("xaxis", i)]]$title <- glue("{Ord1} = {unique(meltdat[, Ord1])[i]}")
p1$x$layoutAttrs[[1]][[paste0("xaxis", i)]]$title <- glue("{Ord1} =\n{levels(meltdat[, Ord1])[i]}")
}
if(!is.null(outfile)){
htmlwidgets::saveWidget(p1, outfile)
}
flog.info('Finish...')
return(p1)
}
}
#' Create Phyloseq object from external datas.
#'
#'
#' @param otutable Tabulated ASVtable file path.
#' @param taxtable Tabulated taxonomy table file path.
#' @param seq Tabulated sequence file path.
#' @param metadata Tabulated metadata file path.
#' @param otutable Path to the tabulated ASV/OTU table or data.frame object.
#' @param taxtable Path to the tabulated taxonomy table or data.frame object.
#' @param Path to sequence in fasta format or DNAStringSet object.
#' @param metadata Path to the tabulated metadata table or data.frame object.
#' @param generateTree Boolean to generate the phylogenetic tree.
#' @param output Output directory
#' @param returnval Boolean to return values in console or not.
......@@ -19,18 +19,41 @@ csv2phyloseq_fun <- function(otutable = NULL, taxtable = NULL, seq = NULL, metad
flog.info("Input...")
flog.info("...OTUtable...")
otable <- read.table(otutable, sep="\t", h=TRUE)
row.names(otable)=otable[,1]
if(is.data.frame(otutable)){
otable <- otutable
}else{
otable <- read.table(otutable, sep="\t", h=TRUE, check.names=FALSE)
row.names(otable)=otable[,1]
otable = otable[,-1]
}
flog.info("...TAXtable...")
ttable <- as.matrix(read.table(taxtable, sep="\t", h=TRUE, stringsAsFactors=FALSE))
row.names(ttable)=ttable[,1]
if(is.data.frame(taxtable)){
print("dF")
head(taxtable)
ttable <- taxtable
}else{
ttable <- as.matrix(read.table(taxtable, sep="\t", h=TRUE, stringsAsFactors=FALSE))
row.names(ttable)=ttable[,1] # first columns is otu ids.
ttable = ttable[,-1]
}
flog.info("...Metadatas...")
mtable <- read.table(metadata, sep="\t", h=TRUE)
if(is.data.frame(metadata)){
mtable <- metadata
}else{
mtable <- read.table(metadata, sep="\t", h=TRUE)
}
row.names(mtable)=mtable[,1]
flog.info("...Sequences...")
stable <- read.table(seq, sep="\t", h=TRUE)
sequences1 = DNAStringSet(stable[,-1])
names(sequences1)=stable[,1]
if(class(seq) == "DNAStringSet"){
sequences1 = seq
} else {
sequences1 <- readDNAStringSet(seq)