Mercurial > repos > qfabrepo > metadegalaxy_symmetric_plot
comparison symmetric_plot.r @ 0:0cac08094b86 draft
"planemo upload for repository https://github.com/QFAB-Bioinformatics/metaDEGalaxy/tree/master/symmetric_plot commit a68579d7bdde7420b8f04346d3b6e361588acf50"
author | qfabrepo |
---|---|
date | Mon, 14 Sep 2020 06:04:44 +0000 |
parents | |
children | 91fb94d203df |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:0cac08094b86 |
---|---|
1 library('getopt') | |
2 library('tidyr') | |
3 suppressPackageStartupMessages(library('dplyr')) | |
4 suppressPackageStartupMessages(library('phyloseq')) | |
5 suppressPackageStartupMessages(library('DESeq2')) | |
6 suppressPackageStartupMessages(library('ggplot2')) | |
7 suppressPackageStartupMessages(library('data.table')) | |
8 Sys.setenv("DISPLAY"=":1") | |
9 | |
10 | |
11 options(warn= -1) | |
12 option_specification = matrix(c( | |
13 'input.data','i',2,'character', | |
14 'meta.data','m',2,'character', | |
15 'obs.data','t',2,'character', | |
16 'record','r',2,'numeric', | |
17 'taxrank','x',2,'character', | |
18 'norm','n',2,'logical', | |
19 'n.column','c',2,'numeric', | |
20 'g.group','g',2,'character', | |
21 'outdir','o',2,'character', | |
22 'htmlfile','h',2,'character' | |
23 ),byrow=TRUE,ncol=4); | |
24 | |
25 | |
26 options <- getopt(option_specification); | |
27 options(bitmapType="cairo") | |
28 | |
29 | |
30 if (!is.null(options$outdir)) { | |
31 # Create the directory | |
32 dir.create(options$outdir,FALSE) | |
33 } | |
34 | |
35 | |
36 input.table<-read.table(options$input.data,header=T,sep="\t",stringsAsFactors = F) | |
37 metadata.table<-read.table(options$meta.data,header=T,sep="\t",stringsAsFactors = F,comment.char="") | |
38 obs.table<-read.table(options$obs.data,header=F,sep="\t",stringsAsFactors = F,comment.char="") | |
39 | |
40 colnames(obs.table)<-c("OTUID","taxonomy") | |
41 | |
42 | |
43 tax_col <- c("OTUID","Kingdom","Phylum","Class","Order","Family","Genus","Species") | |
44 tax_col_extra <- c("OTUID","None","Kingdom","Phylum","Class","Order","Family","Genus","Species") | |
45 | |
46 ### remove the leading #sign in column name if the column name begins with number | |
47 colnames(input.table)<-gsub("^X","",colnames(input.table)) | |
48 colnames(metadata.table)<-gsub("^X.","",colnames(metadata.table)) | |
49 | |
50 | |
51 column.name<-colnames(metadata.table)[options$n.column] | |
52 #in.column<-options$g.column | |
53 in.group<-options$g.group | |
54 nrecord<-options$record | |
55 ranking<-options$taxrank | |
56 | |
57 ### create data frame for group | |
58 group.df<-data.frame(group=unlist(strsplit(in.group,",")),stringsAsFactors = F) | |
59 | |
60 ### get the number of group | |
61 number.of.group<-dim(group.df)[1] | |
62 | |
63 if(number.of.group != 2){ | |
64 print(paste("Number of group for comparision is",number.of.group,sep="")) | |
65 quit("yes") | |
66 } | |
67 | |
68 | |
69 | |
70 group1<-metadata.table[which(metadata.table[,column.name] %in% group.df$group[1]),]$SampleID | |
71 group2<-metadata.table[which(metadata.table[,column.name] %in% group.df$group[2]),]$SampleID | |
72 | |
73 sample2group.map <- data.frame(sample=colnames(input.table[,c(group1,group2)]), | |
74 groups=c(rep(group.df$group[1],length(group1)), | |
75 rep(group.df$group[2],length(group2)))) | |
76 | |
77 | |
78 if(options$norm =="false"){ | |
79 ### raw count table | |
80 count.table<-input.table | |
81 rownames(count.table)<-count.table[,1] | |
82 count.table<-count.table[,-1] | |
83 | |
84 suppressMessages(raw.count.deseq.obj<-DESeqDataSetFromMatrix(countData = count.table,colData=metadata.table, as.formula(paste('~',column.name,sep="")))) | |
85 | |
86 gm_mean = function(x, na.rm=TRUE){ | |
87 exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x)) | |
88 } | |
89 | |
90 geoMeans = apply(counts(raw.count.deseq.obj), 1, gm_mean) | |
91 deseq_obj = estimateSizeFactors(raw.count.deseq.obj, geoMeans = geoMeans) | |
92 deseq_obj_norm<-counts(deseq_obj,normalized=T) | |
93 } else { | |
94 deseq_obj_norm<-input.table | |
95 rownames(deseq_obj_norm)<-deseq_obj_norm[,1] | |
96 deseq_obj_norm<-deseq_obj_norm[,-1] | |
97 } | |
98 | |
99 | |
100 #first.50.otu.id<-rownames(counts(raw.count.deseq.obj))[1:50] | |
101 first.50.otu.id<-rownames(deseq_obj_norm)[1:nrecord] | |
102 | |
103 | |
104 #filtered.data<-counts(raw.count.deseq.obj)[first.50.otu.id,c(group1,group2)] | |
105 filtered.data<-deseq_obj_norm[first.50.otu.id,c(group1,group2)] | |
106 filtered.data<-as.data.frame(cbind(OTUID=rownames(filtered.data),filtered.data),stringsAsFactors=F) | |
107 | |
108 nc <- match('OTUID',colnames(filtered.data)) | |
109 filtered.data[,-nc]<-sapply(filtered.data[,-nc],as.integer) | |
110 stopifnot(min(range(filtered.data[,-nc]))>=0) | |
111 | |
112 long <- gather(filtered.data,sample,expr,-OTUID) | |
113 suppressMessages((long <- left_join(long, sample2group.map))) | |
114 long$expr[long$groups == group.df$group[1]] <- long$expr*-1 | |
115 | |
116 sorted.OTU <- rev(sort(unique(long$OTUID))) | |
117 long$OTUID <- factor(long$OTUID, levels=sorted.OTU) | |
118 long$sample <- as.factor(long$sample) | |
119 | |
120 | |
121 tax.table.new<-as.data.frame(cbind(obs.table[,1],t(as.data.table(strsplit(obs.table[,2],";"))))) | |
122 | |
123 if(length(colnames(tax.table.new)) != length(tax_col_extra)) | |
124 { | |
125 colnames(tax.table.new)<-tax_col | |
126 } else { | |
127 colnames(tax.table.new)<-tax_col_extra | |
128 } | |
129 | |
130 long<-cbind(long,tax.table.new[match(long$OTUID,tax.table.new$OTUID),-1]) | |
131 | |
132 comparison<-paste(column.name,paste(group.df$group[1],group.df$group[2],sep="-"),sep=" ") | |
133 | |
134 | |
135 p<-ggplot(long,aes(x=reorder(OTUID,expr),y=expr, fill=groups)) + | |
136 geom_bar(stat='identity') + theme_bw() + | |
137 xlab("OTU ID") + | |
138 labs(title=comparison) + facet_grid( as.formula(paste(ranking,"~ .",sep="")),scales = "free", space = "free" ) + theme(strip.text.y = element_text(angle = 0)) + | |
139 coord_flip() | |
140 q<-ggplot_build(p) | |
141 q$layout$panel_ranges[[1]]$x.labels <- gsub("-","",q$layout$panel_ranges[[1]]$x.labels) | |
142 #Reassemble the plot using ggplot_gtable() | |
143 q<-ggplot_gtable(q) | |
144 | |
145 | |
146 | |
147 | |
148 pdffile <- gsub("[ ]+", "", paste(options$outdir,"/symmetric.pdf")) | |
149 pngfile_symmetric <- gsub("[ ]+", "", paste(options$outdir,"/symmetric.png")) | |
150 htmlfile <- gsub("[ ]+", "", paste(options$htmlfile)) | |
151 | |
152 | |
153 # Produce PDF file | |
154 pdf(pdffile); | |
155 plot(q) | |
156 garbage<-dev.off(); | |
157 | |
158 #png('richness.png') | |
159 bitmap(pngfile_symmetric,"png16m",height=10,width=10,res=100) | |
160 plot(q) | |
161 garbage<-dev.off() | |
162 | |
163 # Produce the HTML file | |
164 htmlfile_handle <- file(htmlfile) | |
165 html_output = c('<html><body>', | |
166 '<table align="center">', | |
167 '<tr>', | |
168 '<td valign="middle" style="vertical-align:middle;">', | |
169 '<a href="pdffile.pdf"><img src="symmetric.png"/></a>', | |
170 '</td>', | |
171 '</tr>', | |
172 '</table>', | |
173 '<table align="center>', | |
174 '<tr>', | |
175 '<td valign="middle" style="vertical-align:middle;">', | |
176 '<a href="pdffile.pdf"><img src="symmetric.png"/></a>', | |
177 '</td>', | |
178 '</tr>', | |
179 '</table>', | |
180 '</html></body>'); | |
181 writeLines(html_output, htmlfile_handle); | |
182 close(htmlfile_handle); |