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