|
0
|
1 library('getopt')
|
|
|
2 library('ape')
|
|
|
3 suppressPackageStartupMessages(library('phyloseq'))
|
|
|
4 Sys.setenv("DISPLAY"=":1")
|
|
|
5 library("ggplot2")
|
|
|
6
|
|
|
7 options(warn=-1)
|
|
|
8
|
|
|
9 theme_set(theme_bw())
|
|
|
10
|
|
|
11 #http://saml.rilspace.com/creating-a-galaxy-tool-for-r-scripts-that-output-images-and-pdfs
|
|
|
12 #http://joey711.github.io/phyloseq-demo/phyloseq-demo.html
|
|
|
13 option_specification = matrix(c(
|
|
|
14 'otu_table','o',2,'character',
|
|
|
15 'tax_table','t',2,'character',
|
|
|
16 'meta_table','s',2,'character',
|
|
|
17 'biom','i',2,'character',
|
|
|
18 'subset','x',2,'numeric',
|
|
|
19 'method','n',2,'character',
|
|
|
20 'distance','d',2,'character',
|
|
|
21 'kingdom','k',2,'character',
|
|
|
22 'cutoff','v',2,'numeric',
|
|
|
23 'category','c',2,'numeric',
|
|
|
24 'keep','p',2,'numeric',
|
|
|
25 'outdir','r',2,'character',
|
|
|
26 'htmlfile','h',2,'character'
|
|
|
27 ),byrow=TRUE,ncol=4);
|
|
|
28
|
|
|
29
|
|
|
30 options <- getopt(option_specification);
|
|
|
31 options(bitmapType="cairo")
|
|
|
32
|
|
|
33 if (!is.null(options$outdir)) {
|
|
|
34 # Create the directory
|
|
|
35 dir.create(options$outdir,FALSE)
|
|
|
36 }
|
|
|
37
|
|
|
38
|
|
|
39
|
|
|
40 method<-options$method
|
|
|
41 cutoff_value<-options$cutoff
|
|
|
42 ### select a kingdom for phyloseq plot (e.g., "phylum")
|
|
|
43 #kingdom_str<-colnames(tax_table)[options$kingdom]
|
|
|
44 kingdom_str<-options$kingdom
|
|
|
45 distance<-options$distance
|
|
|
46 keep<-options$keep
|
|
|
47
|
|
|
48
|
|
|
49 if(!is.null(options$biom)){
|
|
|
50
|
|
|
51 physeq<-import_biom(options$biom)
|
|
|
52
|
|
|
53 if(length(rank_names(physeq)) == 8){
|
|
|
54 tax_table(physeq) <- tax_table(physeq)[,-1]
|
|
|
55 colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
|
|
|
56 }
|
|
|
57
|
|
|
58 ### select column name from sample table for nmds plot
|
|
|
59 category_type<-colnames(sample_data(physeq))[options$subset]
|
|
|
60
|
|
|
61 ### obtain the unique value in the selected column from sample table
|
|
|
62 category_option<-unique(sample_data(physeq))[,options$subset]
|
|
|
63
|
|
|
64 }else{
|
|
|
65
|
|
|
66 ### read the data into correct data type to create phyloseq object
|
|
|
67 otu_table<-as.matrix(read.table(options$otu_table,header=T,sep="\t"))
|
|
|
68 tax_table<-as.matrix(read.table(options$tax_table,header=T,sep="\t"))
|
|
|
69 sample_table<-read.table(options$meta_table,header=T,sep="\t",stringsAsFactors=F)
|
|
|
70
|
|
|
71
|
|
|
72 ### select column name from sample table for nmds plot
|
|
|
73 category_type<-colnames(sample_table)[options$category]
|
|
|
74
|
|
|
75 ### obtain the unique value in the selected column from sample table
|
|
|
76 category_option<-unique(sample_table[,options$category])
|
|
|
77
|
|
|
78
|
|
|
79 ### create a sample object for phyloseq
|
|
|
80 sample_object<-sample_data(sample_table)
|
|
|
81
|
|
|
82 ### create otu object for phyloseq
|
|
|
83 OTU<-otu_table(otu_table, taxa_are_rows = TRUE)
|
|
|
84
|
|
|
85 ### create tax object for phyloseq
|
|
|
86 TAX<-tax_table(tax_table)
|
|
|
87
|
|
|
88 ### create a phyloseq object
|
|
|
89 physeq = phyloseq(OTU,TAX,sample_object)
|
|
|
90 }
|
|
|
91 ### select a kingdom for phyloseq plot (e.g., "phylum")
|
|
|
92 #kingdom_str<-colnames(tax_table)[options$kingdom]
|
|
|
93 #kingdom_str<-options$kingdom
|
|
|
94
|
|
|
95
|
|
|
96 ### Remove OTUs that do not appear more than 5 times in more than half the samples
|
|
|
97 ### filtering OTUs based on cutoff value (e.g., 5)
|
|
|
98 physeq_temp =genefilter_sample(physeq, filterfun_sample(function(x) x > cutoff_value), A=0.5*nsamples(physeq))
|
|
|
99
|
|
|
100 ### phyloseq object after filtered
|
|
|
101 physeq_filter = prune_taxa(physeq_temp, physeq)
|
|
|
102
|
|
|
103 ## Transform to even sampling depth
|
|
|
104 physeq_filter = transform_sample_counts(physeq_filter, function(x) 1E6 * x/sum(x))
|
|
|
105
|
|
|
106 ## Keep only the most abundant five phyla
|
|
|
107 phylum.sum = tapply(taxa_sums(physeq_filter), tax_table(physeq_filter)[, kingdom_str], sum, na.rm=TRUE)
|
|
|
108
|
|
|
109 ### number of most abundance phyla to keep
|
|
|
110 topphyla = names(sort(phylum.sum, TRUE))[1:keep]
|
|
|
111 physeq_filter = prune_taxa((tax_table(physeq_filter)[, kingdom_str] %in% topphyla), physeq_filter)
|
|
|
112
|
|
|
113 ### select category to plot NMDS
|
|
|
114 category_input = get_variable(physeq_filter, category_type) %in% category_option
|
|
|
115 sample_data(physeq_filter)$category_input <- factor(category_input)
|
|
|
116
|
|
|
117 ### prepare the directory and file name
|
|
|
118 pdffile <- gsub("[ ]+", "", paste(options$outdir,"/pdffile.pdf"))
|
|
|
119 pngfile_nmds <- gsub("[ ]+", "", paste(options$outdir,"/nmds.png"))
|
|
|
120 pngfile_nmds_facet <- gsub("[ ]+", "", paste(options$outdir,"/nmds_facet.png"))
|
|
|
121 htmlfile <- gsub("[ ]+", "", paste(options$htmlfile))
|
|
|
122
|
|
|
123
|
|
|
124 # Produce PDF file
|
|
|
125 pdf(pdffile);
|
|
|
126 physeq_ord<-ordinate(physeq_filter,method,distance)
|
|
|
127 plot_ordination(physeq,physeq_ord,type="taxa",color="Phylum",title="taxa")
|
|
|
128 plot_ordination(physeq,physeq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3)
|
|
|
129 garbage<-dev.off();
|
|
|
130
|
|
|
131
|
|
|
132 #png('nmds.png')
|
|
|
133 bitmap(pngfile_nmds,"png16m")
|
|
|
134 plot_ordination(physeq,physeq_ord,type="taxa",color="Phylum",title="taxa")
|
|
|
135 garbage<-dev.off()
|
|
|
136
|
|
|
137 #png('nmds_facet.png')
|
|
|
138 bitmap(pngfile_nmds_facet,"png16m")
|
|
|
139 plot_ordination(physeq,physeq_ord,type="taxa",color="Phylum",title="taxa") + facet_wrap(formula(paste('~',kingdom_str)),3)
|
|
|
140 garbage<-dev.off()
|
|
|
141
|
|
|
142 # Produce the HTML file
|
|
|
143 htmlfile_handle <- file(htmlfile)
|
|
|
144 html_output = c('<html><body>',
|
|
|
145 '<table align="center>',
|
|
|
146 '<tr>',
|
|
|
147 '<td valign="middle" style="vertical-align:middle;">',
|
|
|
148 '<a href="pdffile.pdf"><img src="nmds.png"/></a>',
|
|
|
149 '</td>',
|
|
|
150 '</tr>',
|
|
|
151 '<tr>',
|
|
|
152 '<td valign="middle" style="vertical-align:middle;">',
|
|
|
153 '<a href="pdffile.pdf"><img src="nmds_facet.png"/></a>',
|
|
|
154 '</td>',
|
|
|
155 '</tr>',
|
|
|
156 '</table>',
|
|
|
157 '</html></body>');
|
|
|
158 writeLines(html_output, htmlfile_handle);
|
|
|
159 close(htmlfile_handle);
|
