Mercurial > repos > blankenberglab > gpca

comparison r-gpca.R @ 1:d93629e79f5e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/blankenberglab/galaxy-tools-blankenberg/tree/master/tools/r-gpca commit 3b791d44853d8928046914cd93112bdcf89d4965
author blankenberglab
date Tue, 27 Nov 2018 11:02:46 -0500
parents ed24bacaa4d0
children 4865e7f4e9c5
comparison
equal deleted inserted replaced
0:ed24bacaa4d0 1:d93629e79f5e
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 library(optparse) 2 library(optparse)
3 library(gPCA) 3 library(gPCA)
4 options(bitmapType='cairo')
4 5
5 6
6 # parse options 7 # parse options
7 option_list = list( 8 option_list = list(
8 make_option( 9 make_option(
103 } 104 }
104 if ( ! file.exists(opt$data)){ 105 if ( ! file.exists(opt$data)){
105 stop((paste('File for data', opt$data, 'does not exist'))) 106 stop((paste('File for data', opt$data, 'does not exist')))
106 } 107 }
107 108
108 batch<-as.numeric(factor(as.matrix(read.delim(opt$batch, header=TRUE, row.names=1, sep="\t")))) 109 dir.create(opt$path);
110 data<-t(as.matrix(read.delim(opt$data, row.names=1, header=TRUE, sep="\t")))
111 batch<-as.matrix(read.delim(opt$batch, header=TRUE, row.names=1, sep="\t"))
109 112
110 data<-t(as.matrix(read.delim(opt$data, row.names=1, header=TRUE, sep="\t"))) 113 write('<html><body><a name="top"></a>',file = opt$html_outfile, append=TRUE)
114 write('#Batch\tDelta\tP-value\tCumulative Variance',file = opt$numbers_outfile, append=TRUE)
115 write('<div align="left><ul id =""><p><b>Table of Contents</b></p></ul></div>', file=opt$html_outfile, append=TRUE)
111 116
112 out<-gPCA.batchdetect(x=data, batch=batch, center=opt$center, scaleY=opt$scaleY, nperm=opt$nperm, filt=opt$filt, seed=opt$seed) 117 for (row in 1:nrow(batch)) {
118 write(paste0('<li><a href="#',row.names(batch)[row],'">',row.names(batch)[row],'</a></li>'), file=opt$html_outfile, append=TRUE)
119 }
120 for (row in 1:nrow(batch)) {
121 batch1<-as.numeric(factor(batch[row,]))
122
123 out<-gPCA.batchdetect(x=data, batch=batch1, center=opt$center, scaleY=opt$scaleY, nperm=opt$nperm, filt=opt$filt, seed=opt$seed)
113 124
114 out$varPCg1<-((out$varPCg1-out$varPCu1)/out$varPCg1)*100 125 out$varPCg1<-((out$varPCg1-out$varPCu1)/out$varPCg1)*100
115 126
116 dir.create(opt$path); 127 write(paste(row.names(batch)[row],out$delta,out$p.val,out$varPCg1,sep="\t"),file = opt$numbers_outfile, append=TRUE)
117
118 write.table(data.frame(out$delta, out$p.val, out$varPCg1),file=opt$numbers_outfile, sep="\t", row.name=FALSE, quote=FALSE)
119 128
120 # General Distribution 129 # General Distribution
121 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8) 130 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8)
122 png(paste(opt$path,'/gDist.png', sep=""), width=1020, height=800, units='px') 131 png(paste(opt$path,'/gDist_',row,'.png', sep=""), width=1020, height=800, units='px')
123 gDist(out) 132 gDist(out)
124 invisible(dev.off()) 133 invisible(dev.off())
125 134
126 135
127 # Guided/Unguided 1v2 136 # Guided/Unguided 1v2
128 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8) 137 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8)
129 png(paste(opt$path,'/guided_1v2.png', sep=""), width=1020, height=800, units='px') 138 png(paste(opt$path,'/guided_1v2_',row,'.png', sep=""), width=1020, height=800, units='px')
130 PCplot(out,ug="unguided",type="1v2") 139 PCplot(out,ug="unguided",type="1v2")
131 invisible(dev.off()) 140 invisible(dev.off())
132 141
133 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8) 142 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8)
134 png(paste(opt$path,'/unguided_1v2.png', sep=""), width=1020, height=800, units='px') 143 png(paste(opt$path,'/unguided_1v2_',row,'.png', sep=""), width=1020, height=800, units='px')
135 PCplot(out,ug="unguided",type="1v2") 144 PCplot(out,ug="unguided",type="1v2")
136 invisible(dev.off()) 145 invisible(dev.off())
137 146
138 147
139 # Guided/Unguided comp,3 148 # Guided/Unguided comp,3
140 par(mai=c(0.65,0.65,0.1,0.1),cex=0.8) 149 par(mai=c(0.65,0.65,0.1,0.1),cex=0.8)
141 png(paste(opt$path,'/guided_npcs.png', sep=""), width=1020, height=800, units='px') 150 png(paste(opt$path,'/guided_npcs_',row,'.png', sep=""), width=1020, height=800, units='px')
142 PCplot(out,ug="guided",type="comp",npcs=opt$npcs) 151 PCplot(out,ug="guided",type="comp",npcs=opt$npcs)
143 invisible(dev.off()) 152 invisible(dev.off())
144 153
145 par(mai=c(0.65,0.65,0.1,0.1),cex=0.8) 154 par(mai=c(0.65,0.65,0.1,0.1),cex=0.8)
146 png(paste(opt$path,'/unguided_npcs.png', sep=""), width=1020, height=800, units='px') 155 png(paste(opt$path,'/unguided_npcs_',row,'.png', sep=""), width=1020, height=800, units='px')
147 PCplot(out,ug="unguided",type="comp",npcs=opt$npcs) 156 PCplot(out,ug="unguided",type="comp",npcs=opt$npcs)
148 invisible(dev.off()) 157 invisible(dev.off())
149 158
150 159
151 # Guided/Unguided CumlativeVarPlot 160 # Guided/Unguided CumlativeVarPlot
152 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8) 161 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8)
153 png(paste(opt$path,'/guided_var.png', sep=""), width=1020, height=800, units='px') 162 png(paste(opt$path,'/guided_var_',row,'.png', sep=""), width=1020, height=800, units='px')
154 CumulativeVarPlot(out,ug="guided",col="blue") 163 CumulativeVarPlot(out,ug="guided",col="blue")
155 invisible(dev.off()) 164 invisible(dev.off())
156 165
157 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8) 166 par(mai=c(0.8,0.8,0.1,0.1),cex=0.8)
158 png(paste(opt$path,'/unguided_var.png', sep=""), width=1020, height=800, units='px') 167 png(paste(opt$path,'/unguided_var_',row,'.png', sep=""), width=1020, height=800, units='px')
159 CumulativeVarPlot(out,ug="unguided",col="blue") 168 CumulativeVarPlot(out,ug="unguided",col="blue")
160 invisible(dev.off()) 169 invisible(dev.off())
161 170
162 171
163 write(paste('<html> 172 write(paste0('
164 <table id="r-gpca-wrap" align="center" border="1"> 173 <header id=',row.names(batch)[row],'><center><font size="11">',row.names(batch)[row],'</font><center></header>
174 <br>
175 <table align="center" border="1">
165 <tr> 176 <tr>
166 <th>Delta</th><th>P-value</th><th>varPCg1</th> 177 <th>Delta</th><th>P-value</th><th>varPCg1</th>
167 </tr> 178 </tr>
168 <tr> 179 <tr>
169 <td id=delta>',out$delta,'</td><td id=p.val>',out$p.val,'</td><td id=varPCg1>',out$varPCg1,'</td> 180 <td id="#">',out$delta,'</td><td id=p.val>',out$p.val,'</td><td id=varPCg1>',out$varPCg1,'</td>
170 </tr> 181 </tr>
171 </table> 182 </table>
183 <center><img src="gDist_',row,'.png"/></center><br>
184 <center><title>Guided 1v2</title><br><img src="guided_1v2_',row,'.png"/></center><br>
185 <center><title>Unguided 1v2</title><br><img src="unguided_1v2_',row,'.png"/></center><br>
172 186
173 <center><img src="gDist.png"/></center><br> 187 <center><title>Guided Compare to ',opt$npcs,'</title><br><img src="guided_npcs_',row,'.png"/></center><br>
174 <center><title>Guided 1v2</title><br><img src="guided_1v2.png"/></center><br> 188 <center><title>Unguided Compare to ',opt$npcs,'</title><br><img src="unguided_npcs_',row,'.png"/></center><br>
175 <center><title>Unguided 1v2</title><br><img src="unguided_1v2.png"/></center><br>
176 189
177 <center><title>Guided Compare to ',opt$npcs,'</title><br><img src="guided_npcs.png"/></center><br> 190 <center><title>Guided Cumulative Variance</title><br><img src="guided_var_',row,'.png"/></center><br>
178 <center><title>Unguided Compare to ',opt$npcs,'</title><br><img src="unguided_npcs.png"/></center><br> 191 <center><title>Unguided Cumulative Variance</title><br><img src="unguided_var_',row,'.png"/></center><br>
179 192
180 <center><title>Guided Cumulative Variance</title><br><img src="guided_var.png"/></center><br> 193 '), file = opt$html_outfile, append=TRUE)
181 <center><title>Unguided Cumulative Variance</title><br><img src="unguided_var.png"/></center><br> 194 }
182 195 write(paste('<p><center><a href="#top">Back to Top</a></center></p></body></html>'),file = opt$html_outfile, append=TRUE)
183 </html>'
184 ),file = opt$html_outfile)
185 196
186 sessionInfo() 197 sessionInfo()