Mercurial > repos > simon-gladman > phyloseq_ordination_plot

comparison phyloseq_ordinate_plot.R @ 0:ae9cd53b7760 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial upload
author simon-gladman
date Wed, 05 Sep 2018 02:07:22 -0400
parents
children 52f009b255a1
comparison
equal deleted inserted replaced
-1:000000000000 0:ae9cd53b7760
1 library('getopt')
2 library('ape')
3 suppressPackageStartupMessages(library('phyloseq'))
4 library(biomformat)
5 library(plyr)
6 Sys.setenv("DISPLAY"=":1")
7 library("ggplot2")
8 suppressPackageStartupMessages(library("doParallel"))
9 ncores = ceiling(detectCores() * 0.8)
10 registerDoParallel(cores=ncores)
11
12 options(warn=-1)
13
14 theme_set(theme_bw())
15
16 #http://saml.rilspace.com/creating-a-galaxy-tool-for-r-scripts-that-output-images-and-pdfs
17 #http://joey711.github.io/phyloseq-demo/phyloseq-demo.html
18 option_specification = matrix(c(
19 'otu_table','o',2,'character',
20 'tax_table','t',2,'character',
21 'meta_table','m',2,'character',
22 'biom','b',2,'character',
23 'subset','s',2,'character',
24 'method','n',2,'character',
25 'distance','d',2,'character',
26 'kingdom','k',2,'character',
27 'plottype','e',2,'numeric',
28 'category','g',2,'numeric',
29 'outdir','r',2,'character',
30 'htmlfile','h',2,'character'
31 ),byrow=TRUE,ncol=4);
32
33
34 options <- getopt(option_specification);
35 options(bitmapType="cairo")
36
37 if (!is.null(options$outdir)) {
38 # Create the directory
39 dir.create(options$outdir,FALSE)
40 }
41
42
43
44 method<-options$method
45 ### select a kingdom for phyloseq plot (e.g., "phylum")
46 #kingdom_str<-colnames(tax_table)[options$kingdom]
47 kingdom_str<-options$kingdom
48 distance<-options$distance
49 plottype<-options$plottype
50
51 ### prepare the directory and file name
52 pdffile <- gsub("[ ]+", "", paste(options$outdir,"/pdffile.pdf"))
53 pngfile_nmds <- gsub("[ ]+", "", paste(options$outdir,"/nmds.png"))
54 pngfile_nmds_facet <- gsub("[ ]+", "", paste(options$outdir,"/nmds_facet.png"))
55 htmlfile <- gsub("[ ]+", "", paste(options$htmlfile))
56
57
58 ### This function accepts different two different type of BIOM file format
59 readBIOM<-function(inBiom){
60 tryCatch({
61 phyloseq_obj<-import_biom(inBiom,parallel=TRUE)
62 return(phyloseq_obj)
63 },
64 error=function(e){
65 biom_obj<-read_biom(inBiom)
66
67 otu_matrix = as(biom_data(biom_obj), "matrix")
68 OTU_TABLE = otu_table(otu_matrix, taxa_are_rows=TRUE)
69
70 taxonomy_matrix = as.matrix(observation_metadata(biom_obj), rownames.force=TRUE)
71 TAXONOMY_TABLE = tax_table(taxonomy_matrix)
72
73 metadata.temp<-sample_metadata(biom_obj)
74 METADATA_TABLE<-plyr::ldply(metadata.temp, rbind)
75 rownames(METADATA_TABLE)<-as.character(METADATA_TABLE$.id)
76
77 phyloseq_obj = phyloseq(OTU_TABLE, TAXONOMY_TABLE,sample_data(METADATA_TABLE))
78 return(phyloseq_obj)
79 }
80 )
81 }
82
83
84 create_OTU_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet){
85 pdf(pdf_file);
86 p1<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa")
87 print(p1)
88 p2<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3)
89 print(p2)
90 garbage<-dev.off();
91
92 #png('nmds.png')
93 bitmap(pngfile_nmds,"png16m")
94 p3<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa")
95 print(p3)
96 garbage<-dev.off()
97
98 #png('nmds_facet.png')
99 bitmap(pngfile_nmds_facet,"png16m")
100 p4<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3)
101 print(p4)
102 garbage<-dev.off()
103
104 create_HTML_1(htmlfile)
105 }
106
107 create_SAMPLE_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,htmlfile,pngfile_nmds,category_type){
108 pdf(pdf_file);
109 p <- plot_ordination(phyloseq_obj, phyloseq_ord, type="samples", color=category_type)
110 p <- p + geom_point(aes(fill=category_type)) + geom_point(size=5) + ggtitle("samples")
111 print(p)
112 garbage<-dev.off();
113
114 #png('nmds.png')
115 bitmap(pngfile_nmds,"png16m")
116 p1 <- plot_ordination(phyloseq_obj, phyloseq_ord, type="samples", color=category_type)
117 p1 <- p1 + geom_point(aes(fill=category_type)) + geom_point(size=5) + ggtitle("samples")
118 print(p1)
119 garbage<-dev.off();
120
121 create_HTML_2(htmlfile)
122 }
123
124 create_BIPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){
125 pdf(pdf_file);
126 print(category_type)
127 p_biplot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT")
128 print(p_biplot)
129 garbage<-dev.off();
130
131 bitmap(pngfile_nmds,"png16m")
132 p_biplot_png <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT")
133 print(p_biplot_png)
134 garbage<-dev.off();
135
136 create_HTML_2(htmlfile)
137 }
138
139 create_SPLITPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){
140 pdf(pdf_file,width=10, height=6);
141 split_plot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="split", color=kingdom_str, shape=kingdom_str, label=category_type, title="SPLIT PLOT")
142 split_plot <- split_plot + theme(plot.margin = unit(c(12,18,12,18),"pt"))
143 print(split_plot)
144 garbage<-dev.off();
145
146 bitmap(pngfile_nmds,"png16m", width=6, height=4, units="in",res=200)
147 split_plot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="split", color=kingdom_str, shape=kingdom_str, label=category_type, title="SPLIT PLOT")
148 split_plot <- split_plot + theme(plot.margin = unit(c(12,18,12,18),"pt"))
149 print(split_plot)
150 garbage<-dev.off();
151 create_HTML_2(htmlfile)
152 }
153
154 create_HTML_1<-function(htmlfile){
155 htmlfile_handle <- file(htmlfile)
156 html_output = c('<html><body>',
157 '<table align="center>',
158 '<tr>',
159 '<td valign="middle" style="vertical-align:middle;">',
160 '<a href="pdffile.pdf"><img src="nmds.png"/></a>',
161 '</td>',
162 '</tr>',
163 '<tr>',
164 '<td valign="middle" style="vertical-align:middle;">',
165 '<a href="pdffile.pdf"><img src="nmds_facet.png"/></a>',
166 '</td>',
167 '</tr>',
168 '</table>',
169 '</html></body>');
170 writeLines(html_output, htmlfile_handle);
171 close(htmlfile_handle);
172 }
173
174 create_HTML_2<-function(htmlfile){
175 htmlfile_handle <- file(htmlfile)
176 html_output = c('<html><body>',
177 '<table align="center>',
178 '<tr>',
179 '<td valign="middle" style="vertical-align:middle;">',
180 '<a href="pdffile.pdf"><img src="nmds.png"/></a>',
181 '</td>',
182 '</tr>',
183 '</table>',
184 '</html></body>');
185 writeLines(html_output, htmlfile_handle);
186 close(htmlfile_handle);
187 }
188
189 if(!is.null(options$biom)){
190
191 #physeq<-import_biom(options$biom)
192 physeq<-readBIOM(options$biom)
193
194 if(length(rank_names(physeq)) == 8){
195 tax_table(physeq) <- tax_table(physeq)[,-1]
196 colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
197 } else {
198 colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
199 }
200
201 ### select column name from sample table for nmds plot
202 ## which(colnames(sample_data(biom)) == "vegetation_type_id")
203 #category_type<-colnames(sample_data(physeq))[options$subset]
204 category_type <- options$subset
205
206 ### obtain the unique value in the selected column from sample table
207 category_option<-unique(sample_data(physeq))[,options$subset]
208
209 }else{
210
211 ### read the data into correct data type to create phyloseq object
212 otu_table<-as.matrix(read.table(options$otu_table,header=T,sep="\t"))
213 tax_table<-as.matrix(read.table(options$tax_table,header=T,sep="\t"))
214 sample_table<-read.table(options$meta_table,header=T,sep="\t",stringsAsFactors=F)
215
216
217 ### select column name from sample table for nmds plot
218 category_type<-colnames(sample_table)[options$category]
219
220 ### obtain the unique value in the selected column from sample table
221 category_option<-unique(sample_table[,options$category])
222
223
224 ### create a sample object for phyloseq
225 sample_object<-sample_data(sample_table)
226
227 ### create otu object for phyloseq
228 OTU<-otu_table(otu_table, taxa_are_rows = TRUE)
229
230 ### create tax object for phyloseq
231 TAX<-tax_table(tax_table)
232
233 ### create a phyloseq object
234 physeq = phyloseq(OTU,TAX,sample_object)
235 }
236 # sample_data(physeq_filter)$category_input <- factor(category_input)
237 category_input = get_variable(physeq, category_type) %in% category_option
238 sample_data(physeq)$category_input <- factor(category_input)
239
240 # physeq_ord<-ordinate(physeq_filter,method,distance)
241 physeq_ord<-ordinate(physeq,method,distance)
242
243 if(plottype == 1){
244 #kingdom_str = colnames(tax_table)[2]
245 create_OTU_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet)
246 }else if(plottype == 2){
247 create_SAMPLE_PDF(pdffile,physeq,physeq_ord,htmlfile,pngfile_nmds,category_type)
248 }else if(plottype == 3){
249 create_BIPLOT_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type)
250 }else{
251 create_SPLITPLOT_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type)
252 }
253