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