Mercurial > repos > qfabrepo > metadegalaxy_phyloseq_net

comparison phyloseq_net.r @ 0:af6d9ad14a0f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/QFAB-Bioinformatics/metaDEGalaxy/tree/master/phyloseq_net commit 8bd68662b72404f6291e9628327dcb109b5fa55e"
author qfabrepo
date Mon, 14 Sep 2020 08:13:43 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:af6d9ad14a0f
1 library('getopt')
2 library('data.table')
3 suppressPackageStartupMessages(library('phyloseq'))
4 suppressPackageStartupMessages(library('DESeq2'))
5
6 Sys.setenv("DISPLAY"=":1")
7
8 options(warn= -1)
9 option_specification = matrix(c(
10 'infile','i',2,'character',
11 'metafile','m',2,'character',
12 'biom','b',2,'character',
13 'obsfile','o',2,'character',
14 'norm','n',2,'logical',
15 'xcolumn','x',2,'numeric',
16 'lcolumn','l',2,'numeric',
17 'outdir','d',2,'character',
18 'htmlfile','h',2,'character'
19 ),byrow=TRUE,ncol=4);
20
21
22 options <- getopt(option_specification)
23 options(bitmapType="cairo")
24
25 matrix.format<-function(x) {
26 m<-as.matrix(x[,-1])
27 rownames(m)<-x[,1]
28 m
29 }
30
31
32 gm_mean = function(x, na.rm=TRUE){ exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))}
33
34
35 tax_col_norm <- c("Kingdom","Phylum","Class","Order","Family","Genus","Species")
36 tax_col_extra <- c("None","Kingdom","Phylum","Class","Order","Family","Genus","Species")
37
38 tax_col_norm_otu <- c("OTUID","Kingdom","Phylum","Class","Order","Family","Genus","Species")
39 tax_col_extra_otu <- c("OTUID","None","Kingdom","Phylum","Class","Order","Family","Genus","Species")
40
41 if (!is.null(options$outdir)) {
42 # Create the directory
43 dir.create(options$outdir,FALSE)
44 }
45
46 is.biom<-options$biom
47
48
49 pdffile <- gsub("[ ]+", "", paste(options$outdir,"/pdffile.pdf"))
50 pngfile_net <- gsub("[ ]+", "", paste(options$outdir,"/net.png"))
51 htmlfile <- gsub("[ ]+", "", paste(options$htmlfile))
52
53 if(is.biom=="set_biom"){
54
55 galaxy_biom <- import_biom(options$infile)
56 galaxy_map <- import_qiime_sample_data(options$metafile)
57
58 number.of.tax.rank<-length(colnames(tax_table(galaxy_biom)))
59
60 if(number.of.tax.rank == 7){
61 colnames(tax_table(galaxy_biom)) <- tax_col_norm
62 } else {
63 colnames(tax_table(galaxy_biom)) <- tax_col_extra
64 }
65
66 physeq_galaxy <- merge_phyloseq(galaxy_biom,galaxy_map)
67
68 } else {
69
70 count.table<-read.table(options$infile,header=T,sep="\t",comment.char="",stringsAsFactors = F)
71 meta.table<-read.table(options$metafile,header=T,sep="\t",comment.char="",stringsAsFactors = F)
72 tax.table<-read.table(options$obsfile,header=T,sep="\t",comment.char="",stringsAsFactors = F)
73
74
75 colnames(count.table)<-gsub("^X","",colnames(count.table))
76 colnames(meta.table)<-gsub("^X.","",colnames(meta.table))
77 colnames(tax.table)<-gsub("^X.","",colnames(tax.table))
78
79 count.table.formatted<-matrix.format(count.table)
80 OTU<-otu_table(count.table.formatted,taxa_are_rows = TRUE)
81
82 tax.table.new<-as.data.frame(cbind(tax.table[,1],t(as.data.table(strsplit(tax.table[,2],";")))))
83
84
85 if(length(colnames(tax.table.new)) != length(tax_col_extra_otu)){
86 colnames(tax.table.new)<-tax_col_norm_otu
87 }else{
88 colnames(tax.table.new)<-tax_col_extra_otu
89 }
90
91 tax.table.formatted<-matrix.format(tax.table.new)
92
93 TAX<-tax_table(tax.table.formatted)
94
95 physeq_galaxy <- phyloseq(OTU, TAX)
96
97
98 galaxy_map<-meta.table
99
100 rownames(galaxy_map)<-meta.table[,1]
101
102 sampledata<-sample_data(as.data.frame(galaxy_map,row.names=sample_names(galaxy_map),stringsAsFactos=F))
103
104 sample_data(physeq_galaxy)<-sampledata
105
106 }
107
108
109 x.selectedColumn<-colnames(galaxy_map)[options$xcolumn]
110 l.selectedColumn<-colnames(galaxy_map)[options$lcolumn]
111
112 ### normalisation
113 if(is.null(options$norm) || options$norm =="false"){
114 suppressMessages(raw.count.deseq2.obj<-phyloseq_to_deseq2(physeq_galaxy,as.formula(paste('~',x.selectedColumn,sep=""))))
115 geoMeans = apply(counts(raw.count.deseq2.obj), 1, gm_mean)
116
117 deseq.obj = estimateSizeFactors(raw.count.deseq2.obj, geoMeans = geoMeans)
118 deseq.obj.norm<-otu_table(as.matrix(counts(deseq.obj,normalized=T)),taxa_are_rows=TRUE)
119
120 otu_table(physeq_galaxy)<-deseq.obj.norm
121 }
122
123
124 # Produce PDF file
125
126 pdf(pdffile);
127 plot_net(physeq_galaxy,point_label=x.selectedColumn,color=l.selectedColumn)
128 garbage<-dev.off();
129
130 #Cairo(pngfile_net, type="png", bg="white",pointsize=12,dpi=100,units="in",width=6,height=6)
131 png(pngfile_net,units="in",width=6,height=6,pointsize=12,res=100,bg="white")
132 plot_net(physeq_galaxy,point_label=x.selectedColumn,color=l.selectedColumn)
133 garbage<-dev.off()
134
135 # Produce the HTML file
136 htmlfile_handle <- file(htmlfile)
137 html_output = c('<html><body>',
138 '<table align="center">',
139 '<tr>',
140 '<td valign="middle" style="vertical-align:middle;">',
141 '<a href="pdffile.pdf"><img src="net.png"/></a>',
142 '</td>',
143 '</tr>',
144 '</table>',
145 '</html></body>');
146 writeLines(html_output, htmlfile_handle);
147 close(htmlfile_handle);