changeset 1:e376a618eb9f draft

Updated phyloseq_nmds.R to allow for non standard BIOM files.

--- a/phyloseq_nmds.R	Thu Jun 14 09:08:23 2018 -0400
+++ b/phyloseq_nmds.R	Sat Jun 16 05:03:43 2018 -0400
@@ -1,6 +1,8 @@
 library('getopt')
 library('ape')
 suppressPackageStartupMessages(library('phyloseq'))
+library(biomformat)
+library(plyr)
 Sys.setenv("DISPLAY"=":1")
 library("ggplot2")
 
@@ -46,15 +48,43 @@
 keep<-options$keep
 
 
+
+### This function accepts different two different type of BIOM file format
+readBIOM<-function(inBiom){
+	tryCatch({
+		phyloseq_obj<-import_biom(inBiom)
+		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)
+		}
+	)
+}
+
+
+
 if(!is.null(options$biom)){
     
-    physeq<-import_biom(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")
     }
-    
     ### select column name from sample table for nmds plot
     category_type<-colnames(sample_data(physeq))[options$subset]
 
@@ -95,7 +125,7 @@
 
     ### Remove OTUs that do not appear more than 5 times in more than half the samples
     ### filtering OTUs based on cutoff value (e.g., 5)
-    physeq_temp =genefilter_sample(physeq, filterfun_sample(function(x) x > cutoff_value), A=0.5*nsamples(physeq))
+    physeq_temp =genefilter_sample(physeq, filterfun_sample(function(x) x > cutoff_value), A=0.1*nsamples(physeq))
 
     ### phyloseq object after filtered
     physeq_filter = prune_taxa(physeq_temp, physeq)