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