Mercurial > repos > simon-gladman > phyloseq_filter
diff phyloseq_filter.R @ 1:9fbb104e16d9 draft
update
| author | simon-gladman |
|---|---|
| date | Mon, 12 Nov 2018 22:42:14 -0500 |
| parents | 910739de7b80 |
| children |
line wrap: on
line diff
--- a/phyloseq_filter.R Wed Sep 05 22:21:43 2018 -0400 +++ b/phyloseq_filter.R Mon Nov 12 22:42:14 2018 -0500 @@ -2,10 +2,10 @@ library('ape') library('ggplot2') suppressPackageStartupMessages(library('phyloseq')) -library(biomformat) library(plyr) Sys.setenv("DISPLAY"=":1") library(biomformat) +library(jsonlite) suppressPackageStartupMessages(library(metagenomeSeq)) suppressPackageStartupMessages(library("doParallel")) ncores = ceiling(detectCores() * 0.8) @@ -65,6 +65,12 @@ pngfile_post_phyla_filtering<- gsub("[ ]+", "", paste(options$outdir,"/barplot_after_phyla_filtering.png")) htmlfile <- gsub("[ ]+", "", paste(options$htmlfile)) +### overwrite the write_biom function for proper BIOM format +### https://github.com/smdabdoub/biomformat/blob/master/R/IO-methods.R#L124 +write_biom <- function(x, biom_file){ + cat(toJSON(x, always_decimal=TRUE, auto_unbox=TRUE), file=biom_file) +} + ### This function accepts different two different type of BIOM file format readBIOM<-function(inBiom){ tryCatch({ @@ -152,7 +158,7 @@ #png('barplot_pre_phyla_filtering.png') bitmap(pngfile_pre_phyla_filtering,"png16m") - print(sample_data(physeq_pre_phyla_filtering)) + #print(sample_data(physeq_pre_phyla_filtering)) barplot_pre_phyla_filtering<-plot_bar(physeq_pre_phyla_filtering, x=colnames(sample_data(physeq_pre_phyla_filtering))[1], fill=kingdom_str) + geom_bar(stat="identity", position="stack") + labs(title=paste("Sample Depth Bar Chart",kingdom_str,sep=":"),subtitle="Sample Vs Abundance (Pre Phyla Filtering)",caption="source: Input Biom") + @@ -181,7 +187,7 @@ create_HTML<-function(htmlfile){ htmlfile_handle <- file(htmlfile) html_output = c('<html><body>', - '<table align="center>', + '<table align="center">', '<tr>', '<td valign="middle" style="vertical-align:middle;">', '<a href="pdffile.pdf"><img src="barplot_before_filtering.png"/></a>', @@ -203,7 +209,7 @@ '</td>', '</tr>', '</table>', - '</html></body>'); + '</body></html>'); writeLines(html_output, htmlfile_handle); close(htmlfile_handle); } @@ -212,7 +218,7 @@ temp.df<-data.frame(otu_table(physeq_obj)) temp.df.counts<-as.data.frame(colSums(temp.df)) colnames(temp.df.counts)<-"Abundance" - print(temp.df.counts) + #print(temp.df.counts) return(temp.df.counts) } @@ -317,12 +323,19 @@ physeq_filter_post_transform = prune_taxa((tax_table(physeq_filter_post_transform)[, kingdom_str] %in% topphyla_post_transform), physeq_filter_post_transform) + ### remove samples with zero value + otu_table(physeq_filter_post_transform)<-otu_table(physeq_filter_post_transform)[,colSums(otu_table(physeq_filter_post_transform)) > 0] create_PDF(pdffile,before_filtering_dataframe_sampleCounts,after_filtering_dataframe_sampleCounts,physeq_filter_pre_transform,physeq_filter_post_transform,kingdom_str,htmlfile,pngfile_before_filtering,pngfile_after_filtering,pngfile_pre_phyla_filtering,pngfile_post_phyla_filtering) ### convert phyloseq object to metagenomeSeq - preparing for BIOM output - metagenomeSeq_obj <- phyloseq_to_metagenomeSeq(physeq_filter_post_transform) - metagenomeSeq_biom <- MRexperiment2biom(metagenomeSeq_obj) + #metagenomeSeq_obj <- phyloseq_to_metagenomeSeq(physeq_filter_post_transform) + #metagenomeSeq_biom <- MRexperiment2biom(metagenomeSeq_obj) + + biom_obj=make_biom(otu_table(physeq_filter_post_transform),sample_metadata=sample_data(physeq_filter_post_transform),observation_metadata=tax_table(physeq_filter_post_transform),matrix_element_type="float") + biom_obj_2_metagenomeSeq_obj<-biom2MRexperiment(biom_obj) + metagenomeSeq_biom <- MRexperiment2biom(biom_obj_2_metagenomeSeq_obj) + ## write biom file write_biom(metagenomeSeq_biom, biom_file=options$outbiom)