Mercurial > repos > stef > qdnaseq

comparison QDNAseq.R @ 2:336697c6f7fa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author stef
date Fri, 13 Jun 2014 09:42:10 -0400
parents
children 8509c112abaa
comparison
equal deleted inserted replaced
1:75d96e0555d1 2:336697c6f7fa
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 ## --------------------
10 ## return the directory this script exists
11 ## --------------------
12 getScriptPath <- function(){
13 cmd.args <- commandArgs()
14 m <- regexpr("(?<=^--file=).+", cmd.args, perl=TRUE)
15 script.dir <- dirname(regmatches(cmd.args, m))
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")
18 return(script.dir)
19 }
20 ## --------------------
21 ## Some html helper functions
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 ## --------------------
34 ## constructs a list with bam file info
35 ## --------------------
36 makeBamFileList <- function( paths, names ){
37 tmp <- list()
38 l1 <- length(paths)
39 l2 <- length(names)
40 if ( l1 != l2 ) stop( "Unequal amount of bam-paths (",l1,") and -names (",l2,") in makeBamFileList!!!\n" )
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 ## --------------------
55 ## took a function from Matias/MarkVdWuel? for extracting the regions by call
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 ){
85 cat( MAIN_NAME, " ...found called/segmented region (", chr, ':', start, ' call=', prev.call, ' segment=', prev.segm, ")\n", sep="" )
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...?' )
110 if ( file.exists( outputDir ) ) cat( MAIN_NAME, " Using dir ", outputDir, " for output\n", sep="")
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 systemUser <- system("whoami",T)
120 bedgraphColumns <- c( 'chr', 'start', 'end', 'segmentval' )
121
122 cat( MAIN_NAME, " Hello ", systemUser, "!!\n", sep="")
123 cat( MAIN_NAME, " There are ", sampleCount, " samples found in input list...\n", sep="")
124
125 for ( sample in sampleNames ){
126 cat( MAIN_NAME, " Working on sample ", sample, "\n", sep="")
127 regionCount <- nrow( regionsList[[sample]] )
128
129 outSampleBase <- paste( outputBasename, '_', sample, '_QDNAseqRegions', sep='')
130 outBedFile <- paste( outSampleBase, '.bed', sep="" )
131 outBedPath <- paste( outputDir, '/', outBedFile, sep="" )
132 outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
133 outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )
134
135 ## ---------- BED ----------
136 txt <- "#"
137 sink( outBedPath )
138 cat( txt )
139 sink()
140 write.table( regionsList[[sample]], outBedPath, quote=F, sep="\t", row.names=F, append=T)
141
142 ## ---------- BEDGRAPH ----------
143 txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample,"_segmReg description=segmented_regions_from_QDNAseq\n", sep="")
144 sink( outBedgraphPath )
145 cat( txt )
146 sink()
147 write.table( regionsList[[sample]][,bedgraphColumns], outBedgraphPath, quote=F, sep="\t", row.names=F, append=T, col.names=F)
148 outFiles[[sample]] <- c( outBedFile, outBedgraphFile )
149 }
150 outFiles
151 }
152
153 printIgvFile <- function( dat, filename ){
154
155 if ( 'calls' %in% assayDataElementNames(dat) ) {
156 #output <- paste(output, '-called.igv', sep="")
157 cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
158 df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), calls(dat), check.names=FALSE, stringsAsFactors=FALSE)
159 }else {
160 #output <- paste(output, '-normalized.igv', sep="")
161 cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
162 df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), round(copynumber(dat), digits=2), check.names=FALSE, stringsAsFactors=FALSE)
163 }
164
165 df$chromosome[df$chromosome == '23'] <- 'X'
166 df$chromosome[df$chromosome == '24'] <- 'Y'
167 df$chromosome[df$chromosome == '25'] <- 'MT'
168 #return( df )
169 write.table( df, file=filename, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )
170 }
171
172
173 ## ==================================================
174 ## Start of analysis
175 ## ==================================================
176 TOOL_PATH <- getScriptPath()
177 MAIN_NAME <- '[INFO] '
178 CSS_FILE <- paste( TOOL_PATH, '/QDNAseq.css', sep="" )
179 DECIMALS <- 3
180 WEB_LINK <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
181
182 catMsg( "Starting QDNAseq wrapper" )
183 catMsg( "Loading R libraries" )
184 suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
185 suppressWarnings( suppressMessages( library( CGHcall, quietly = TRUE ) ) )
186 suppressWarnings( suppressMessages( library( matrixStats, quietly = TRUE ) ) )
187
188 ## only one param: the tmp config file
189 cmdLineArgs <- commandArgs(TRUE)
190 config <- cmdLineArgs[1]
191
192 ## sourcing the config file will load all input params
193 source( config )
194
195 ## if call output requested, set doCall such that we will segment and call
196 if ( doOutputCallsRds == TRUE ){ doCall <- TRUE }
197
198 ## get the comma separated list of chromosomes to exclude
199 excludeChrs <- unlist( strsplit( excludeChrsString, ",") )
200
201 ## ------------------------
202 ## DEBUG
203 #catMsg( "PARAM" )
204 #catMsg( galaxy_path )
205 #catMsg( repos_path )
206 #catMsg( instal_path )
207 ## /DEBUG
208 ## ------------------------
209
210 ## setup bam filelist for easy retrieval later
211 catMsg( "Setting up input bam list" )
212 cat( bamsPaths, "\n" )
213 catMsg( "Namews" )
214 cat( bamsNames, "\n" )
215 fileList <- makeBamFileList( bamsPaths, bamsNames )
216 bamCount <- length( fileList[[ 'all_paths' ]] )
217
218 ## help msg still needs work!
219 if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){
220 cat( paste( MAIN_NAME, "Usage: ", params_help, sep=''), "\n" )
221 quit( save = 'no', status=0, runLast=F )
222 }
223
224 if ( !file.exists( outputPath) ){
225 dir.create( outputPath )
226 }
227
228 ## because we circumvent params that galaxy can save, we want to
229 ## copy source config file to output dir to include it in output zip
230 file.copy( config, paste(outputPath, 'qdnaseq_config_file.R', sep='/') )
231
232 ## ------------------------
233 ## construct output file-names and -paths
234 ## ------------------------
235 htmlOutputName <- 'index.html'
236 gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
237 robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
238 robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
239 robjCalledName <- paste( binSize, 'kbp_QDNAseqCalledSegments.rds', sep='')
240 regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
241 igvCalledName <- paste( binSize, 'kbp_QDNAseq-calls.igv', sep='')
242 igvCopyNrName <- paste( binSize, 'kbp_QDNAseq-normalized.igv', sep='')
243
244 gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep="")
245 htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep="")
246 robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep="")
247 robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep="")
248 robjCalledPath <- paste( outputPath, '/', robjCalledName, sep="")
249 robjRegionPath <- paste( outputPath, '/', regiOutputName, sep="")
250 igv_calledPath <- paste( outputPath, '/', igvCalledName, sep="")
251 igv_copyNrPath <- paste( outputPath, '/', igvCopyNrName, sep="")
252
253
254 ## ------------------------
255 ## performing QDNAseq analysis steps
256 ## ------------------------
257 if ( debug ){
258 ## in case of debug just use inbuilt LGG data for speedup
259 data(LGG150)
260 readCounts <- LGG150
261 }else{
262 if ( nchar(binAnnotations) == 0 ){
263 binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
264 }else{
265 ## if user provided file, check if correct class
266 if ( class(binAnnotations)[1] != 'AnnotatedDataFrame' ){
267 stop( "Provided binAnnotations file is not of class 'AnnotatedDataFrame'\n" )
268 }
269 binAnnotations <- readRDS( binAnnotations )
270 }
271 ## provide bamnames because in galaxy everyting is called "dataset_###"
272 readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=bamsNames )
273 }
274
275 readCountsFiltered <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
276 readCountsFiltered <- estimateCorrection( readCountsFiltered )
277 copyNumbers <- correctBins( readCountsFiltered )
278 copyNumbersNormalized <- normalizeBins( copyNumbers )
279 copyNumbersSmooth <- smoothOutlierBins( copyNumbersNormalized )
280
281 ## save objects to output dir
282 saveRDS( readCounts, robjReadCoPath );
283 saveRDS( copyNumbersSmooth, robjCopyNrPath );
284 #printIgvFile( copyNumbersSmooth, igvCopyNrName )
285 #exportBins(copyNumbersSmooth, file=igvCopyNrName, format="igv")
286
287 ## also save objects for galaxy history output if requested
288 if ( doOutputReadcountsRds ){
289 saveRDS( readCountsFiltered, readCountsDatasetFile );
290 }
291 if ( doOutputCopynumbersRds ){
292 saveRDS( copyNumbersSmooth, copyNumbersDatasetFile );
293 }
294
295 ## proceed with calling if requested
296 if ( doCall ){
297 copyNumbersSegmented <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
298 copyNumbersSegmented <- normalizeSegmentedBins( copyNumbersSegmented )
299 copyNumbersCalled <- callBins( copyNumbersSegmented )
300 cgh <- makeCgh( copyNumbersCalled )
301 saveRDS( copyNumbersCalled, robjCalledPath );
302 if ( doOutputCallsRds ){
303 saveRDS( copyNumbersCalled, calledSegmentsDatasetFile );
304 }
305 #df <- getIGVdf( copyNumbersCalled )
306 #printIgvFile( copyNumbersCalled, igvCalledName )
307 #write.table( df, file=igvCalledName, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )
308 }
309
310 sampleNames <- readCountsFiltered@phenoData@data$name
311
312 ## ------------------------
313 ## create output files
314 ## ------------------------
315 plotted_images <- list() # to keep track of images for later linking
316 regions <- list() # will contain the (called) segments
317
318 noise_img_file <- paste( binSize, 'kbp_QDNAseqNoisePlot.png', sep='')
319 noise_img_file_path <- paste( outputPath, '/', noise_img_file, sep='' )
320 png( noise_img_file_path, width=480, height=480 );
321 noisePlot( readCountsFiltered )
322 dev.off()
323
324 for (i in 1:length(sampleNames) ){
325 #for (sample in sampleNames(copyNumbersSmooth) ){
326 sample <- sampleNames[i]
327 usedReads <- readCountsFiltered@phenoData@data$used.reads[i]
328 catMsg( c("Creating plots for sample: ", sample ) )
329
330 type <- 'CopyNumbers'
331 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
332 img_file_path <- paste( outputPath, '/', img_file, sep='' )
333 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT ); plot( copyNumbersSmooth[ ,sample ] ); dev.off()
334 plotted_images[[ sample ]][[ type ]] <- img_file
335
336 if ( doCall ){
337 type <- 'Called'
338 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
339 img_file_path <- paste( outputPath, '/', img_file, sep='' )
340 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT );
341 plot( copyNumbersCalled[ ,sample ] );
342 dev.off()
343 plotted_images[[ sample ]][[ type ]] <- img_file
344
345 cat( MAIN_NAME, " Fusing regions of sample: ", sample, "\n", sep="")
346 regions[[ sample ]] <- fuse.regions_test( cgh[, sample] )
347 region_count <- nrow( data.frame( regions[[ sample ]] ) )
348 cat( MAIN_NAME, ' sample "', sample, '" has ', region_count, " regions\n", sep="" )
349 plotted_images[[ sample ]][[ 'region_count' ]] <- region_count
350 }
351
352 ## add read counts
353 catMsg( "Used")
354 catMsg( usedReads )
355
356 plotted_images[[ sample ]][[ 'usedReads' ]] <- usedReads
357 }
358 #cat( MAIN_NAME, "PLOTTED_IMAGES: ", names(plotted_images), "\n", sep="" )
359
360 if ( doCall ){
361 saveRDS( regions, robjRegionPath )
362 printedFiles <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath )
363 }
364
365 ## ------------------------
366 ## prepare output
367 ## ------------------------
368 cat( MAIN_NAME, "...zipping output\n")
369 zip_cmd <- paste( "zip -j", gzipOutputPath, paste(outputPath,'/*',sep='') ) ## -j is for removing dirs from the tree
370 system( zip_cmd )
371
372 ## ------------------------
373 ## get filesizes for report
374 ## ------------------------
375 zippedSize <- paste( round( file.info( gzipOutputPath )[["size"]] / 1000000, digits=2 ), 'MB' )
376 readCoSize <- paste( round( file.info( robjReadCoPath )[["size"]] / 1000000, digits=2 ), 'MB' )
377 copyNrSize <- paste( round( file.info( robjCopyNrPath )[["size"]] / 1000000, digits=2 ), 'MB' )
378 calledSize <- paste( round( file.info( robjCalledPath )[["size"]] / 1000000, digits=2 ), 'MB' )
379 regionSize <- paste( round( file.info( robjRegionPath )[["size"]] / 1000000, digits=2 ), 'MB' )
380
381 ## ------------------------
382 ## creating html output to be linked to from the middle galaxy pane
383 ## ------------------------
384
385 #sink( file = outputHtml, type = "output" )
386 sink( file = htmlOutputPath, type = "output" )
387 cat( "<html>\n")
388 cat( "<head>\n")
389 #cat( '<link rel="stylesheet" type="text/css" href="test.css" media="screen" />', "\n" )
390 #cat( '<link rel="stylesheet" type="text/css" href="../../../../static/style/test.css" media="screen" />',
391 cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
392 cat( "\t<style>\n")
393 ## have to include CSS into html file, because css referencing outside own dir doesn't seem to work...
394 cat( paste( "\t\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
395 cat( paste( "\t\t", readLines( CSS_FILE, n = -1)), sep="\n" )
396 #cat( "\t\th1 {color:red;}", "\n")
397
398 cat( "\t</style>\n" )
399
400 cat( "\n</head>\n")
401 cat( "\n<body>\n")
402
403 cat( "<h1>QDNAseq Report</h1>", "\n")
404
405 cat( '<h3 class="qdnaseq">About this analysis</h3>', "\n")
406 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='')
407 #cat( '<a href="#" class="button-success button-xsmall pure-button">test</a>' )
408
409
410 ## ------------------------
411 ## table with general info
412 ## ------------------------
413 cat( '<h3 class="qdnaseq">Settings</h3><p>', "\n")
414 cat( '<table class="pure-table pure-table-striped"><thead><tr><th>setting</th><th>value</th></tr></thead><tbody>' )
415 cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
416 cat( htmlTableRow( c( "AnalysisDate", 'todo' ) ) )
417 cat( htmlTableRow( c( "BinSize (kb)", binSize ) ) )
418 #cat( htmlTableRow( c( "undoSD", undoSD ) ) )
419 #cat( htmlTableRow( c( "useBlacklist", filterBlacklistedBins ) ) )
420
421 for ( galaxyName in fileList[[ 'all_files' ]] ){
422 sampleName <- fileList[[ galaxyName ]]
423 cat( htmlTableRow( c( "InputBam", paste( galaxyName, ' (', sampleName, ')', sep='' ) ) ) )
424 }
425 cat( "</tbody></table></p>", "\n")
426
427 ## ------------------------
428 ## list with links to all output files
429 ## ------------------------
430 r_code <- '<p>'
431 r_code <- paste( r_code, '<code class="code">## R code to load files</code><br />', sep="\n" )
432 r_code <- paste( r_code, '<code class="code">library( QDNAseq )</code><br />', sep="\n")
433
434 cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
435 cat( '<dl>', "\n" )
436 #cat( '<dt>Definition term</dt>', "\n", '<dd>Definition term</dd>', "\n" )
437 cat( '<dt>', htmlLink( path=robjReadCoName, robjReadCoName ), '</dt>', "\n" )
438 cat( '<dd>QDNAseq object with read counts per bin, ', readCoSize,'</dd>', "\n" )
439 r_code <- paste( r_code, '<code class="code">readCounts <- readRDS(', robjReadCoName, ")</code><br />", sep="")
440
441 cat( '<dt>', htmlLink( path=robjCopyNrName, robjCopyNrName ), '</dt>', "\n" )
442 cat( '<dd>QDNAseq object with copy numbers per bin, ', copyNrSize,'</dd>', "\n" )
443 r_code <- paste( r_code, '<code class="code">copyNumbersSmooth <- readRDS(', robjCopyNrName, ")</code><br />", sep="")
444
445 if ( doCall ){
446 cat( '<dt>', htmlLink( path=robjCalledName, robjCalledName ), '</dt>', "\n" )
447 cat( '<dd>QDNAseq object with segment and call status per bin, ', calledSize,'</dd>', "\n" )
448 r_code <- paste( r_code, '<code class="code">copyNumbersCalled <- readRDS(', robjCalledName, ")</code><br />", sep="")
449
450 cat( '<dt>', htmlLink( path=regiOutputName, regiOutputName ), '</dt>', "\n" )
451 cat( '<dd>list with segmented/called regions for each sample, ', regionSize, '</dd>', "\n" )
452 r_code <- paste( r_code, '<code class="code">calledRegions <- readRDS(', regiOutputName, ")</code><br />", sep="")
453
454 }
455 cat( '</dl></p>', "\n" )
456
457 cat( r_code, "</p>\n", sep="\n")
458 cat( '<p>See ', htmlLink( WEB_LINK, 'the bioconductor QDNAseq documentation' ), ' for more information on how to work with these files</p>', "\n", sep='' )
459
460 ## ------------------------
461 ## table with links to files
462 ## ------------------------
463 cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
464 plots_html <- ''
465
466 cat( '<table class="pure-table pure-table-striped"><thead><tr><th>Sample</th><th>CopyNumber</th><th>Called</th><th>ReadCount</th><th>RegionCount</th><th>Files</th></tr></thead><tbody>' )
467
468 for ( bam_file in sampleNames ){
469
470 #width <- 600; height <- 240
471 width <- PLOT_WIDTH; height <- PLOT_HEIGHT
472 width_t <- 100; height_t <- 40
473
474 ## add thumbnails to table with links to anchors on html page
475 copy_img <- plotted_images[[ bam_file ]][[ 'CopyNumbers' ]]
476 usedReads <- plotted_images[[ bam_file ]][[ 'usedReads' ]]
477 usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )
478
479 html_copy_thumb <- htmlLink( path=paste('#', copy_img, sep=''), paste('<img src="',copy_img,'" alt="', bam_file, '" width="', width_t, '" height="', height_t, '">', sep='') )
480 html_copy_img <- htmlLink( path=copy_img, paste('<img id="', copy_img,'" src="',copy_img,'" alt="',bam_file, '" width="', width, '" height="', height, '">', sep='') )
481 html_call_thumb <- 'NA'
482 html_call_img <- ''
483 html_bed <- 'NA'
484 html_bedGraph <- 'NA'
485 region_count <- 'NA'
486
487 if ( doCall ){
488 call_img <- plotted_images[[ bam_file ]][[ 'Called' ]]
489 region_count <- plotted_images[[ bam_file ]][[ 'region_count' ]]
490 html_call_thumb <- htmlLink( path=paste('#', call_img, sep=''), paste('<img src="', call_img, '" alt="', bam_file, '" width="', width_t,'" height="', height_t,'">', sep='') )
491
492 files <- printedFiles[[ bam_file ]]
493 html_bed <- htmlLink( path=files[1], 'bed' )
494 html_bedGraph <- htmlLink( path=files[2], 'bedGraph' )
495 html_call_img <- htmlLink( path=call_img, paste('<img id="', call_img,'" src="', call_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
496 }
497
498 ## add info to overview table, including small thumbnails
499 cat( htmlTableRow( c(bam_file, html_copy_thumb, html_call_thumb, usedReads, region_count, paste( html_bed, html_bedGraph, sep=", ")) ), "\n" )
500 ## now include (large) images in html page
501 plots_html <- paste( plots_html, html_copy_img, "\n", html_call_img, "\n<hr \\>\n", sep='' )
502 }
503 cat( "</tbody></table></p>", "\n")
504
505 ## add noise plot
506 html_noise_img <- htmlLink( path=noise_img_file, paste('<img id="', noise_img_file,'" src="',noise_img_file,'" alt="NoisePlot">', sep='') )
507 plots_html <- paste( plots_html, html_noise_img, "\n<hr \\>\n", sep='' )
508
509 ## ------------------------
510 ## section with various output shown
511 ## ------------------------
512 cat( '<h3 class="qdnaseq">Results: plots</h3><p>', "\n")
513 cat( plots_html, "\n")
514 cat( "\n</p></body>\n")
515 cat( "\n</html>\n")
516 sink()
517
518 ## ------------------------
519 ## creating html output to be viewed in middle galaxy pane
520 ## ------------------------
521 #sink( file = htmlOutputPath, type = "output" )
522 sink( file = outputHtml, type = "output" )
523
524 cat( "<head>", "\n")
525 cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
526
527 cat( "<style>", "\n")
528 ## have to include CSS into html file, because css referencing outside own dir doesn't seem to work...makes it more portable anyway :P
529 cat( paste( "\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
530 cat( paste( "\t", readLines( CSS_FILE, n = -1)), sep="\n" )
531 cat( "</style>", "\n")
532 cat( "</head>", "\n")
533
534
535 cat( '<h1>QDNAseq Results (', outputName,')</h1>', "\n", sep="")
536 cat( '<p>Explore <a href="', htmlOutputName, '" class="button-success button-small pure-button">the results</a> directly within galaxy</p>', "\n", sep="")
537 cat( '<p>Or download a <a href="', gzipOutputName, '" class="button-success button-small pure-button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )
538
539 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="")
540
541 sink()
542
543
544 cat( MAIN_NAME, "...zipping output\n")
545 zip_cmd <- paste( "zip -j ", gzipOutputPath, paste(outputPath,'/', htmlOutputName, sep='') ) ## -j is for removing dirs from the tree
546 system( zip_cmd )
547 cat( MAIN_NAME, "done...\n", sep="")
548 q(status=0)