diff phyloseq_ordinate_plot.R @ 0:ae9cd53b7760 draft

Initial upload

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/phyloseq_ordinate_plot.R	Wed Sep 05 02:07:22 2018 -0400
@@ -0,0 +1,253 @@
+library('getopt')
+library('ape')
+suppressPackageStartupMessages(library('phyloseq'))
+library(biomformat)
+library(plyr)
+Sys.setenv("DISPLAY"=":1")
+library("ggplot2")
+suppressPackageStartupMessages(library("doParallel"))
+ncores = ceiling(detectCores() * 0.8)
+registerDoParallel(cores=ncores)
+
+options(warn=-1)
+
+theme_set(theme_bw())
+
+#http://saml.rilspace.com/creating-a-galaxy-tool-for-r-scripts-that-output-images-and-pdfs
+#http://joey711.github.io/phyloseq-demo/phyloseq-demo.html
+option_specification = matrix(c(
+   'otu_table','o',2,'character',
+   'tax_table','t',2,'character',
+  'meta_table','m',2,'character',
+        'biom','b',2,'character',
+      'subset','s',2,'character',
+      'method','n',2,'character',
+    'distance','d',2,'character',
+     'kingdom','k',2,'character',
+    'plottype','e',2,'numeric',
+    'category','g',2,'numeric',
+      'outdir','r',2,'character',
+    'htmlfile','h',2,'character'
+),byrow=TRUE,ncol=4);
+
+
+options <- getopt(option_specification);
+options(bitmapType="cairo")
+
+if (!is.null(options$outdir)) {
+  # Create the directory
+  dir.create(options$outdir,FALSE)
+}
+
+
+
+method<-options$method
+### select a kingdom for phyloseq plot (e.g., "phylum")
+#kingdom_str<-colnames(tax_table)[options$kingdom]
+kingdom_str<-options$kingdom
+distance<-options$distance
+plottype<-options$plottype
+
+### prepare the directory and file name
+pdffile <- gsub("[ ]+", "", paste(options$outdir,"/pdffile.pdf"))
+pngfile_nmds <- gsub("[ ]+", "", paste(options$outdir,"/nmds.png"))
+pngfile_nmds_facet <- gsub("[ ]+", "", paste(options$outdir,"/nmds_facet.png"))
+htmlfile <- gsub("[ ]+", "", paste(options$htmlfile))
+
+
+### This function accepts different two different type of BIOM file format
+readBIOM<-function(inBiom){
+	tryCatch({
+		phyloseq_obj<-import_biom(inBiom,parallel=TRUE)
+		return(phyloseq_obj)
+		},
+		error=function(e){
+		biom_obj<-read_biom(inBiom)
+
+		otu_matrix = as(biom_data(biom_obj), "matrix")
+		OTU_TABLE = otu_table(otu_matrix, taxa_are_rows=TRUE)
+		
+		taxonomy_matrix = as.matrix(observation_metadata(biom_obj), rownames.force=TRUE)	
+		TAXONOMY_TABLE = tax_table(taxonomy_matrix)	
+		
+		metadata.temp<-sample_metadata(biom_obj)
+		METADATA_TABLE<-plyr::ldply(metadata.temp, rbind)
+        rownames(METADATA_TABLE)<-as.character(METADATA_TABLE$.id)
+
+		phyloseq_obj = phyloseq(OTU_TABLE, TAXONOMY_TABLE,sample_data(METADATA_TABLE))
+		return(phyloseq_obj)
+		}
+	)
+}
+
+
+create_OTU_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet){
+    pdf(pdf_file);
+    p1<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa")
+    print(p1)
+    p2<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3)
+    print(p2)
+    garbage<-dev.off();
+
+    #png('nmds.png')
+    bitmap(pngfile_nmds,"png16m")
+    p3<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa")
+    print(p3)
+    garbage<-dev.off()
+
+    #png('nmds_facet.png')
+    bitmap(pngfile_nmds_facet,"png16m")
+    p4<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3)
+    print(p4)
+    garbage<-dev.off()
+
+    create_HTML_1(htmlfile)
+}
+
+create_SAMPLE_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,htmlfile,pngfile_nmds,category_type){
+    pdf(pdf_file);
+    p <- plot_ordination(phyloseq_obj, phyloseq_ord, type="samples", color=category_type)
+    p <- p + geom_point(aes(fill=category_type)) + geom_point(size=5) + ggtitle("samples")
+    print(p)
+    garbage<-dev.off();
+
+    #png('nmds.png')
+    bitmap(pngfile_nmds,"png16m")
+    p1 <- plot_ordination(phyloseq_obj, phyloseq_ord, type="samples", color=category_type)
+    p1 <- p1 + geom_point(aes(fill=category_type)) + geom_point(size=5) + ggtitle("samples")
+    print(p1)
+    garbage<-dev.off();
+
+    create_HTML_2(htmlfile)
+}
+
+create_BIPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){
+    pdf(pdf_file);
+    print(category_type)
+    p_biplot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT")
+    print(p_biplot)
+    garbage<-dev.off();
+
+    bitmap(pngfile_nmds,"png16m")
+    p_biplot_png <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT")
+    print(p_biplot_png)
+    garbage<-dev.off();
+
+    create_HTML_2(htmlfile)
+}
+
+create_SPLITPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){
+   pdf(pdf_file,width=10, height=6);
+   split_plot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="split", color=kingdom_str, shape=kingdom_str, label=category_type, title="SPLIT PLOT")
+   split_plot <- split_plot + theme(plot.margin = unit(c(12,18,12,18),"pt"))
+   print(split_plot)
+   garbage<-dev.off();
+
+   bitmap(pngfile_nmds,"png16m", width=6, height=4, units="in",res=200)
+   split_plot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="split", color=kingdom_str, shape=kingdom_str, label=category_type, title="SPLIT PLOT") 
+   split_plot <- split_plot + theme(plot.margin = unit(c(12,18,12,18),"pt"))
+   print(split_plot)
+   garbage<-dev.off();
+   create_HTML_2(htmlfile)
+}
+
+create_HTML_1<-function(htmlfile){
+    htmlfile_handle <- file(htmlfile)
+    html_output = c('<html><body>',
+                    '<table align="center>',
+                    '<tr>',
+                    '<td valign="middle" style="vertical-align:middle;">',
+                        '<a href="pdffile.pdf"><img src="nmds.png"/></a>',
+                    '</td>',
+                    '</tr>',
+                    '<tr>',
+                    '<td valign="middle" style="vertical-align:middle;">',
+                    '<a href="pdffile.pdf"><img src="nmds_facet.png"/></a>',
+                    '</td>',
+                    '</tr>',
+                    '</table>',
+                    '</html></body>');
+     writeLines(html_output, htmlfile_handle);
+     close(htmlfile_handle);
+}
+
+create_HTML_2<-function(htmlfile){
+    htmlfile_handle <- file(htmlfile)
+    html_output = c('<html><body>',
+                    '<table align="center>',
+                    '<tr>',
+                    '<td valign="middle" style="vertical-align:middle;">',
+                        '<a href="pdffile.pdf"><img src="nmds.png"/></a>',
+                    '</td>',
+                    '</tr>',
+                    '</table>',
+                    '</html></body>');
+    writeLines(html_output, htmlfile_handle);
+    close(htmlfile_handle);
+}
+
+if(!is.null(options$biom)){
+    
+    #physeq<-import_biom(options$biom)
+    physeq<-readBIOM(options$biom)
+    
+    if(length(rank_names(physeq)) == 8){
+        tax_table(physeq) <- tax_table(physeq)[,-1]
+        colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
+    } else {
+        colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
+    }
+    
+    ### select column name from sample table for nmds plot
+    ## which(colnames(sample_data(biom)) == "vegetation_type_id")
+    #category_type<-colnames(sample_data(physeq))[options$subset]
+    category_type <- options$subset
+
+    ### obtain the unique value in the selected column from sample table
+    category_option<-unique(sample_data(physeq))[,options$subset]
+
+}else{
+
+    ### read the data into correct data type to create phyloseq object
+    otu_table<-as.matrix(read.table(options$otu_table,header=T,sep="\t"))
+    tax_table<-as.matrix(read.table(options$tax_table,header=T,sep="\t"))
+    sample_table<-read.table(options$meta_table,header=T,sep="\t",stringsAsFactors=F)
+
+
+    ### select column name from sample table for nmds plot
+    category_type<-colnames(sample_table)[options$category]
+
+    ### obtain the unique value in the selected column from sample table
+    category_option<-unique(sample_table[,options$category])
+
+
+    ### create a sample object for phyloseq
+    sample_object<-sample_data(sample_table)
+
+    ### create otu object for phyloseq
+    OTU<-otu_table(otu_table, taxa_are_rows = TRUE)
+
+    ### create tax object for phyloseq
+    TAX<-tax_table(tax_table)
+
+    ### create a phyloseq object
+    physeq = phyloseq(OTU,TAX,sample_object)
+}
+#   sample_data(physeq_filter)$category_input <- factor(category_input)
+    category_input = get_variable(physeq, category_type) %in% category_option
+    sample_data(physeq)$category_input <- factor(category_input)
+
+#   physeq_ord<-ordinate(physeq_filter,method,distance)
+    physeq_ord<-ordinate(physeq,method,distance)
+
+if(plottype == 1){
+    #kingdom_str = colnames(tax_table)[2]
+    create_OTU_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet) 
+}else if(plottype == 2){
+    create_SAMPLE_PDF(pdffile,physeq,physeq_ord,htmlfile,pngfile_nmds,category_type)
+}else if(plottype == 3){
+    create_BIPLOT_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type)
+}else{
+    create_SPLITPLOT_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type)
+}
+