diff QDNAseq.R @ 0:cc5254dac2b6 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/QDNAseq.R	Wed Aug 27 04:59:55 2014 -0400
@@ -0,0 +1,520 @@
+#!/usr/bin/Rscript
+
+## --------------------
+## prints all arguments as msg
+## --------------------
+catMsg <- function( msg=c() ){	
+	cat( MAIN_NAME, paste( msg, collapse="" ), "\n", sep='')
+}
+## --------------------
+## return the location of this script
+## --------------------
+getScriptPath <- function(){
+    cmd.args <- commandArgs()
+    m <- regexpr("(?<=^--file=).+", cmd.args, perl=TRUE)
+    script.dir <- dirname(regmatches(cmd.args, m))
+    if( length(script.dir) == 0 ) stop("[ERR] Can't determine script dir: please call the script with Rscript\n")
+    if( length(script.dir) > 1 ) stop("[ERR] Can't determine script dir: more than one '--file' argument detected\n")
+    return(script.dir)
+}
+## --------------------
+## Some html creation functions
+## --------------------
+htmlTableRow <- function( string_array=c() ){
+	td_cells <- ''
+	for ( i in string_array ){ 
+		td_cells <- paste( td_cells, '<td>', i, '</td>', sep='' )	
+	}
+	return( paste( "<tr>", td_cells, "</tr>") )
+}
+htmlLink <- function( path, desc="LINK" ){
+	return( paste( '<a href="', path, '">', desc, "</a>", sep='') )
+}
+## --------------------
+## constructs a list with input bam file info
+## --------------------
+makeBamFileList <- function( paths, names ){	
+	tmp <- list()
+	l1 <- length(paths)
+	l2 <- length(names)
+	if ( l1 != l2 ) stop( "Unequal amount of bam-paths (", l1, ") and -names (", l2, ") in makeBamFileList!!!\n" )
+	for ( i in 1:length(paths) ){
+		path <- paths[i]
+		name <- names[i]
+		file <- basename(path)
+
+		tmp[[ file ]] <- name
+		tmp[[ 'all_paths' ]] <- c( tmp[[ 'all_paths' ]], path )
+		tmp[[ 'all_files' ]] <- c( tmp[[ 'all_files' ]], file )
+		tmp[[ 'all_names' ]] <- c( tmp[[ 'all_names' ]], name )
+	}
+	return( tmp )
+}
+
+## --------------------
+## copied code for extracting the regions by segment call status
+## --------------------
+fuseRegions <- function( obj, minRatio=0 ) {
+	if ( ncol(obj) > 1 ) stop('Please specify which sample...')
+
+	data <- data.frame( obj@featureData@data[,1:3], copynumber(obj), segmented(obj), check.names=FALSE, stringsAsFactors=FALSE)
+	colnames( data ) <- c( "chr", "start", "end", "log2", "segmentval" )
+	
+	fused.data <- data.frame()
+	curr.bin <- 1
+	for ( chr in unique( data$chr ) ) {
+		chr.data  <- data[ data$chr == chr, ]
+		prev.bin  <- curr.bin
+		prev.log2 <- chr.data[ 1, 'log2' ]
+		prev.segm <- chr.data[ 1, 'segmentval' ]
+		start     <- chr.data[ 1, 'start' ]
+
+		if ( nrow(chr.data) > 1) {
+			for ( i in 2:nrow(chr.data) ) {
+				curr.bin  <- curr.bin + 1
+				curr.segm <- chr.data[ i, 'segmentval']
+
+				if ( curr.segm != prev.segm ) {
+					fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
+					prev.segm <- curr.segm
+					prev.bin  <- curr.bin
+					start     <- chr.data[ i, 'start']
+				}
+			}
+			fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
+		}else{
+			fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
+		}
+	}
+	## remove regions with low amplitude
+	fused.data <- fused.data[ abs(fused.data$segmentval) >= minRatio, ]
+	fused.data
+}
+
+## DESC: takes the output of fuse.regions and outputs a txt file per sample
+outputRegionsFromList <- function ( regionsList, outputBasename, outputDir="./", binSize, storeList ){
+	if ( missing(regionsList) ) stop( 'Please provide regionsList...' )
+	if ( missing(outputBasename) ) stop( 'Please provide outputBasename...' )
+	if ( !is.list(regionsList) ) stop( 'Input not a list...?' )
+	if ( length(regionsList) < 1 ) stop( 'List seems empty...?' )
+	if ( file.exists( outputDir ) ) catMsg( c(" Using dir ", outputDir, " for output") )
+	else dir.create( outputDir )
+
+	## have to set R output options otherwise scientific method is used at some point
+	options( "scipen"=100 )
+
+	sampleCount <- length( regionsList )
+	sampleNames <- names( regionsList )
+	bedgraphColumns <- c( 'chr', 'start', 'end', 'segmentval' )
+	
+	catMsg( c( " There are ", sampleCount, " samples found in input list") )
+
+	for ( sample in sampleNames ){		
+		catMsg( c(" Working on sample ", sample ) )
+		regionCount <- nrow( regionsList[[sample]] )
+		
+		outSampleBase   <- paste( outputBasename, '_', sample, '_', binSize, 'kbp', sep='')
+		outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
+		outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )
+
+		## ---------- BEDGRAPH ----------
+		txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample," description=segmented_regions_from_QDNAseq_",binSize,"kbp\n", sep="")
+		sink( outBedgraphPath )
+			cat( txt )
+		sink()
+		write.table( regionsList[[sample]][,bedgraphColumns], outBedgraphPath, quote=F, sep="\t", row.names=F, append=T, col.names=F)
+		#outFiles[[sample]] <- c( outBedgraphFile )
+		storeList[[ paste( binSize, sample, 'bedgraph', sep="_")]] <- outBedgraphFile
+	}
+	return(storeList)
+}
+
+
+## ==================================================
+## Start of analysis
+## ==================================================
+MAIN_NAME <- '[INFO] '
+TOOL_PATH <- getScriptPath()
+CSS_FILE  <- paste( TOOL_PATH, '/static/css/QDNAseq.css', sep="" )
+DECIMALS  <- 3
+WEB_LINK  <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
+PURE_CSS  <- 'http://yui.yahooapis.com/pure/0.5.0/pure-min.css'
+
+catMsg( "Starting QDNAseq wrapper" )	
+catMsg( "Loading R libraries" )
+suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
+suppressWarnings( suppressMessages( library( CGHcall, quietly = TRUE ) ) )
+
+## only one param: the tmp config file
+cmdLineArgs <- commandArgs(TRUE)
+config      <- cmdLineArgs[1]
+
+## sourcing the config file will load all input params
+## many variables are imported via sourced "config"
+source( config )
+
+## desparate tries to make png text scale well, damn you R...!
+PLOT_RES  <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 6.3 
+#PLOT_PS <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 40
+#TXT_FACTOR <- PLOT_WIDTH/480
+#PAR_SET   <- list( pch=22, cex.lab=TXT_FACTOR, cex.main=TXT_FACTOR )
+PAR_SET   <- list( pch=22 )
+
+systemUser <- system("whoami",T)
+qdnaseqVersion <- packageDescription( "QDNAseq" )$Version
+rVersion <- R.version.string
+startTime <- Sys.time()
+analysisStart <- as.character( startTime )
+catMsg( c("QDNAseq version: ", qdnaseqVersion) )
+catMsg( c( rVersion ) )
+
+## get the comma separated list of chromosomes to exclude
+excludeChrs <- unlist( strsplit( excludeChrsString, ",") )
+binSizes <- as.numeric( unlist( strsplit( binSizesString, ",") ) )
+
+
+## ------------------------
+## DEBUG
+if ( debug ){
+	catMsg( c("Analysis run under user: ", systemUser ) )
+	catMsg( c("Here comes sessionInfo: " ) )
+	print( sessionInfo() )
+}
+## /DEBUG
+## ------------------------
+
+## prepare output dir
+if ( !file.exists( outputPath) ){
+	dir.create( outputPath )
+}
+
+## copy source config file to output dir to include it in output zip
+file.copy( config, paste(outputPath, 'galaxyConfigFile.R', sep='/') )
+
+## setup bam filelist for easy retrieval later
+fileList    <- makeBamFileList( bamsPaths, bamsNames )
+bamCount    <- length( fileList[[ 'all_paths' ]] )
+
+gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
+gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep='')
+htmlOutputName <- 'index.html'
+htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep='')
+
+plotted_images <- list() # to keep track of images for later linking
+regions <- list() # will contain the segments
+outputFiles <- list()
+
+## ------------------------
+## in case of debug just use inbuilt LGG data for speedup
+if ( debug ){ 
+	binSizes <- c(15)
+	bamsPaths  <- c( "BUILD_IN_DATA")
+	bamsNames  <- c( "LGG150")
+	fileList   <- makeBamFileList( bamsPaths, bamsNames )
+	bamCount   <- length( fileList[[ 'all_paths' ]] )
+}
+
+for ( binSize in binSizes ){
+
+	## ------------------------
+	## construct output file-names and -paths
+	## ------------------------
+	robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
+	robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
+	igvCopyNrName  <- paste( binSize, 'kbp_QDNAseqCopyNumbers.igv', sep='')
+	robjSegmntName <- paste( binSize, 'kbp_QDNAseqCopyNumbersSegmented.rds', sep='')
+	regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
+	noiseImgName   <- paste( binSize, 'kbp_QDNAseqNoiseplot.png', sep='')
+	
+	robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep='')
+	robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep='')
+	robjSegmntPath <- paste( outputPath, '/', robjSegmntName, sep='')
+	robjRegionPath <- paste( outputPath, '/', regiOutputName, sep='')
+	igvCopyNrPath  <- paste( outputPath, '/', igvCopyNrName, sep='')
+	noiseImgPath   <- paste( outputPath, '/', noiseImgName, sep='')
+
+	binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
+
+	## in case of debug just use inbuilt LGG data for speedup
+	if ( debug ){
+		data(LGG150)
+		readCounts <- LGG150
+	}else{
+		## provide bamnames because in galaxy everyting is called "dataset_###"
+		readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=fileList[[ 'all_names' ]] )	
+	}
+
+	readCountsFiltered    <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
+	readCountsFiltered    <- estimateCorrection( readCountsFiltered )
+	copyNumbers           <- correctBins( readCountsFiltered )
+	copyNumbersNormalized <- normalizeBins( copyNumbers )
+	copyNumbersSmooth     <- smoothOutlierBins( copyNumbersNormalized )
+	sampleNames           <- readCountsFiltered@phenoData@data$name
+
+	## save objects to output dir
+	saveRDS( readCountsFiltered, robjReadCoPath );
+	saveRDS( copyNumbersSmooth, robjCopyNrPath );
+	exportBins( copyNumbersSmooth, file=igvCopyNrPath, format="igv" )
+
+	## also save objects for galaxy history output if requested
+	if ( doOutputCopynumbersIgv ){
+		#@ a bit hacky galaxy way to allow an unknown number of output files based on param selection
+		#@ see: https://wiki.galaxyproject.org/Admin/Tools/Multiple%20Output%20Files
+		historyName <- paste(binSize, 'kbp-IGV', sep="")
+		igvFile <- paste( newFilePath, "/primary_", outputId, "_", historyName, "_visible_txt", sep="" )
+		exportBins( copyNumbersSmooth, file=igvFile, format="igv" )
+		catMsg( c("Exported igv file to history for ", binSize, "kbp bin") )
+	}
+
+	## proceed with calling if requested
+	if ( doSegment ){
+		copyNumbersSegmented  <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
+		copyNumbersSegmented  <- normalizeSegmentedBins( copyNumbersSegmented )
+		cgh <- makeCgh( copyNumbersSegmented )
+		saveRDS( copyNumbersSegmented, robjSegmntPath );
+	}
+
+	## ------------------------
+	## create output files
+	## ------------------------
+	png( noiseImgPath, width=PLOT_HEIGHT, height=PLOT_HEIGHT, res=PLOT_RES );
+		par( PAR_SET )
+		noisePlot( readCountsFiltered, main=paste( "Noise Plot ", binSize, "kbp", sep=''), col="darkgreen" )
+	dev.off()
+
+	binSize <- as.character( binSize ) # to avoid R using it as array index... (*#$^@ you R!)
+	binSizeString <- paste( binSize, 'kbp', sep='')
+
+	for (i in 1:length(sampleNames) ){
+		
+		sample <- sampleNames[i]
+		usedReads  <- readCountsFiltered@phenoData@data$used.reads[i]
+		catMsg( c("Creating plots for sample: ", sample, " (", binSizeString, ")" ) )	
+
+		type <- 'CopyNumbers'
+		img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png',  sep='')
+		img_file_path <- paste( outputPath, '/', img_file, sep='' )
+
+		## COPYNUMBER PLOT
+		png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES ); 
+			par( PAR_SET )
+			plot( copyNumbersSmooth[ ,sample ], main=paste(sample, ": CopyNumbers", sep="") ) 
+			mtext( paste( binSizeString, " bins", sep=""), 3 )
+			abline( h=c(-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
+		dev.off()
+		
+		plotted_images[[ paste(binSize, sample, type, sep="_" ) ]] <- img_file
+		
+		if ( doSegment ){
+			type <- 'Segmented'
+			img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png',  sep='')
+			img_file_path <- paste( outputPath, '/', img_file, sep='' )
+			
+			## COPYNUMBER PLOT
+			png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES ); 
+				par( PAR_SET )
+				plot( copyNumbersSegmented[ ,sample ], main=paste(sample, ": CopyNumbers and Segments", sep="") )
+				mtext( paste( "(", binSizeString, " bins)", sep=""), 3 )
+				abline( h=c(-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
+			dev.off()
+
+			plotted_images[[ paste(binSize, sample, type, sep="_" ) ]] <- img_file
+
+			catMsg( c(" Fusing regions of sample: ", sample) )
+			regions[[ sample ]] <- fuseRegions( cgh[, sample] )
+
+			region_count <- nrow( data.frame( regions[[ sample ]] ) )
+			catMsg( c( ' sample "', sample, '" has ', region_count, " regions" ) )
+			plotted_images[[ paste(binSize, sample, 'region_count', sep="_" ) ]] <- region_count
+		}
+		
+		## add USED read counts
+		plotted_images[[ paste(binSize, sample, 'usedReads', sep="_" ) ]] <- usedReads
+	}
+
+	if ( doSegment ){
+		saveRDS( regions, robjRegionPath )
+		plotted_images <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath, binSize=binSize, storeList=plotted_images )	
+	}
+}# end bin
+
+
+## ----- debug -----
+#catMsg( "done" )
+#q(status=0)
+## ---- /debug -----
+
+
+## ------------------------
+## prepare output
+## ------------------------
+catMsg( "...zipping output")
+zip_cmd <- paste( "zip -j", gzipOutputPath, paste(outputPath,'/*',sep='') ) ## -j is for removing dirs from the tree
+system( zip_cmd )
+
+## ------------------------
+## get filesizes for report
+## ------------------------
+zippedSize <- paste( round( file.info( gzipOutputPath )[["size"]] / 1e+6, digits=2 ), 'MB' )
+endTime <- Sys.time()
+timeDiff <- format( round( endTime - startTime, 3 ) )
+analysisEnd <- as.character( endTime )
+
+## ------------------------
+## creating html output to be linked to from the middle galaxy pane
+## ------------------------
+sink( file = htmlOutputPath, type = "output" )
+		cat( "<html>\n")
+		cat( "<head>\n")
+
+			cat( "\t", '<title>QDNAseq Report | ', outputName,'</title>', "\n", sep='' )
+			cat( "\t", '<link rel="stylesheet" href="', PURE_CSS, '">', "\n", sep='' )
+			cat( "\t<style>\n", sep='')
+				## include CSS into html file, makes it more portable
+				cat( "\t\t", readLines( CSS_FILE ), sep="\n\t\t" )
+				#cat( "\t\th1 {color:red;}", "\n")
+			cat( "\n\t</style>\n" )
+			
+		cat( "\n</head>\n")
+		cat( "\n<body>\n")
+
+		cat( "<h1>QDNAseq Report</h1>", "\n")
+		
+		cat( '<h3 class="qdnaseq">About this analysis</h3>', "\n")
+		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='')		
+		
+		## ------------------------
+		## table with general info
+		## ------------------------
+		cat( '<h3 class="qdnaseq">Settings</h3><p>', "\n")
+		cat( '<table class="pure-table pure-table-striped"><tbody>' )
+			cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
+			cat( htmlTableRow( c( "AnalysisStart", analysisStart ) ) )
+			cat( htmlTableRow( c( "AnalysisEnd", analysisEnd ) ) )
+			cat( htmlTableRow( c( "AnalysisTime", timeDiff ) ) )
+			cat( htmlTableRow( c( "BinSizes (kbp)", paste(binSizes,collapse=", ") ) ) )
+			cat( htmlTableRow( c( "R info", rVersion ) ) )
+			cat( htmlTableRow( c( "QDNAseq info", qdnaseqVersion ) ) )
+			
+			sampleStrings <- c()
+			for ( galaxyName in fileList[[ 'all_files' ]] ){
+				sampleName <- fileList[[ galaxyName ]]
+				sampleStrings <- c( sampleStrings, paste( galaxyName, ' (', sampleName, ')', sep='' ) )
+			}
+			cat( htmlTableRow( c( "InputBams", paste( sampleStrings, collapse=", ") ) ) )
+
+		cat( "</tbody></table></p>", "\n")
+		
+		## ------------------------
+		## list with links to all output files
+		## ------------------------
+		cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
+		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='')
+		cat( '<table class="pure-table pure-table-striped">', "\n" )
+		cat( '<thead><th>Type</th>', as.vector( mapply( paste, "<th>", binSizes, "kbp</th>", sep="" ) ),'</thead>', "\n" )
+		cat( "<tbody>", "\n")
+			files <- list()
+			fileTypes <- c( 'ReadCounts.rds', 'CopyNumbers.rds' )
+			if ( doSegment ){ fileTypes <- c( fileTypes, 'CopyNumbersSegmented.rds') }
+
+			for ( fileType in fileTypes ){
+				fileNames <- mapply( paste, binSizes, paste( 'kbp_QDNAseq', fileType, sep=''), sep='')
+				fileLinks <- mapply( htmlLink, fileNames, paste( binSizes, "kbp", sep="" ) )
+				cat( htmlTableRow( c( fileType, fileLinks ) ) )	
+			}
+		cat( "\n</tbody></table></p>", "\n")
+
+		## ------------------------
+		## table with links to files	
+		## ------------------------
+		ratio <- PLOT_WIDTH / PLOT_HEIGHT
+		width <- 960; height <- width / ratio ## bigger img
+		width_t <- 100; height_t <- 40 ## thumb img
+
+		cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
+		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='')
+		plots_html <- ''
+
+		colspan <- 1
+		binHeader <- "<th>Image</th>"
+		if ( doSegment ){ # extra column with segment info
+			colspan <- 2 
+			binHeader <- "<th>Image</th><th>Segments</th>"
+		} 
+		cat( '<table class="pure-table pure-table-striped">', "\n" )
+		cat( '<thead><tr><th></th><th></th>', as.vector( mapply( paste, "<th colspan=\"", colspan,"\">", binSizes, "kbp</th>", sep="" ) ), '</tr></thead>' )
+		cat( '<thead><tr><th>Sample / File</th><th>Reads</th>', rep( binHeader, length(binSizes) ), '</tr></thead>' )
+		cat( '<tbody>' )
+
+			for ( bam_file in bamsNames ){
+				
+				usedReads <- plotted_images[[ paste(binSize, bam_file, 'usedReads', sep="_" ) ]]
+				usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )
+				htmlRow <- paste( '<tr><td>', bam_file, '</td><td>', usedReads, '</td>', sep='' )
+
+				for ( binSize in binSizes ){
+					
+					## add thumbnails to table with links to anchors on html page
+					copy_img <- plotted_images[[ paste(binSize, bam_file, 'CopyNumbers', sep="_" ) ]]
+					html_copy_thumb <- htmlLink( path=paste('#', copy_img, sep=''), paste('<img src="',copy_img,'" alt="', bam_file, '" width="', width_t, '" height="', height_t, '">', sep='') )
+					html_copy_img <- htmlLink( path=copy_img, paste('<img id="', copy_img,'" src="',copy_img,'" alt="',bam_file, '" width="', width, '" height="', height, '">', sep='') )
+					html_segm_img <- ''
+					html_bedGraph <- ''
+					region_count <- ''
+					htmlRow <- paste( htmlRow, '<td>', html_copy_thumb, '</td>' )
+
+					if ( doSegment ){
+						segm_img <- plotted_images[[ paste(binSize, bam_file, 'Segmented', sep="_" ) ]]
+						region_count <- plotted_images[[ paste(binSize, bam_file, 'region_count', sep="_" ) ]]
+
+						html_bedGraph <- htmlLink( path=plotted_images[[ paste(binSize, bam_file, 'bedgraph', sep="_" ) ]], 'bedGraph' )
+						html_segm_img <- htmlLink( path=segm_img, paste('<img id="', segm_img,'" src="', segm_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
+						htmlRow <- paste( htmlRow, '<td>', region_count, ' (', html_bedGraph, ')</td>', sep="" )
+					}
+					plots_html <- paste( plots_html, html_copy_img, "\n", html_segm_img, "\n<br \\>\n", sep='' )	
+				}
+				plots_html <- paste( plots_html, "\n<hr \\>\n", sep='' )
+				## add info to overview table, including small thumbnails
+				htmlRow <- paste( htmlRow, '</tr>', sep='' )
+				cat( htmlRow, "\n" )				
+			}
+		cat( "</tbody></table></p>", "\n")
+		
+		## ------------------------
+		## section with various output shown
+		## ------------------------
+		cat( '<h3 class="qdnaseq">Results: Sample plots</h3><p>', "\n")
+		## now include (large) images in html page
+		cat( plots_html, "\n")
+		cat( "\n</p></body>\n")
+		cat( "\n</html>\n")
+sink()
+
+## ------------------------
+## creating main html output for galaxy history
+## ------------------------
+sink( file = outputHtml, type = "output" )
+		
+		cat( "<head>", "\n")
+			cat( "\t", '<link rel="stylesheet" href="', PURE_CSS, '">', "\n", sep='' )
+
+			cat( "<style>", "\n")
+				## include CSS directly into html file
+				cat( paste( "\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
+				cat( paste( "\t", readLines( CSS_FILE, n = -1)), sep="\n" )
+			cat( "</style>", "\n")
+		cat( "</head>", "\n")
+
+		cat( '<h1>QDNAseq Results (', outputName,')</h1>', "\n", sep="")
+		cat( '<p>Explore <a href="', htmlOutputName, '" class="button">the results</a> directly within galaxy</p>', "\n", sep="")
+		cat( '<p>Or download a <a href="', gzipOutputName, '" class="button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )
+
+sink()
+
+## ------------------------
+## create final zip and quit with status 0 to tell galaxy all was fine
+## ------------------------
+catMsg( "zipping all output")
+system( paste( "zip -j ", gzipOutputPath, paste(outputPath,'/', htmlOutputName, sep='') ) )
+catMsg( "done" )
+q(status=0)