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