Mercurial > repos > simon-gladman > phyloseq_ordination_plot

comparison phyloseq_ordinate_plot.R @ 1:52f009b255a1 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Updated tool
author simon-gladman
date Thu, 22 Nov 2018 07:07:27 -0500
parents ae9cd53b7760
children
comparison
equal deleted inserted replaced
0:ae9cd53b7760 1:52f009b255a1
24 'method','n',2,'character', 24 'method','n',2,'character',
25 'distance','d',2,'character', 25 'distance','d',2,'character',
26 'kingdom','k',2,'character', 26 'kingdom','k',2,'character',
27 'plottype','e',2,'numeric', 27 'plottype','e',2,'numeric',
28 'category','g',2,'numeric', 28 'category','g',2,'numeric',
29 'log','l',2,'character',
29 'outdir','r',2,'character', 30 'outdir','r',2,'character',
30 'htmlfile','h',2,'character' 31 'htmlfile','h',2,'character'
31 ),byrow=TRUE,ncol=4); 32 ),byrow=TRUE,ncol=4);
32 33
33 34
51 ### prepare the directory and file name 52 ### prepare the directory and file name
52 pdffile <- gsub("[ ]+", "", paste(options$outdir,"/pdffile.pdf")) 53 pdffile <- gsub("[ ]+", "", paste(options$outdir,"/pdffile.pdf"))
53 pngfile_nmds <- gsub("[ ]+", "", paste(options$outdir,"/nmds.png")) 54 pngfile_nmds <- gsub("[ ]+", "", paste(options$outdir,"/nmds.png"))
54 pngfile_nmds_facet <- gsub("[ ]+", "", paste(options$outdir,"/nmds_facet.png")) 55 pngfile_nmds_facet <- gsub("[ ]+", "", paste(options$outdir,"/nmds_facet.png"))
55 htmlfile <- gsub("[ ]+", "", paste(options$htmlfile)) 56 htmlfile <- gsub("[ ]+", "", paste(options$htmlfile))
56 57 output_summary <- gsub("[ ]+", "", paste(options$log))
57 58
58 ### This function accepts different two different type of BIOM file format 59 ### This function accepts different two different type of BIOM file format
59 readBIOM<-function(inBiom){ 60 readBIOM<-function(inBiom){
60 tryCatch({ 61 tryCatch({
61 phyloseq_obj<-import_biom(inBiom,parallel=TRUE) 62 phyloseq_obj<-import_biom(inBiom,parallel=TRUE)
81 } 82 }
82 83
83 84
84 create_OTU_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet){ 85 create_OTU_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet){
85 pdf(pdf_file); 86 pdf(pdf_file);
86 p1<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") 87 p1<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color=kingdom_str,title="taxa") + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
87 print(p1) 88 print(p1)
88 p2<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3) 89 p2<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color=kingdom_str,title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3) + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
89 print(p2) 90 print(p2)
90 garbage<-dev.off(); 91 garbage<-dev.off();
91 92
92 #png('nmds.png') 93 #png('nmds.png')
93 bitmap(pngfile_nmds,"png16m") 94 bitmap(pngfile_nmds,"png16m",res=120)
94 p3<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") 95 p3<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color=kingdom_str,title="taxa") + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
95 print(p3) 96 print(p3)
96 garbage<-dev.off() 97 garbage<-dev.off()
97 98
98 #png('nmds_facet.png') 99 #png('nmds_facet.png')
99 bitmap(pngfile_nmds_facet,"png16m") 100 bitmap(pngfile_nmds_facet,"png16m",res=120)
100 p4<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3) 101 p4<-plot_ordination(phyloseq_obj,phyloseq_ord,type="taxa",color=kingdom_str,title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3) + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
101 print(p4) 102 print(p4)
102 garbage<-dev.off() 103 garbage<-dev.off()
103 104
104 create_HTML_1(htmlfile) 105 create_HTML_1(htmlfile)
105 } 106 }
106 107
107 create_SAMPLE_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,htmlfile,pngfile_nmds,category_type){ 108 create_SAMPLE_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,htmlfile,pngfile_nmds,category_type){
108 pdf(pdf_file); 109 pdf(pdf_file);
109 p <- plot_ordination(phyloseq_obj, phyloseq_ord, type="samples", color=category_type) 110 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 p <- p + geom_point(aes(fill=category_type)) + geom_point(size=5) + ggtitle(paste("Samples - Stress value",formatC(phyloseq_ord$stress,digits=4,format="f"),sep=":")) + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
111 print(p) 112 print(p)
112 garbage<-dev.off(); 113 garbage<-dev.off();
113 114
114 #png('nmds.png') 115 #png('nmds.png')
115 bitmap(pngfile_nmds,"png16m") 116 bitmap(pngfile_nmds,"png16m",res=120)
116 p1 <- plot_ordination(phyloseq_obj, phyloseq_ord, type="samples", color=category_type) 117 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 p1 <- p1 + geom_point(aes(fill=category_type)) + geom_point(size=5) + ggtitle(paste("Samples - Stress value",formatC(phyloseq_ord$stress,digits=4,format="f"),sep=":")) + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
118 print(p1) 119 print(p1)
119 garbage<-dev.off(); 120 garbage<-dev.off();
120 121
121 create_HTML_2(htmlfile) 122 create_HTML_2(htmlfile)
122 } 123 }
123 124
124 create_BIPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){ 125 create_BIPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){
125 pdf(pdf_file); 126 pdf(pdf_file);
126 print(category_type) 127 print(category_type)
127 p_biplot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT") 128 p_biplot <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT") + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
128 print(p_biplot) 129 print(p_biplot)
129 garbage<-dev.off(); 130 garbage<-dev.off();
130 131
131 bitmap(pngfile_nmds,"png16m") 132 bitmap(pngfile_nmds,"png16m",res=120)
132 p_biplot_png <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT") 133 p_biplot_png <- plot_ordination(phyloseq_obj, phyloseq_ord, type="biplot", color=category_type, shape=kingdom_str,title="BIPLOT") + theme(legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
133 print(p_biplot_png) 134 print(p_biplot_png)
134 garbage<-dev.off(); 135 garbage<-dev.off();
135 136
136 create_HTML_2(htmlfile) 137 create_HTML_2(htmlfile)
137 } 138 }
138 139
139 create_SPLITPLOT_PDF<-function(pdf_file,phyloseq_obj,phyloseq_ord,kingdom_str,htmlfile,pngfile_nmds,category_type){ 140 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 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 <- 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 split_plot <- split_plot + theme(plot.margin = unit(c(12,18,12,18),"pt"),legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
143 print(split_plot) 144 print(split_plot)
144 garbage<-dev.off(); 145 garbage<-dev.off();
145 146
146 bitmap(pngfile_nmds,"png16m", width=6, height=4, units="in",res=200) 147 bitmap(pngfile_nmds,"png16m", res=120)
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 <- 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 split_plot <- split_plot + theme(plot.margin = unit(c(12,18,12,18),"pt"),legend.position="bottom",legend.box="vertical",legend.direction="horizontal")
149 print(split_plot) 150 print(split_plot)
150 garbage<-dev.off(); 151 garbage<-dev.off();
151 create_HTML_2(htmlfile) 152 create_HTML_2(htmlfile)
152 } 153 }
153 154
154 create_HTML_1<-function(htmlfile){ 155 create_HTML_1<-function(htmlfile){
155 htmlfile_handle <- file(htmlfile) 156 htmlfile_handle <- file(htmlfile)
156 html_output = c('<html><body>', 157 html_output = c('<html><body>',
157 '<table align="center>', 158 '<table align="center">',
158 '<tr>', 159 '<tr>',
159 '<td valign="middle" style="vertical-align:middle;">', 160 '<td valign="middle" style="vertical-align:middle;">',
160 '<a href="pdffile.pdf"><img src="nmds.png"/></a>', 161 '<a href="pdffile.pdf"><img src="nmds.png" width="800" height="800"/></a>',
161 '</td>', 162 '</td>',
162 '</tr>', 163 '</tr>',
163 '<tr>', 164 '<tr>',
164 '<td valign="middle" style="vertical-align:middle;">', 165 '<td valign="middle" style="vertical-align:middle;">',
165 '<a href="pdffile.pdf"><img src="nmds_facet.png"/></a>', 166 '<a href="pdffile.pdf"><img src="nmds_facet.png" width="800" height="800"/></a>',
166 '</td>', 167 '</td>',
167 '</tr>', 168 '</tr>',
168 '</table>', 169 '</table>',
169 '</html></body>'); 170 '</html></body>');
170 writeLines(html_output, htmlfile_handle); 171 writeLines(html_output, htmlfile_handle);
172 } 173 }
173 174
174 create_HTML_2<-function(htmlfile){ 175 create_HTML_2<-function(htmlfile){
175 htmlfile_handle <- file(htmlfile) 176 htmlfile_handle <- file(htmlfile)
176 html_output = c('<html><body>', 177 html_output = c('<html><body>',
177 '<table align="center>', 178 '<table align="center">',
178 '<tr>', 179 '<tr>',
179 '<td valign="middle" style="vertical-align:middle;">', 180 '<td valign="middle" style="vertical-align:middle;">',
180 '<a href="pdffile.pdf"><img src="nmds.png"/></a>', 181 '<a href="pdffile.pdf"><img src="nmds.png" width="800" height="800"/></a>',
181 '</td>', 182 '</td>',
182 '</tr>', 183 '</tr>',
183 '</table>', 184 '</table>',
184 '</html></body>'); 185 '</html></body>');
185 writeLines(html_output, htmlfile_handle); 186 writeLines(html_output, htmlfile_handle);
231 TAX<-tax_table(tax_table) 232 TAX<-tax_table(tax_table)
232 233
233 ### create a phyloseq object 234 ### create a phyloseq object
234 physeq = phyloseq(OTU,TAX,sample_object) 235 physeq = phyloseq(OTU,TAX,sample_object)
235 } 236 }
236 # sample_data(physeq_filter)$category_input <- factor(category_input) 237
237 category_input = get_variable(physeq, category_type) %in% category_option 238 category_input = get_variable(physeq, category_type) %in% category_option
238 sample_data(physeq)$category_input <- factor(category_input) 239 sample_data(physeq)$category_input <- factor(category_input)
239 240
240 # physeq_ord<-ordinate(physeq_filter,method,distance) 241 # compute distance matrix
241 physeq_ord<-ordinate(physeq,method,distance) 242 physeq_ord<-ordinate(physeq,method,distance)
243
244 # get column sum
245 sum_table<-data.frame(column_sum=as.matrix(colSums(otu_table(physeq))))
246
247 rowname_table<-data.frame(sample=rownames(sum_table))
248
249 output_table<-as.data.frame(cbind(rowname_table,sum_table))
250
251 output_table<-output_table[order(output_table$column_sum),]
252
253 # Reformat distance matrix
254 distance_matrix<-as.data.frame(physeq_ord$points)
255 distance_matrix<-cbind(sample=rownames(distance_matrix),distance_matrix)
256
257 sink(output_summary)
258 cat('--------------------------------------')
259 cat('\n')
260 cat('Stress value')
261 cat('\n')
262 cat(formatC(physeq_ord$stress,digits=4,format="f"))
263 cat('\n')
264 cat('--------------------------------------')
265 cat('\n')
266 cat('Sample - Column Sum')
267 cat('\n')
268 cat('--------------------------------------')
269 cat('\n')
270 write.table(output_table,row.names=F,quote=F)
271 cat('\n')
272 cat('--------------------------------------')
273 cat('\n')
274 cat('Distance Matrix')
275 cat('\n')
276 cat('--------------------------------------')
277 cat('\n')
278 write.table(distance_matrix,row.names=F,quote=F)
279 cat('\n')
280 cat('--------------------------------------')
281 sink()
242 282
243 if(plottype == 1){ 283 if(plottype == 1){
244 #kingdom_str = colnames(tax_table)[2] 284 #kingdom_str = colnames(tax_table)[2]
245 create_OTU_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet) 285 create_OTU_PDF(pdffile,physeq,physeq_ord,kingdom_str,htmlfile,pngfile_nmds,pngfile_nmds_facet)
246 }else if(plottype == 2){ 286 }else if(plottype == 2){