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