|
2
|
1 #!/usr/bin/Rscript
|
|
|
2
|
|
|
3 ## --------------------
|
|
|
4 ## prints all arguments as msg
|
|
|
5 ## --------------------
|
|
|
6 catMsg <- function( msg=c() ){
|
|
|
7 cat( MAIN_NAME, paste( msg, collapse="" ), "\n", sep='')
|
|
|
8 }
|
|
|
9 ## --------------------
|
|
28
|
10 ## return the location of this script
|
|
2
|
11 ## --------------------
|
|
|
12 getScriptPath <- function(){
|
|
|
13 cmd.args <- commandArgs()
|
|
|
14 m <- regexpr("(?<=^--file=).+", cmd.args, perl=TRUE)
|
|
|
15 script.dir <- dirname(regmatches(cmd.args, m))
|
|
42
|
16 if( length(script.dir) == 0 ) stop("[ERR] Can't determine script dir: please call the script with Rscript\n")
|
|
|
17 if( length(script.dir) > 1 ) stop("[ERR] Can't determine script dir: more than one '--file' argument detected\n")
|
|
2
|
18 return(script.dir)
|
|
|
19 }
|
|
|
20 ## --------------------
|
|
28
|
21 ## Some html creation functions
|
|
2
|
22 ## --------------------
|
|
|
23 htmlTableRow <- function( string_array=c() ){
|
|
|
24 td_cells <- ''
|
|
|
25 for ( i in string_array ){
|
|
|
26 td_cells <- paste( td_cells, '<td>', i, '</td>', sep='' )
|
|
|
27 }
|
|
|
28 return( paste( "<tr>", td_cells, "</tr>") )
|
|
|
29 }
|
|
|
30 htmlLink <- function( path, desc="LINK" ){
|
|
|
31 return( paste( '<a href="', path, '">', desc, "</a>", sep='') )
|
|
|
32 }
|
|
|
33 ## --------------------
|
|
28
|
34 ## constructs a list with input bam file info
|
|
2
|
35 ## --------------------
|
|
|
36 makeBamFileList <- function( paths, names ){
|
|
|
37 tmp <- list()
|
|
|
38 l1 <- length(paths)
|
|
|
39 l2 <- length(names)
|
|
42
|
40 if ( l1 != l2 ) stop( "Unequal amount of bam-paths (", l1, ") and -names (", l2, ") in makeBamFileList!!!\n" )
|
|
2
|
41 for ( i in 1:length(paths) ){
|
|
|
42 path <- paths[i]
|
|
|
43 name <- names[i]
|
|
|
44 file <- basename(path)
|
|
|
45
|
|
|
46 tmp[[ file ]] <- name
|
|
|
47 tmp[[ 'all_paths' ]] <- c( tmp[[ 'all_paths' ]], path )
|
|
|
48 tmp[[ 'all_files' ]] <- c( tmp[[ 'all_files' ]], file )
|
|
|
49 tmp[[ 'all_names' ]] <- c( tmp[[ 'all_names' ]], name )
|
|
|
50 }
|
|
|
51 return( tmp )
|
|
|
52 }
|
|
|
53
|
|
|
54 ## --------------------
|
|
28
|
55 ## copied code for extracting the regions by segment call status
|
|
2
|
56 ## --------------------
|
|
|
57 fuse.regions_test <- function(cgh, onlyCalled=T) {
|
|
|
58 if (ncol(cgh) > 1) stop('Please specify which sample...')
|
|
|
59
|
|
|
60 x <- data.frame(cgh@featureData@data[,1:3], calls(cgh), copynumber(cgh), segmented(cgh), check.names=FALSE, stringsAsFactors=FALSE)
|
|
|
61 colnames( x ) <- c( "chr", "start", "end", "call", "log2", "segmentval" )
|
|
|
62
|
|
|
63 fused.data <- data.frame()
|
|
|
64 curr.bin <- 1
|
|
|
65 for (chr in unique(x$chr)) {
|
|
|
66
|
|
|
67 chr.data <- x[x$chr == chr,]
|
|
|
68 prev.bin <- curr.bin
|
|
|
69 prev.call <- chr.data[1, 'call']
|
|
|
70 prev.log2 <- chr.data[1, 'log2']
|
|
|
71 prev.segm <- chr.data[1, 'segmentval']
|
|
|
72 start <- chr.data[1, 'start']
|
|
|
73
|
|
|
74 if ( nrow(chr.data) > 1) {
|
|
|
75 for ( i in 2:nrow(chr.data) ) {
|
|
|
76
|
|
|
77 curr.bin <- curr.bin + 1
|
|
|
78 curr.call <- chr.data[ i, 'call']
|
|
|
79 curr.segm <- chr.data[ i, 'segmentval']
|
|
|
80
|
|
|
81 if ( curr.segm != prev.segm ) {
|
|
|
82
|
|
|
83 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
84 if ( prev.call != 0 ){
|
|
42
|
85 catMsg( c(" ...found called/segmented region (", chr, ':', start, ' call=', prev.call, ' segment=', prev.segm, ")" ) )
|
|
2
|
86 }
|
|
|
87 prev.call <- curr.call
|
|
|
88 prev.segm <- curr.segm
|
|
|
89 prev.bin <- curr.bin
|
|
|
90 start <- chr.data[ i, 'start']
|
|
|
91 }
|
|
|
92 }
|
|
|
93 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
94 }else {
|
|
|
95 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
96 }
|
|
|
97 }
|
|
|
98 if ( onlyCalled == T ){
|
|
|
99 fused.data <- fused.data[ fused.data$call != 0, ]
|
|
|
100 }
|
|
|
101 fused.data
|
|
|
102 }
|
|
|
103
|
|
|
104 ## DESC: takes the output of fuse.regions and outputs a txt file per sample
|
|
|
105 outputRegionsFromList <- function ( regionsList, outputBasename, outputDir="./" ){
|
|
|
106 if ( missing(regionsList) ) stop( 'Please provide regionsList...' )
|
|
|
107 if ( missing(outputBasename) ) stop( 'Please provide outputBasename...' )
|
|
|
108 if ( !is.list(regionsList) ) stop( 'Input not a list...?' )
|
|
|
109 if ( length(regionsList) < 1 ) stop( 'List seems empty...?' )
|
|
42
|
110 if ( file.exists( outputDir ) ) catMsg( c(" Using dir ", outputDir, " for output") )
|
|
2
|
111 else dir.create( outputDir )
|
|
|
112 outFiles <- list()
|
|
|
113
|
|
|
114 ## have to set R output options otherwise scientific method is used at some point
|
|
|
115 options( "scipen"=100 )
|
|
|
116
|
|
|
117 sampleCount <- length( regionsList )
|
|
|
118 sampleNames <- names( regionsList )
|
|
|
119 bedgraphColumns <- c( 'chr', 'start', 'end', 'segmentval' )
|
|
30
|
120
|
|
42
|
121 catMsg( c( " There are ", sampleCount, " samples found in input list") )
|
|
2
|
122
|
|
|
123 for ( sample in sampleNames ){
|
|
42
|
124 catMsg( c(" Working on sample ", sample ) )
|
|
2
|
125 regionCount <- nrow( regionsList[[sample]] )
|
|
|
126
|
|
|
127 outSampleBase <- paste( outputBasename, '_', sample, '_QDNAseqRegions', sep='')
|
|
|
128 outBedFile <- paste( outSampleBase, '.bed', sep="" )
|
|
|
129 outBedPath <- paste( outputDir, '/', outBedFile, sep="" )
|
|
|
130 outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
|
|
|
131 outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )
|
|
|
132
|
|
|
133 ## ---------- BED ----------
|
|
|
134 txt <- "#"
|
|
|
135 sink( outBedPath )
|
|
|
136 cat( txt )
|
|
|
137 sink()
|
|
|
138 write.table( regionsList[[sample]], outBedPath, quote=F, sep="\t", row.names=F, append=T)
|
|
|
139
|
|
|
140 ## ---------- BEDGRAPH ----------
|
|
|
141 txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample,"_segmReg description=segmented_regions_from_QDNAseq\n", sep="")
|
|
|
142 sink( outBedgraphPath )
|
|
|
143 cat( txt )
|
|
|
144 sink()
|
|
|
145 write.table( regionsList[[sample]][,bedgraphColumns], outBedgraphPath, quote=F, sep="\t", row.names=F, append=T, col.names=F)
|
|
|
146 outFiles[[sample]] <- c( outBedFile, outBedgraphFile )
|
|
|
147 }
|
|
|
148 outFiles
|
|
|
149 }
|
|
|
150
|
|
|
151 ## ==================================================
|
|
|
152 ## Start of analysis
|
|
|
153 ## ==================================================
|
|
42
|
154 MAIN_NAME <- '[INFO] '
|
|
2
|
155 TOOL_PATH <- getScriptPath()
|
|
|
156 CSS_FILE <- paste( TOOL_PATH, '/QDNAseq.css', sep="" )
|
|
|
157 DECIMALS <- 3
|
|
|
158 WEB_LINK <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
|
|
54
|
159 PAR_SET <- list( pch=22 )
|
|
2
|
160
|
|
|
161 catMsg( "Starting QDNAseq wrapper" )
|
|
|
162 catMsg( "Loading R libraries" )
|
|
|
163 suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
|
|
|
164 suppressWarnings( suppressMessages( library( CGHcall, quietly = TRUE ) ) )
|
|
25
|
165
|
|
2
|
166 ## only one param: the tmp config file
|
|
|
167 cmdLineArgs <- commandArgs(TRUE)
|
|
|
168 config <- cmdLineArgs[1]
|
|
|
169
|
|
|
170 ## sourcing the config file will load all input params
|
|
42
|
171 ## many variables are created in sourced "config"
|
|
2
|
172 source( config )
|
|
|
173
|
|
42
|
174 PLOT_RES <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 6.3 # desparate try to make png text scale well, damn you R...!
|
|
2
|
175 if ( doOutputCallsRds == TRUE ){ doCall <- TRUE }
|
|
|
176
|
|
42
|
177 systemUser <- system("whoami",T)
|
|
|
178 qdnaseqVersion <- packageDescription( "QDNAseq" )$Version
|
|
|
179 catMsg( c("Analysis run with user: ", systemUser ) )
|
|
|
180 catMsg( c("QDNAseq version loaded: ", qdnaseqVersion) )
|
|
45
|
181 catMsg( c( sessionInfo()$R.version$version.string ) )
|
|
42
|
182
|
|
2
|
183 ## get the comma separated list of chromosomes to exclude
|
|
|
184 excludeChrs <- unlist( strsplit( excludeChrsString, ",") )
|
|
|
185
|
|
|
186 ## ------------------------
|
|
|
187 ## DEBUG
|
|
|
188 #catMsg( "PARAM" )
|
|
|
189 #catMsg( galaxy_path )
|
|
|
190 #catMsg( repos_path )
|
|
|
191 #catMsg( instal_path )
|
|
|
192 ## /DEBUG
|
|
|
193 ## ------------------------
|
|
|
194
|
|
|
195 ## setup bam filelist for easy retrieval later
|
|
|
196 fileList <- makeBamFileList( bamsPaths, bamsNames )
|
|
|
197 bamCount <- length( fileList[[ 'all_paths' ]] )
|
|
|
198
|
|
42
|
199 ## help msg for running script separate still needs work!
|
|
28
|
200 #if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){
|
|
42
|
201 # catMsg( c( "Usage: ", paste(params_help,sep=",") ) )
|
|
28
|
202 # quit( save = 'no', status=0, runLast=F )
|
|
|
203 #}
|
|
2
|
204
|
|
|
205 if ( !file.exists( outputPath) ){
|
|
|
206 dir.create( outputPath )
|
|
|
207 }
|
|
|
208
|
|
|
209 ## because we circumvent params that galaxy can save, we want to
|
|
|
210 ## copy source config file to output dir to include it in output zip
|
|
|
211 file.copy( config, paste(outputPath, 'qdnaseq_config_file.R', sep='/') )
|
|
|
212
|
|
|
213 ## ------------------------
|
|
|
214 ## construct output file-names and -paths
|
|
|
215 ## ------------------------
|
|
|
216 htmlOutputName <- 'index.html'
|
|
|
217 gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
|
|
|
218 robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
|
|
|
219 robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
|
|
|
220 robjCalledName <- paste( binSize, 'kbp_QDNAseqCalledSegments.rds', sep='')
|
|
|
221 regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
|
|
|
222 igvCalledName <- paste( binSize, 'kbp_QDNAseq-calls.igv', sep='')
|
|
|
223 igvCopyNrName <- paste( binSize, 'kbp_QDNAseq-normalized.igv', sep='')
|
|
42
|
224 noiseImgName <- paste( 'QDNAseqNoisePlot.png', sep='')
|
|
|
225 freqImgName <- paste( 'QDNAseqFrequencyPlot.png', sep='')
|
|
2
|
226
|
|
42
|
227
|
|
|
228 gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep='')
|
|
|
229 htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep='')
|
|
|
230 robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep='')
|
|
|
231 robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep='')
|
|
|
232 robjCalledPath <- paste( outputPath, '/', robjCalledName, sep='')
|
|
|
233 robjRegionPath <- paste( outputPath, '/', regiOutputName, sep='')
|
|
|
234 igvCalledPath <- paste( outputPath, '/', igvCalledName, sep='')
|
|
|
235 igvCopyNrPath <- paste( outputPath, '/', igvCopyNrName, sep='')
|
|
|
236 noiseImgPath <- paste( outputPath, '/', noiseImgName, sep='' )
|
|
|
237 freqImgPath <- paste( outputPath, '/', freqImgName, sep='' )
|
|
2
|
238
|
|
|
239
|
|
|
240 ## ------------------------
|
|
|
241 ## performing QDNAseq analysis steps
|
|
|
242 ## ------------------------
|
|
|
243 if ( debug ){
|
|
|
244 ## in case of debug just use inbuilt LGG data for speedup
|
|
|
245 data(LGG150)
|
|
|
246 readCounts <- LGG150
|
|
|
247 }else{
|
|
|
248 if ( nchar(binAnnotations) == 0 ){
|
|
|
249 binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
|
|
|
250 }else{
|
|
|
251 ## if user provided file, check if correct class
|
|
|
252 if ( class(binAnnotations)[1] != 'AnnotatedDataFrame' ){
|
|
|
253 stop( "Provided binAnnotations file is not of class 'AnnotatedDataFrame'\n" )
|
|
|
254 }
|
|
|
255 binAnnotations <- readRDS( binAnnotations )
|
|
|
256 }
|
|
|
257 ## provide bamnames because in galaxy everyting is called "dataset_###"
|
|
|
258 readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=bamsNames )
|
|
|
259 }
|
|
|
260
|
|
|
261 readCountsFiltered <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
|
|
|
262 readCountsFiltered <- estimateCorrection( readCountsFiltered )
|
|
|
263 copyNumbers <- correctBins( readCountsFiltered )
|
|
|
264 copyNumbersNormalized <- normalizeBins( copyNumbers )
|
|
|
265 copyNumbersSmooth <- smoothOutlierBins( copyNumbersNormalized )
|
|
28
|
266 sampleNames <- readCountsFiltered@phenoData@data$name
|
|
2
|
267
|
|
|
268 ## save objects to output dir
|
|
42
|
269 saveRDS( readCountsFiltered, robjReadCoPath );
|
|
29
|
270 saveRDS( copyNumbersSmooth, robjCopyNrPath );
|
|
42
|
271 exportBins( copyNumbersSmooth, file=igvCopyNrPath, format="igv" )
|
|
2
|
272
|
|
|
273 ## also save objects for galaxy history output if requested
|
|
|
274 if ( doOutputReadcountsRds ){
|
|
|
275 saveRDS( readCountsFiltered, readCountsDatasetFile );
|
|
|
276 }
|
|
|
277 if ( doOutputCopynumbersRds ){
|
|
|
278 saveRDS( copyNumbersSmooth, copyNumbersDatasetFile );
|
|
|
279 }
|
|
|
280
|
|
|
281 ## proceed with calling if requested
|
|
|
282 if ( doCall ){
|
|
|
283 copyNumbersSegmented <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
|
|
|
284 copyNumbersSegmented <- normalizeSegmentedBins( copyNumbersSegmented )
|
|
|
285 copyNumbersCalled <- callBins( copyNumbersSegmented )
|
|
|
286 cgh <- makeCgh( copyNumbersCalled )
|
|
|
287 saveRDS( copyNumbersCalled, robjCalledPath );
|
|
|
288 if ( doOutputCallsRds ){
|
|
|
289 saveRDS( copyNumbersCalled, calledSegmentsDatasetFile );
|
|
|
290 }
|
|
31
|
291 exportBins( copyNumbersCalled, file=igvCalledPath, format="igv")
|
|
42
|
292
|
|
2
|
293 }
|
|
|
294
|
|
28
|
295
|
|
2
|
296
|
|
|
297 ## ------------------------
|
|
|
298 ## create output files
|
|
|
299 ## ------------------------
|
|
|
300 plotted_images <- list() # to keep track of images for later linking
|
|
|
301 regions <- list() # will contain the (called) segments
|
|
|
302
|
|
53
|
303 png( noiseImgPath, width=PLOT_HEIGHT, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
304 par( PAR_SET )
|
|
42
|
305 noisePlot( readCountsFiltered, col="darkgreen" )
|
|
2
|
306 dev.off()
|
|
|
307
|
|
42
|
308 if ( doCall ){
|
|
|
309 png( freqImgPath, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
53
|
310 par( PAR_SET )
|
|
42
|
311 frequencyPlot( copyNumbersCalled )
|
|
|
312 dev.off()
|
|
|
313 }
|
|
|
314
|
|
|
315
|
|
2
|
316 for (i in 1:length(sampleNames) ){
|
|
|
317 #for (sample in sampleNames(copyNumbersSmooth) ){
|
|
|
318 sample <- sampleNames[i]
|
|
|
319 usedReads <- readCountsFiltered@phenoData@data$used.reads[i]
|
|
|
320 catMsg( c("Creating plots for sample: ", sample ) )
|
|
|
321
|
|
|
322 type <- 'CopyNumbers'
|
|
|
323 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
|
|
|
324 img_file_path <- paste( outputPath, '/', img_file, sep='' )
|
|
42
|
325 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
53
|
326 par( PAR_SET )
|
|
42
|
327 plot( copyNumbersSmooth[ ,sample ] );
|
|
|
328 dev.off()
|
|
2
|
329 plotted_images[[ sample ]][[ type ]] <- img_file
|
|
|
330
|
|
|
331 if ( doCall ){
|
|
|
332 type <- 'Called'
|
|
|
333 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
|
|
|
334 img_file_path <- paste( outputPath, '/', img_file, sep='' )
|
|
42
|
335 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
53
|
336 par( PAR_SET )
|
|
2
|
337 plot( copyNumbersCalled[ ,sample ] );
|
|
|
338 dev.off()
|
|
|
339 plotted_images[[ sample ]][[ type ]] <- img_file
|
|
|
340
|
|
42
|
341 catMsg( c(" Fusing regions of sample: ", sample) )
|
|
2
|
342 regions[[ sample ]] <- fuse.regions_test( cgh[, sample] )
|
|
|
343 region_count <- nrow( data.frame( regions[[ sample ]] ) )
|
|
42
|
344 catMsg( c( ' sample "', sample, '" has ', region_count, " regions" ) )
|
|
2
|
345 plotted_images[[ sample ]][[ 'region_count' ]] <- region_count
|
|
|
346 }
|
|
|
347
|
|
30
|
348 ## add USED read counts
|
|
2
|
349 plotted_images[[ sample ]][[ 'usedReads' ]] <- usedReads
|
|
|
350 }
|
|
|
351
|
|
|
352 if ( doCall ){
|
|
|
353 saveRDS( regions, robjRegionPath )
|
|
|
354 printedFiles <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath )
|
|
|
355 }
|
|
|
356
|
|
|
357 ## ------------------------
|
|
|
358 ## prepare output
|
|
|
359 ## ------------------------
|
|
42
|
360 catMsg( "...zipping output")
|
|
2
|
361 zip_cmd <- paste( "zip -j", gzipOutputPath, paste(outputPath,'/*',sep='') ) ## -j is for removing dirs from the tree
|
|
|
362 system( zip_cmd )
|
|
|
363
|
|
|
364 ## ------------------------
|
|
|
365 ## get filesizes for report
|
|
|
366 ## ------------------------
|
|
|
367 zippedSize <- paste( round( file.info( gzipOutputPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
|
368 readCoSize <- paste( round( file.info( robjReadCoPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
|
369 copyNrSize <- paste( round( file.info( robjCopyNrPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
|
370 calledSize <- paste( round( file.info( robjCalledPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
|
371 regionSize <- paste( round( file.info( robjRegionPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
31
|
372 igvCopyNrSize <- paste( round( file.info( igvCopyNrPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
34
|
373 igvCalledSize <- paste( round( file.info( igvCalledPath )[["size"]] / 1000000, digits=2 ), 'MB' )
|
|
2
|
374
|
|
|
375 ## ------------------------
|
|
|
376 ## creating html output to be linked to from the middle galaxy pane
|
|
|
377 ## ------------------------
|
|
|
378
|
|
|
379 #sink( file = outputHtml, type = "output" )
|
|
|
380 sink( file = htmlOutputPath, type = "output" )
|
|
|
381 cat( "<html>\n")
|
|
|
382 cat( "<head>\n")
|
|
|
383 #cat( '<link rel="stylesheet" type="text/css" href="test.css" media="screen" />', "\n" )
|
|
|
384 cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
|
|
|
385 cat( "\t<style>\n")
|
|
|
386 ## have to include CSS into html file, because css referencing outside own dir doesn't seem to work...
|
|
|
387 cat( paste( "\t\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
|
|
|
388 cat( paste( "\t\t", readLines( CSS_FILE, n = -1)), sep="\n" )
|
|
|
389 #cat( "\t\th1 {color:red;}", "\n")
|
|
|
390
|
|
|
391 cat( "\t</style>\n" )
|
|
|
392
|
|
|
393 cat( "\n</head>\n")
|
|
|
394 cat( "\n<body>\n")
|
|
|
395
|
|
|
396 cat( "<h1>QDNAseq Report</h1>", "\n")
|
|
|
397
|
|
|
398 cat( '<h3 class="qdnaseq">About this analysis</h3>', "\n")
|
|
42
|
399 cat( '<p>This page provides access to all results. To have a local copy of this report just download the <a href="', gzipOutputName, '" class="button-success button-small pure-button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep='')
|
|
2
|
400
|
|
|
401 ## ------------------------
|
|
|
402 ## table with general info
|
|
|
403 ## ------------------------
|
|
|
404 cat( '<h3 class="qdnaseq">Settings</h3><p>', "\n")
|
|
|
405 cat( '<table class="pure-table pure-table-striped"><thead><tr><th>setting</th><th>value</th></tr></thead><tbody>' )
|
|
|
406 cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
|
|
|
407 cat( htmlTableRow( c( "AnalysisDate", 'todo' ) ) )
|
|
|
408 cat( htmlTableRow( c( "BinSize (kb)", binSize ) ) )
|
|
|
409
|
|
|
410 for ( galaxyName in fileList[[ 'all_files' ]] ){
|
|
|
411 sampleName <- fileList[[ galaxyName ]]
|
|
|
412 cat( htmlTableRow( c( "InputBam", paste( galaxyName, ' (', sampleName, ')', sep='' ) ) ) )
|
|
|
413 }
|
|
|
414 cat( "</tbody></table></p>", "\n")
|
|
|
415
|
|
|
416 ## ------------------------
|
|
|
417 ## list with links to all output files
|
|
|
418 ## ------------------------
|
|
|
419 r_code <- '<p>'
|
|
|
420 r_code <- paste( r_code, '<code class="code">## R code to load files</code><br />', sep="\n" )
|
|
|
421 r_code <- paste( r_code, '<code class="code">library( QDNAseq )</code><br />', sep="\n")
|
|
|
422
|
|
|
423 cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
|
|
|
424 cat( '<dl>', "\n" )
|
|
|
425 cat( '<dt>', htmlLink( path=robjReadCoName, robjReadCoName ), '</dt>', "\n" )
|
|
|
426 cat( '<dd>QDNAseq object with read counts per bin, ', readCoSize,'</dd>', "\n" )
|
|
|
427 r_code <- paste( r_code, '<code class="code">readCounts <- readRDS(', robjReadCoName, ")</code><br />", sep="")
|
|
|
428
|
|
|
429 cat( '<dt>', htmlLink( path=robjCopyNrName, robjCopyNrName ), '</dt>', "\n" )
|
|
|
430 cat( '<dd>QDNAseq object with copy numbers per bin, ', copyNrSize,'</dd>', "\n" )
|
|
|
431 r_code <- paste( r_code, '<code class="code">copyNumbersSmooth <- readRDS(', robjCopyNrName, ")</code><br />", sep="")
|
|
|
432
|
|
31
|
433 cat( '<dt>', htmlLink( path=igvCopyNrName, igvCopyNrName ), '</dt>', "\n" )
|
|
|
434 cat( '<dd>IGV formatted text file with copy numbers per bin, ', igvCopyNrSize,'</dd>', "\n" )
|
|
|
435
|
|
2
|
436 if ( doCall ){
|
|
|
437 cat( '<dt>', htmlLink( path=robjCalledName, robjCalledName ), '</dt>', "\n" )
|
|
|
438 cat( '<dd>QDNAseq object with segment and call status per bin, ', calledSize,'</dd>', "\n" )
|
|
|
439 r_code <- paste( r_code, '<code class="code">copyNumbersCalled <- readRDS(', robjCalledName, ")</code><br />", sep="")
|
|
|
440
|
|
|
441 cat( '<dt>', htmlLink( path=regiOutputName, regiOutputName ), '</dt>', "\n" )
|
|
|
442 cat( '<dd>list with segmented/called regions for each sample, ', regionSize, '</dd>', "\n" )
|
|
|
443 r_code <- paste( r_code, '<code class="code">calledRegions <- readRDS(', regiOutputName, ")</code><br />", sep="")
|
|
|
444
|
|
31
|
445 cat( '<dt>', htmlLink( path=igvCalledName, igvCalledName ), '</dt>', "\n" )
|
|
|
446 cat( '<dd>IGV formatted text file with calls per bin , ', igvCalledSize,'</dd>', "\n" )
|
|
|
447
|
|
2
|
448 }
|
|
|
449 cat( '</dl></p>', "\n" )
|
|
|
450
|
|
42
|
451 ## r_code shows example code on how to load output files
|
|
25
|
452 #cat( r_code, "</p>\n", sep="\n")
|
|
2
|
453 cat( '<p>See ', htmlLink( WEB_LINK, 'the bioconductor QDNAseq documentation' ), ' for more information on how to work with these files</p>', "\n", sep='' )
|
|
|
454
|
|
|
455 ## ------------------------
|
|
|
456 ## table with links to files
|
|
|
457 ## ------------------------
|
|
|
458 cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
|
|
|
459 plots_html <- ''
|
|
|
460
|
|
40
|
461 cat( '<table class="pure-table pure-table-striped"><thead><tr><th>Sample / File</th><th>CopyNr</th><th>Calls</th><th>Reads</th><th>Regions</th><th>Files</th></tr></thead><tbody>' )
|
|
2
|
462
|
|
|
463 for ( bam_file in sampleNames ){
|
|
|
464
|
|
|
465 #width <- 600; height <- 240
|
|
|
466 width <- PLOT_WIDTH; height <- PLOT_HEIGHT
|
|
|
467 width_t <- 100; height_t <- 40
|
|
|
468
|
|
|
469 ## add thumbnails to table with links to anchors on html page
|
|
|
470 copy_img <- plotted_images[[ bam_file ]][[ 'CopyNumbers' ]]
|
|
|
471 usedReads <- plotted_images[[ bam_file ]][[ 'usedReads' ]]
|
|
|
472 usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )
|
|
|
473
|
|
|
474 html_copy_thumb <- htmlLink( path=paste('#', copy_img, sep=''), paste('<img src="',copy_img,'" alt="', bam_file, '" width="', width_t, '" height="', height_t, '">', sep='') )
|
|
|
475 html_copy_img <- htmlLink( path=copy_img, paste('<img id="', copy_img,'" src="',copy_img,'" alt="',bam_file, '" width="', width, '" height="', height, '">', sep='') )
|
|
|
476 html_call_thumb <- 'NA'
|
|
|
477 html_call_img <- ''
|
|
|
478 html_bed <- 'NA'
|
|
|
479 html_bedGraph <- 'NA'
|
|
|
480 region_count <- 'NA'
|
|
|
481
|
|
|
482 if ( doCall ){
|
|
|
483 call_img <- plotted_images[[ bam_file ]][[ 'Called' ]]
|
|
|
484 region_count <- plotted_images[[ bam_file ]][[ 'region_count' ]]
|
|
|
485 html_call_thumb <- htmlLink( path=paste('#', call_img, sep=''), paste('<img src="', call_img, '" alt="', bam_file, '" width="', width_t,'" height="', height_t,'">', sep='') )
|
|
|
486
|
|
|
487 files <- printedFiles[[ bam_file ]]
|
|
|
488 html_bed <- htmlLink( path=files[1], 'bed' )
|
|
|
489 html_bedGraph <- htmlLink( path=files[2], 'bedGraph' )
|
|
|
490 html_call_img <- htmlLink( path=call_img, paste('<img id="', call_img,'" src="', call_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
|
|
|
491 }
|
|
|
492
|
|
|
493 ## add info to overview table, including small thumbnails
|
|
|
494 cat( htmlTableRow( c(bam_file, html_copy_thumb, html_call_thumb, usedReads, region_count, paste( html_bed, html_bedGraph, sep=", ")) ), "\n" )
|
|
|
495 ## now include (large) images in html page
|
|
|
496 plots_html <- paste( plots_html, html_copy_img, "\n", html_call_img, "\n<hr \\>\n", sep='' )
|
|
|
497 }
|
|
|
498 cat( "</tbody></table></p>", "\n")
|
|
|
499
|
|
42
|
500 ## ------------------------
|
|
|
501 ## create the noise and frequency plots
|
|
|
502 ## ------------------------
|
|
|
503 html_noise_img <- htmlLink(
|
|
|
504 path=noiseImgName,
|
|
55
|
505 paste('<img id="', noiseImgName,'" src="', noiseImgName, '" width="', height/2, '" height="', height/2, '" alt="NoisePlot">', sep='')
|
|
42
|
506 )
|
|
|
507 html_freq_img <- htmlLink(
|
|
|
508 path=freqImgName,
|
|
55
|
509 paste('<img id="', freqImgName,'" src="', freqImgName, '" width="', width/2, '" height="', height/2, ' alt="FrequenceyPlot">', sep='')
|
|
42
|
510 )
|
|
|
511 extra_plots_html <- paste(
|
|
55
|
512 html_noise_img, " \n",
|
|
42
|
513 html_freq_img, "\n", sep=''
|
|
|
514 )
|
|
2
|
515
|
|
|
516 ## ------------------------
|
|
|
517 ## section with various output shown
|
|
|
518 ## ------------------------
|
|
42
|
519 cat( '<h3 class="qdnaseq">Results: Analysis plots</h3><p>', "\n")
|
|
|
520 cat( extra_plots_html, "\n")
|
|
|
521 cat( '<h3 class="qdnaseq">Results: Sample plots</h3><p>', "\n")
|
|
2
|
522 cat( plots_html, "\n")
|
|
|
523 cat( "\n</p></body>\n")
|
|
|
524 cat( "\n</html>\n")
|
|
|
525 sink()
|
|
|
526
|
|
|
527 ## ------------------------
|
|
42
|
528 ## creating main html output for galaxy history
|
|
2
|
529 ## ------------------------
|
|
|
530 sink( file = outputHtml, type = "output" )
|
|
|
531
|
|
|
532 cat( "<head>", "\n")
|
|
|
533 cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
|
|
|
534
|
|
|
535 cat( "<style>", "\n")
|
|
42
|
536 ## include CSS directly into html file
|
|
2
|
537 cat( paste( "\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
|
|
|
538 cat( paste( "\t", readLines( CSS_FILE, n = -1)), sep="\n" )
|
|
|
539 cat( "</style>", "\n")
|
|
|
540 cat( "</head>", "\n")
|
|
|
541
|
|
|
542 cat( '<h1>QDNAseq Results (', outputName,')</h1>', "\n", sep="")
|
|
|
543 cat( '<p>Explore <a href="', htmlOutputName, '" class="button-success button-small pure-button">the results</a> directly within galaxy</p>', "\n", sep="")
|
|
|
544 cat( '<p>Or download a <a href="', gzipOutputName, '" class="button-success button-small pure-button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )
|
|
|
545
|
|
|
546 cat( '<p>The zip file contains all output files, including *.rds files allowing you to load the R copyNumber object(s) and perform further detailed analysis or create your own output for further processing. You can load the rds file with <code class="code">loadRDS(file.rds)</code></p>', "\n", sep="")
|
|
|
547
|
|
|
548 sink()
|
|
|
549
|
|
42
|
550 ## ------------------------
|
|
|
551 ## create final zip and quit with status 0 to tell galaxy all was fine
|
|
|
552 ## ------------------------
|
|
|
553 catMsg( "zipping all output")
|
|
|
554 system( paste( "zip -j ", gzipOutputPath, paste(outputPath,'/', htmlOutputName, sep='') ) )
|
|
|
555 catMsg( "done" )
|
|
25
|
556 q(status=0)
|
