|
6
|
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 location of this script
|
|
|
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 creation 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 input 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 if ( l1 == 0 ){ return(tmp) } # empty list in debug mode
|
|
|
42
|
|
|
43 for ( i in 1:length(paths) ){
|
|
|
44 path <- paths[i]
|
|
|
45 name <- names[i]
|
|
|
46 file <- basename(path)
|
|
|
47
|
|
|
48 tmp[[ file ]] <- name
|
|
|
49 tmp[[ 'all_paths' ]] <- c( tmp[[ 'all_paths' ]], path )
|
|
|
50 tmp[[ 'all_files' ]] <- c( tmp[[ 'all_files' ]], file )
|
|
|
51 tmp[[ 'all_names' ]] <- c( tmp[[ 'all_names' ]], name )
|
|
|
52 }
|
|
|
53 return( tmp )
|
|
|
54 }
|
|
|
55
|
|
|
56 ## --------------------
|
|
|
57 ## copied code for extracting the regions by segment call status
|
|
|
58 ## --------------------
|
|
|
59 fuseRegions <- function( obj, minRatio=0 ) {
|
|
|
60 if ( ncol(obj) > 1 ) stop('Please specify which sample...')
|
|
|
61
|
|
|
62 data <- data.frame( obj@featureData@data[,1:3], copynumber(obj), segmented(obj), check.names=FALSE, stringsAsFactors=FALSE)
|
|
|
63 colnames( data ) <- c( "chr", "start", "end", "log2", "segmentval" )
|
|
|
64
|
|
|
65 fused.data <- data.frame()
|
|
|
66 curr.bin <- 1
|
|
|
67 for ( chr in unique( data$chr ) ) {
|
|
|
68 chr.data <- data[ data$chr == chr, ]
|
|
|
69 prev.bin <- curr.bin
|
|
|
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 curr.bin <- curr.bin + 1
|
|
|
77 curr.segm <- chr.data[ i, 'segmentval']
|
|
|
78
|
|
|
79 if ( curr.segm != prev.segm ) {
|
|
|
80 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
81 prev.segm <- curr.segm
|
|
|
82 prev.bin <- curr.bin
|
|
|
83 start <- chr.data[ i, 'start']
|
|
|
84 }
|
|
|
85 }
|
|
|
86 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
87 }else{
|
|
|
88 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
89 }
|
|
|
90 }
|
|
|
91 ## remove regions with low amplitude
|
|
|
92 fused.data <- fused.data[ abs(fused.data$segmentval) >= minRatio, ]
|
|
|
93 fused.data
|
|
|
94 }
|
|
|
95
|
|
|
96 ## DESC: takes the output of fuse.regions and outputs a txt file per sample
|
|
|
97 outputRegionsFromList <- function ( regionsList, outputBasename, outputDir="./", binSize, storeList ){
|
|
|
98 if ( missing(regionsList) ) stop( 'Please provide regionsList...' )
|
|
|
99 if ( missing(outputBasename) ) stop( 'Please provide outputBasename...' )
|
|
|
100 if ( !is.list(regionsList) ) stop( 'Input not a list...?' )
|
|
|
101 if ( length(regionsList) < 1 ) stop( 'List seems empty...?' )
|
|
|
102 if ( file.exists( outputDir ) ) catMsg( c(" Using dir ", outputDir, " for output") )
|
|
|
103 else dir.create( outputDir )
|
|
|
104
|
|
|
105 ## have to set R output options otherwise scientific method is used at some point
|
|
|
106 options( "scipen"=100 )
|
|
|
107
|
|
|
108 sampleCount <- length( regionsList )
|
|
|
109 sampleNames <- names( regionsList )
|
|
|
110 bedgraphColumns <- c( 'chr', 'start', 'end', 'segmentval' )
|
|
|
111
|
|
|
112 catMsg( c( " There are ", sampleCount, " samples found in input list") )
|
|
|
113
|
|
|
114 for ( sample in sampleNames ){
|
|
|
115 catMsg( c(" Working on sample ", sample ) )
|
|
|
116 regionCount <- nrow( regionsList[[sample]] )
|
|
|
117
|
|
|
118 outSampleBase <- paste( outputBasename, '_', sample, '_', binSize, 'kbp', sep='')
|
|
|
119 outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
|
|
|
120 outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )
|
|
|
121
|
|
|
122 ## ---------- BEDGRAPH ----------
|
|
|
123 txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample," description=segmented_regions_from_QDNAseq_",binSize,"kbp\n", sep="")
|
|
|
124 sink( outBedgraphPath )
|
|
|
125 cat( txt )
|
|
|
126 sink()
|
|
|
127 write.table( regionsList[[sample]][,bedgraphColumns], outBedgraphPath, quote=F, sep="\t", row.names=F, append=T, col.names=F)
|
|
|
128 #outFiles[[sample]] <- c( outBedgraphFile )
|
|
|
129 storeList[[ paste( binSize, sample, 'bedgraph', sep="_")]] <- outBedgraphFile
|
|
|
130 }
|
|
|
131 return(storeList)
|
|
|
132 }
|
|
|
133
|
|
|
134
|
|
|
135 ## ==================================================
|
|
|
136 ## Start of analysis
|
|
|
137 ## ==================================================
|
|
|
138 MAIN_NAME <- '[INFO] '
|
|
|
139 TOOL_PATH <- getScriptPath()
|
|
|
140 CSS_FILE <- paste( TOOL_PATH, '/static/css/QDNAseq.css', sep="" )
|
|
|
141 DECIMALS <- 3
|
|
|
142 WEB_LINK <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
|
|
|
143 PURE_CSS <- 'http://yui.yahooapis.com/pure/0.5.0/pure-min.css'
|
|
|
144
|
|
|
145 catMsg( "Starting QDNAseq wrapper" )
|
|
|
146 catMsg( "Loading R libraries" )
|
|
7
|
147 catMsg( R.version.string )
|
|
|
148
|
|
8
|
149 packs <- c( "limma", "marray", "CGHbase", "impute", "DNAcopy", "CGHcall" )
|
|
|
150 for ( pack in packs ){
|
|
|
151 catMsg( c("Loading ", pack) )
|
|
|
152 library( pack )
|
|
|
153 }
|
|
7
|
154
|
|
|
155 ## supress msg to allow R to finish with non-error msg
|
|
6
|
156 suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
|
|
|
157 suppressWarnings( suppressMessages( library( CGHcall, quietly = TRUE ) ) )
|
|
|
158
|
|
|
159 ## only one param: the tmp config file
|
|
|
160 cmdLineArgs <- commandArgs(TRUE)
|
|
|
161 config <- cmdLineArgs[1]
|
|
|
162
|
|
|
163 ## sourcing the config file will load all input params
|
|
|
164 ## many variables are imported via sourced "config"
|
|
|
165 source( config )
|
|
|
166
|
|
|
167 ## desparate tries to make png text scale well, damn you R...!
|
|
|
168 PLOT_RES <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 6.3
|
|
|
169 PAR_SET <- list( pch=22 )
|
|
|
170
|
|
|
171 systemUser <- system("whoami",T)
|
|
|
172 qdnaseqVersion <- packageDescription( "QDNAseq" )$Version
|
|
|
173 rVersion <- R.version.string
|
|
|
174 startTime <- Sys.time()
|
|
|
175 analysisStart <- as.character( startTime )
|
|
|
176 catMsg( c("QDNAseq version: ", qdnaseqVersion) )
|
|
|
177 catMsg( c( rVersion ) )
|
|
|
178
|
|
|
179 ## get the comma separated list of chromosomes to exclude
|
|
|
180 excludeChrs <- unlist( strsplit( excludeChrsString, ",") )
|
|
|
181 binSizes <- as.numeric( unlist( strsplit( binSizesString, ",") ) )
|
|
|
182
|
|
|
183
|
|
|
184 ## ------------------------
|
|
|
185 ## DEBUG
|
|
|
186 if ( debug ){
|
|
|
187 catMsg( c("Analysis run by user: ", systemUser ) )
|
|
|
188 catMsg( c("DEBUG SessionInfo: " ) )
|
|
|
189 print( sessionInfo() )
|
|
|
190 }
|
|
|
191 ## /DEBUG
|
|
|
192 ## ------------------------
|
|
|
193
|
|
|
194 ## prepare output dir
|
|
|
195 if ( !file.exists( outputPath) ){
|
|
|
196 dir.create( outputPath )
|
|
|
197 }
|
|
|
198
|
|
|
199 ## copy source config file to output dir to include it in output zip
|
|
|
200 if ( inGalaxy ){
|
|
|
201 file.copy( config, paste(outputPath, 'galaxyConfigFile.R', sep='/') )
|
|
|
202 }
|
|
|
203
|
|
|
204 ## setup bam filelist for easy retrieval later
|
|
|
205 fileList <- makeBamFileList( bamsPaths, bamsNames )
|
|
|
206 bamCount <- length( fileList[[ 'all_paths' ]] )
|
|
|
207
|
|
|
208 gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
|
|
|
209 gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep='')
|
|
|
210 htmlOutputName <- 'index.html'
|
|
|
211 htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep='')
|
|
|
212
|
|
|
213 plotted_images <- list() # to keep track of images for later linking
|
|
|
214 regions <- list() # will contain the segments
|
|
|
215 outputFiles <- list()
|
|
|
216
|
|
|
217 ## ------------------------
|
|
|
218 ## in case of debug just use inbuilt LGG data for speedup
|
|
|
219 if ( debug ){
|
|
|
220 binSizes <- c(15)
|
|
|
221 bamsPaths <- c( "BUILD_IN_DATA")
|
|
|
222 bamsNames <- c( "LGG150")
|
|
|
223 fileList <- makeBamFileList( bamsPaths, bamsNames )
|
|
|
224 bamCount <- length( fileList[[ 'all_paths' ]] )
|
|
|
225 }
|
|
|
226
|
|
|
227 for ( binSize in binSizes ){
|
|
|
228
|
|
|
229 ## ------------------------
|
|
|
230 ## construct output file-names and -paths
|
|
|
231 ## ------------------------
|
|
|
232 robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
|
|
|
233 robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
|
|
|
234 igvCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.igv', sep='')
|
|
|
235 robjSegmntName <- paste( binSize, 'kbp_QDNAseqCopyNumbersSegmented.rds', sep='')
|
|
|
236 regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
|
|
|
237 noiseImgName <- paste( binSize, 'kbp_QDNAseqNoiseplot.png', sep='')
|
|
|
238
|
|
|
239 robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep='')
|
|
|
240 robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep='')
|
|
|
241 robjSegmntPath <- paste( outputPath, '/', robjSegmntName, sep='')
|
|
|
242 robjRegionPath <- paste( outputPath, '/', regiOutputName, sep='')
|
|
|
243 igvCopyNrPath <- paste( outputPath, '/', igvCopyNrName, sep='')
|
|
|
244 noiseImgPath <- paste( outputPath, '/', noiseImgName, sep='')
|
|
|
245
|
|
|
246 binAnnFile <- paste( TOOL_PATH, '/static/binannotation/', binSize, 'kbp_binAnnotations.rds', sep="" )
|
|
|
247 if ( file.exists(binAnnFile) ){
|
|
|
248 binAnnotations <- readRDS( binAnnFile )
|
|
|
249 catMsg( c("Using local binAnnotations file" ) )
|
|
|
250 }else{
|
|
|
251 binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
|
|
|
252 }
|
|
|
253
|
|
|
254 ## in case of debug just use inbuilt LGG data for speedup
|
|
|
255 if ( debug ){
|
|
|
256 data(LGG150)
|
|
|
257 readCounts <- LGG150
|
|
|
258 }else{
|
|
|
259 ## provide bamnames because in galaxy everyting is called "dataset_###"
|
|
|
260 readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=fileList[[ 'all_names' ]] )
|
|
|
261 }
|
|
|
262
|
|
|
263 readCountsFiltered <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
|
|
|
264 readCountsFiltered <- estimateCorrection( readCountsFiltered )
|
|
|
265 copyNumbers <- correctBins( readCountsFiltered )
|
|
|
266 copyNumbersNormalized <- normalizeBins( copyNumbers )
|
|
|
267 copyNumbersSmooth <- smoothOutlierBins( copyNumbersNormalized )
|
|
|
268 sampleNames <- readCountsFiltered@phenoData@data$name
|
|
|
269
|
|
|
270 ## save objects to output dir
|
|
|
271 saveRDS( readCountsFiltered, robjReadCoPath );
|
|
|
272 saveRDS( copyNumbersSmooth, robjCopyNrPath );
|
|
|
273 exportBins( copyNumbersSmooth, file=igvCopyNrPath, format="igv" )
|
|
|
274
|
|
|
275 ## also save objects for galaxy history output if requested
|
|
|
276 if ( doOutputCopynumbersIgv ){
|
|
|
277 #@ a bit hacky galaxy way to allow an unknown number of output files based on param selection
|
|
|
278 #@ see: https://wiki.galaxyproject.org/Admin/Tools/Multiple%20Output%20Files
|
|
|
279 historyName <- paste(binSize, 'kbp-IGV', sep="")
|
|
|
280 igvFile <- paste( newFilePath, "/primary_", outputId, "_", historyName, "_visible_txt", sep="" )
|
|
|
281 exportBins( copyNumbersSmooth, file=igvFile, format="igv" )
|
|
|
282 catMsg( c("Exported igv file to history for ", binSize, "kbp bin") )
|
|
|
283 }
|
|
|
284
|
|
|
285 ## proceed with calling if requested
|
|
|
286 if ( doSegment ){
|
|
|
287 copyNumbersSegmented <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
|
|
|
288 copyNumbersSegmented <- normalizeSegmentedBins( copyNumbersSegmented )
|
|
|
289 cgh <- makeCgh( copyNumbersSegmented )
|
|
|
290 saveRDS( copyNumbersSegmented, robjSegmntPath );
|
|
|
291 }
|
|
|
292
|
|
|
293 ## ------------------------
|
|
|
294 ## create output files
|
|
|
295 ## ------------------------
|
|
|
296 png( noiseImgPath, width=PLOT_HEIGHT, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
297 par( PAR_SET )
|
|
|
298 noisePlot( readCountsFiltered, main=paste( "Noise Plot ", binSize, "kbp", sep=''), col="darkgreen" )
|
|
|
299 dev.off()
|
|
|
300
|
|
|
301 binSize <- as.character( binSize ) # to avoid R using it as array index... (*#$^@ you R!)
|
|
|
302 binSizeString <- paste( binSize, 'kbp', sep='')
|
|
|
303
|
|
|
304 for (i in 1:length(sampleNames) ){
|
|
|
305
|
|
|
306 sample <- sampleNames[i]
|
|
|
307 usedReads <- readCountsFiltered@phenoData@data$used.reads[i]
|
|
|
308 catMsg( c("Creating plots for sample: ", sample, " (", binSizeString, ")" ) )
|
|
|
309
|
|
|
310 type <- 'CopyNumbers'
|
|
|
311 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
|
|
|
312 img_file_path <- paste( outputPath, '/', img_file, sep='' )
|
|
|
313
|
|
|
314 ## COPYNUMBER PLOT
|
|
|
315 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
316 par( PAR_SET )
|
|
|
317 plot( copyNumbersSmooth[ ,sample ], main=paste(sample, ": CopyNumbers", sep="") )
|
|
|
318 mtext( paste( binSizeString, " bins", sep=""), 3 )
|
|
|
319 abline( h=c(-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
|
|
|
320 dev.off()
|
|
|
321
|
|
|
322 plotted_images[[ paste(binSize, sample, type, sep="_" ) ]] <- img_file
|
|
|
323
|
|
|
324 if ( doSegment ){
|
|
|
325 type <- 'Segmented'
|
|
|
326 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
|
|
|
327 img_file_path <- paste( outputPath, '/', img_file, sep='' )
|
|
|
328
|
|
|
329 ## COPYNUMBER PLOT
|
|
|
330 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
331 par( PAR_SET )
|
|
|
332 plot( copyNumbersSegmented[ ,sample ], main=paste(sample, ": CopyNumbers and Segments", sep="") )
|
|
|
333 mtext( paste( "(", binSizeString, " bins)", sep=""), 3 )
|
|
|
334 abline( h=c(-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
|
|
|
335 dev.off()
|
|
|
336
|
|
|
337 plotted_images[[ paste(binSize, sample, type, sep="_" ) ]] <- img_file
|
|
|
338
|
|
|
339 catMsg( c(" Fusing regions of sample: ", sample) )
|
|
|
340 regions[[ sample ]] <- fuseRegions( cgh[, sample] )
|
|
|
341
|
|
|
342 region_count <- nrow( data.frame( regions[[ sample ]] ) )
|
|
|
343 catMsg( c( ' sample "', sample, '" has ', region_count, " regions" ) )
|
|
|
344 plotted_images[[ paste(binSize, sample, 'region_count', sep="_" ) ]] <- region_count
|
|
|
345 }
|
|
|
346
|
|
|
347 ## add USED read counts
|
|
|
348 plotted_images[[ paste(binSize, sample, 'usedReads', sep="_" ) ]] <- usedReads
|
|
|
349 }
|
|
|
350
|
|
|
351 if ( doSegment ){
|
|
|
352 saveRDS( regions, robjRegionPath )
|
|
|
353 plotted_images <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath, binSize=binSize, storeList=plotted_images )
|
|
|
354 }
|
|
|
355 }# end bin
|
|
|
356
|
|
|
357
|
|
|
358 ## ----- debug -----
|
|
|
359 #catMsg( "done" )
|
|
|
360 #q(status=0)
|
|
|
361 ## ---- /debug -----
|
|
|
362
|
|
|
363
|
|
|
364 ## ------------------------
|
|
|
365 ## prepare output
|
|
|
366 ## ------------------------
|
|
|
367 catMsg( "...zipping output")
|
|
|
368 zip_cmd <- paste( "zip -j", gzipOutputPath, paste(outputPath,'/*',sep='') ) ## -j is for removing dirs from the tree
|
|
|
369 system( zip_cmd )
|
|
|
370
|
|
|
371 ## ------------------------
|
|
|
372 ## get filesizes for report
|
|
|
373 ## ------------------------
|
|
|
374 zippedSize <- paste( round( file.info( gzipOutputPath )[["size"]] / 1e+6, digits=2 ), 'MB' )
|
|
|
375 endTime <- Sys.time()
|
|
|
376 timeDiff <- format( round( endTime - startTime, 3 ) )
|
|
|
377 analysisEnd <- as.character( endTime )
|
|
|
378
|
|
|
379 ## ------------------------
|
|
|
380 ## creating html output to be linked to from the middle galaxy pane
|
|
|
381 ## ------------------------
|
|
|
382 sink( file = htmlOutputPath, type = "output" )
|
|
|
383 cat( "<html>\n")
|
|
|
384 cat( "<head>\n")
|
|
|
385
|
|
|
386 cat( "\t", '<title>QDNAseq Report | ', outputName,'</title>', "\n", sep='' )
|
|
|
387 cat( "\t", '<link rel="stylesheet" href="', PURE_CSS, '">', "\n", sep='' )
|
|
|
388 cat( "\t<style>\n", sep='')
|
|
|
389 ## include CSS into html file, makes it more portable
|
|
|
390 cat( "\t\t", readLines( CSS_FILE ), sep="\n\t\t" )
|
|
|
391 #cat( "\t\th1 {color:red;}", "\n")
|
|
|
392 cat( "\n\t</style>\n" )
|
|
|
393
|
|
|
394 cat( "\n</head>\n")
|
|
|
395 cat( "\n<body>\n")
|
|
|
396
|
|
|
397 cat( "<h1>QDNAseq Report</h1>", "\n")
|
|
|
398
|
|
|
399 cat( '<h3 class="qdnaseq">About this analysis</h3>', "\n")
|
|
|
400 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">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep='')
|
|
|
401
|
|
|
402 ## ------------------------
|
|
|
403 ## table with general info
|
|
|
404 ## ------------------------
|
|
|
405 cat( '<h3 class="qdnaseq">Settings</h3><p>', "\n")
|
|
|
406 cat( '<table class="pure-table pure-table-striped"><tbody>' )
|
|
|
407 cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
|
|
|
408 cat( htmlTableRow( c( "AnalysisStart", analysisStart ) ) )
|
|
|
409 cat( htmlTableRow( c( "AnalysisEnd", analysisEnd ) ) )
|
|
|
410 cat( htmlTableRow( c( "AnalysisTime", timeDiff ) ) )
|
|
|
411 cat( htmlTableRow( c( "BinSizes (kbp)", paste(binSizes,collapse=", ") ) ) )
|
|
|
412 cat( htmlTableRow( c( "R info", rVersion ) ) )
|
|
|
413 cat( htmlTableRow( c( "QDNAseq info", qdnaseqVersion ) ) )
|
|
|
414
|
|
|
415 sampleStrings <- c()
|
|
|
416 for ( galaxyName in fileList[[ 'all_files' ]] ){
|
|
|
417 sampleName <- fileList[[ galaxyName ]]
|
|
|
418 sampleStrings <- c( sampleStrings, paste( galaxyName, ' (', sampleName, ')', sep='' ) )
|
|
|
419 }
|
|
|
420 cat( htmlTableRow( c( "InputBams", paste( sampleStrings, collapse=", ") ) ) )
|
|
|
421
|
|
|
422 cat( "</tbody></table></p>", "\n")
|
|
|
423
|
|
|
424 ## ------------------------
|
|
|
425 ## list with links to all output files
|
|
|
426 ## ------------------------
|
|
|
427 cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
|
|
|
428 cat( '<p>This table contains output files that can be used for local downstream analysis with the bioconductor QDNAseq package. For each bin-size / data-level there is a R data structure file with data of all samples. See ', htmlLink( WEB_LINK, 'the bioconductor QDNAseq documentation' ), ' for more information on how to work with these files</p>', "\n", sep='')
|
|
|
429 cat( '<table class="pure-table pure-table-striped">', "\n" )
|
|
|
430 cat( '<thead><th>Type</th>', as.vector( mapply( paste, "<th>", binSizes, "kbp</th>", sep="" ) ),'</thead>', "\n" )
|
|
|
431 cat( "<tbody>", "\n")
|
|
|
432 files <- list()
|
|
|
433 fileTypes <- c( 'ReadCounts.rds', 'CopyNumbers.rds' )
|
|
|
434 if ( doSegment ){ fileTypes <- c( fileTypes, 'CopyNumbersSegmented.rds') }
|
|
|
435
|
|
|
436 for ( fileType in fileTypes ){
|
|
|
437 fileNames <- mapply( paste, binSizes, paste( 'kbp_QDNAseq', fileType, sep=''), sep='')
|
|
|
438 fileLinks <- mapply( htmlLink, fileNames, paste( binSizes, "kbp", sep="" ) )
|
|
|
439 cat( htmlTableRow( c( fileType, fileLinks ) ) )
|
|
|
440 }
|
|
|
441 cat( "\n</tbody></table></p>", "\n")
|
|
|
442
|
|
|
443 ## ------------------------
|
|
|
444 ## table with links to files
|
|
|
445 ## ------------------------
|
|
|
446 ratio <- PLOT_WIDTH / PLOT_HEIGHT
|
|
|
447 width <- 960; height <- width / ratio ## bigger img
|
|
|
448 width_t <- 100; height_t <- 40 ## thumb img
|
|
|
449
|
|
|
450 cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
|
|
|
451 cat( '<p>This table contains the visual results of the copy number aberration analysis. You can click on an image to jump to the larger version. If segmentation was performed as well the number of segments is shown and a file with genomic regions can be downloaded (just remember to inspect the results carefully as this is a more exprimental step).</p>', "\n", sep='')
|
|
|
452 plots_html <- ''
|
|
|
453
|
|
|
454 colspan <- 1
|
|
|
455 binHeader <- "<th>Image</th>"
|
|
|
456 if ( doSegment ){ # extra column with segment info
|
|
|
457 colspan <- 2
|
|
|
458 binHeader <- "<th>Image</th><th>Segments</th>"
|
|
|
459 }
|
|
|
460 cat( '<table class="pure-table pure-table-striped">', "\n" )
|
|
|
461 cat( '<thead><tr><th></th><th></th>', as.vector( mapply( paste, "<th colspan=\"", colspan,"\">", binSizes, "kbp</th>", sep="" ) ), '</tr></thead>' )
|
|
|
462 cat( '<thead><tr><th>Sample / File</th><th>Reads</th>', rep( binHeader, length(binSizes) ), '</tr></thead>' )
|
|
|
463 cat( '<tbody>' )
|
|
|
464
|
|
|
465 for ( bam_file in bamsNames ){
|
|
|
466
|
|
|
467 usedReads <- plotted_images[[ paste(binSize, bam_file, 'usedReads', sep="_" ) ]]
|
|
|
468 usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )
|
|
|
469 htmlRow <- paste( '<tr><td>', bam_file, '</td><td>', usedReads, '</td>', sep='' )
|
|
|
470
|
|
|
471 for ( binSize in binSizes ){
|
|
|
472
|
|
|
473 ## add thumbnails to table with links to anchors on html page
|
|
|
474 copy_img <- plotted_images[[ paste(binSize, bam_file, 'CopyNumbers', sep="_" ) ]]
|
|
|
475 html_copy_thumb <- htmlLink( path=paste('#', copy_img, sep=''), paste('<img src="',copy_img,'" alt="', bam_file, '" width="', width_t, '" height="', height_t, '">', sep='') )
|
|
|
476 html_copy_img <- htmlLink( path=copy_img, paste('<img id="', copy_img,'" src="',copy_img,'" alt="',bam_file, '" width="', width, '" height="', height, '">', sep='') )
|
|
|
477 html_segm_img <- ''
|
|
|
478 html_bedGraph <- ''
|
|
|
479 region_count <- ''
|
|
|
480 htmlRow <- paste( htmlRow, '<td>', html_copy_thumb, '</td>' )
|
|
|
481
|
|
|
482 if ( doSegment ){
|
|
|
483 segm_img <- plotted_images[[ paste(binSize, bam_file, 'Segmented', sep="_" ) ]]
|
|
|
484 region_count <- plotted_images[[ paste(binSize, bam_file, 'region_count', sep="_" ) ]]
|
|
|
485
|
|
|
486 html_bedGraph <- htmlLink( path=plotted_images[[ paste(binSize, bam_file, 'bedgraph', sep="_" ) ]], 'bedGraph' )
|
|
|
487 html_segm_img <- htmlLink( path=segm_img, paste('<img id="', segm_img,'" src="', segm_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
|
|
|
488 htmlRow <- paste( htmlRow, '<td>', region_count, ' (', html_bedGraph, ')</td>', sep="" )
|
|
|
489 }
|
|
|
490 plots_html <- paste( plots_html, html_copy_img, "\n", html_segm_img, "\n<br \\>\n", sep='' )
|
|
|
491 }
|
|
|
492 plots_html <- paste( plots_html, "\n<hr \\>\n", sep='' )
|
|
|
493 ## add info to overview table, including small thumbnails
|
|
|
494 htmlRow <- paste( htmlRow, '</tr>', sep='' )
|
|
|
495 cat( htmlRow, "\n" )
|
|
|
496 }
|
|
|
497 cat( "</tbody></table></p>", "\n")
|
|
|
498
|
|
|
499 ## ------------------------
|
|
|
500 ## section with various output shown
|
|
|
501 ## ------------------------
|
|
|
502 cat( '<h3 class="qdnaseq">Results: Sample plots</h3><p>', "\n")
|
|
|
503 ## now include (large) images in html page
|
|
|
504 cat( plots_html, "\n")
|
|
|
505 cat( "\n</p></body>\n")
|
|
|
506 cat( "\n</html>\n")
|
|
|
507 sink()
|
|
|
508
|
|
|
509 ## ------------------------
|
|
|
510 ## creating main html output for galaxy history
|
|
|
511 ## ------------------------
|
|
|
512 if ( inGalaxy ){ # dont create when running outside Galaxy
|
|
|
513 sink( file = outputHtml, type = "output" )
|
|
|
514
|
|
|
515 cat( "<head>", "\n")
|
|
|
516 cat( "\t", '<link rel="stylesheet" href="', PURE_CSS, '">', "\n", sep='' )
|
|
|
517
|
|
|
518 cat( "<style>", "\n")
|
|
|
519 ## include CSS directly into html file
|
|
|
520 cat( paste( "\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
|
|
|
521 cat( paste( "\t", readLines( CSS_FILE, n = -1)), sep="\n" )
|
|
|
522 cat( "</style>", "\n")
|
|
|
523 cat( "</head>", "\n")
|
|
|
524
|
|
|
525 cat( '<h1>QDNAseq Results (', outputName,')</h1>', "\n", sep="")
|
|
|
526 cat( '<p>Explore <a href="', htmlOutputName, '" class="button">the results</a> directly within galaxy</p>', "\n", sep="")
|
|
|
527 cat( '<p>Or download a <a href="', gzipOutputName, '" class="button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )
|
|
|
528
|
|
|
529 sink()
|
|
|
530 }
|
|
|
531
|
|
|
532 ## ------------------------
|
|
|
533 ## create final zip and quit with status 0 to tell galaxy all was fine
|
|
|
534 ## ------------------------
|
|
|
535 catMsg( "zipping all output")
|
|
|
536 system( paste( "zip -j ", gzipOutputPath, paste(outputPath,'/', htmlOutputName, sep='') ) )
|
|
|
537 catMsg( "done" )
|
|
|
538 q(status=0)
|
