view QDNAseq.R @ 55:b85c575ac1e3 draft

Uploaded

#!/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
## --------------------
fuse.regions_test <- function(cgh, onlyCalled=T) {
	if (ncol(cgh) > 1) stop('Please specify which sample...')

	x <- data.frame(cgh@featureData@data[,1:3], calls(cgh), copynumber(cgh), segmented(cgh), check.names=FALSE, stringsAsFactors=FALSE)
	colnames( x ) <- c( "chr", "start", "end", "call", "log2", "segmentval" )
	
	fused.data <- data.frame()
	curr.bin <- 1
	for (chr in unique(x$chr)) {
		
		chr.data <- x[x$chr == chr,]
		prev.bin <- curr.bin
		prev.call <- chr.data[1, 'call']
		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.call <- chr.data[ i, 'call']
				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'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
					if ( prev.call != 0 ){
						catMsg( c(" ...found called/segmented region (", chr, ':', start, ' call=', prev.call, ' segment=', prev.segm, ")" ) )
					}
					prev.call <- curr.call
					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'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
		}else {
			fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
		}
	}
	if ( onlyCalled == T ){
		fused.data <- fused.data[ fused.data$call != 0, ]
	}
	fused.data
}

## DESC: takes the output of fuse.regions and outputs a txt file per sample
outputRegionsFromList <- function ( regionsList, outputBasename, outputDir="./" ){
	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 )
	outFiles <- list()

	## 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, '_QDNAseqRegions', sep='')
		outBedFile      <- paste( outSampleBase, '.bed', sep="" )
		outBedPath      <- paste( outputDir, '/', outBedFile, sep="" )
		outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
		outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )

		## ---------- BED ----------
		txt <- "#"
		sink( outBedPath )
			cat( txt )
		sink()
		write.table( regionsList[[sample]], outBedPath, quote=F, sep="\t", row.names=F, append=T)

		## ---------- BEDGRAPH ----------
		txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample,"_segmReg description=segmented_regions_from_QDNAseq\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( outBedFile, outBedgraphFile )
	}
	outFiles
}

## ==================================================
## Start of analysis
## ==================================================
MAIN_NAME <- '[INFO] '
TOOL_PATH <- getScriptPath()
CSS_FILE  <- paste( TOOL_PATH, '/QDNAseq.css', sep="" )
DECIMALS  <- 3
WEB_LINK  <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
PAR_SET   <- list( pch=22 )

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 created in sourced "config"
source( config )

PLOT_RES  <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 6.3 # desparate try to make png text scale well, damn you R...!
if ( doOutputCallsRds == TRUE ){ doCall <- TRUE }

systemUser <- system("whoami",T)
qdnaseqVersion <- packageDescription( "QDNAseq" )$Version
catMsg( c("Analysis run with user: ", systemUser ) )
catMsg( c("QDNAseq version loaded: ", qdnaseqVersion) )
catMsg( c( sessionInfo()$R.version$version.string ) )

## get the comma separated list of chromosomes to exclude
excludeChrs <- unlist( strsplit( excludeChrsString, ",") )

## ------------------------
## DEBUG
#catMsg( "PARAM" )
#catMsg( galaxy_path )
#catMsg( repos_path )
#catMsg( instal_path )
## /DEBUG
## ------------------------

## setup bam filelist for easy retrieval later
fileList    <- makeBamFileList( bamsPaths, bamsNames )
bamCount    <- length( fileList[[ 'all_paths' ]] )

## help msg for running script separate still needs work!
#if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){ 
#	catMsg( c( "Usage: ", paste(params_help,sep=",") ) )
#	quit( save = 'no', status=0, runLast=F )
#}

if ( !file.exists( outputPath) ){
	dir.create( outputPath )
}

## because we circumvent params that galaxy can save, we want to
## copy source config file to output dir to include it in output zip
file.copy( config, paste(outputPath, 'qdnaseq_config_file.R', sep='/') )

## ------------------------
## construct output file-names and -paths
## ------------------------
htmlOutputName <- 'index.html'
gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
robjCalledName <- paste( binSize, 'kbp_QDNAseqCalledSegments.rds', sep='')
regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
igvCalledName  <- paste( binSize, 'kbp_QDNAseq-calls.igv', sep='')
igvCopyNrName  <- paste( binSize, 'kbp_QDNAseq-normalized.igv', sep='')
noiseImgName   <- paste( 'QDNAseqNoisePlot.png',  sep='')
freqImgName    <- paste( 'QDNAseqFrequencyPlot.png',  sep='')


gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep='')
htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep='')
robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep='')
robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep='')
robjCalledPath <- paste( outputPath, '/', robjCalledName, sep='')
robjRegionPath <- paste( outputPath, '/', regiOutputName, sep='')
igvCalledPath  <- paste( outputPath, '/', igvCalledName, sep='')
igvCopyNrPath  <- paste( outputPath, '/', igvCopyNrName, sep='')
noiseImgPath   <- paste( outputPath, '/', noiseImgName, sep='' )
freqImgPath    <- paste( outputPath, '/', freqImgName, sep='' )


## ------------------------
## performing QDNAseq analysis steps
## ------------------------
if ( debug ){
	## in case of debug just use inbuilt LGG data for speedup
	data(LGG150)
	readCounts <- LGG150
}else{
	if ( nchar(binAnnotations) == 0 ){
	binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
	}else{
		## if user provided file, check if correct class
		if ( class(binAnnotations)[1] != 'AnnotatedDataFrame' ){
			stop( "Provided binAnnotations file is not of class 'AnnotatedDataFrame'\n" )
		}
		binAnnotations <- readRDS( binAnnotations )
	}
	## provide bamnames because in galaxy everyting is called "dataset_###"
	readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=bamsNames )
}

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 ( doOutputReadcountsRds ){
	saveRDS( readCountsFiltered, readCountsDatasetFile );	
}
if ( doOutputCopynumbersRds ){
	saveRDS( copyNumbersSmooth, copyNumbersDatasetFile );
}

## proceed with calling if requested
if ( doCall ){
	copyNumbersSegmented  <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
	copyNumbersSegmented  <- normalizeSegmentedBins( copyNumbersSegmented )
	copyNumbersCalled     <- callBins( copyNumbersSegmented )
	cgh <- makeCgh( copyNumbersCalled )
	saveRDS( copyNumbersCalled, robjCalledPath );
	if ( doOutputCallsRds ){
		saveRDS( copyNumbersCalled, calledSegmentsDatasetFile );
	}
	exportBins( copyNumbersCalled, file=igvCalledPath, format="igv")
	
}



## ------------------------
## create output files
## ------------------------
plotted_images <- list() # to keep track of images for later linking
regions <- list() # will contain the (called) segments

png( noiseImgPath, width=PLOT_HEIGHT, height=PLOT_HEIGHT, res=PLOT_RES );
	par( PAR_SET )
	noisePlot( readCountsFiltered, col="darkgreen" )
dev.off()

if ( doCall ){
	png( freqImgPath, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES ); 
		par( PAR_SET )
		frequencyPlot( copyNumbersCalled )
	dev.off()	
}


for (i in 1:length(sampleNames) ){
#for (sample in sampleNames(copyNumbersSmooth) ){
	sample <- sampleNames[i]
	usedReads  <- readCountsFiltered@phenoData@data$used.reads[i]
	catMsg( c("Creating plots for sample: ", sample ) )	

	type <- 'CopyNumbers'
	img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png',  sep='')
	img_file_path <- paste( outputPath, '/', img_file, sep='' )
	png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES ); 
		par( PAR_SET )
		plot( copyNumbersSmooth[ ,sample ] ); 
	dev.off()
	plotted_images[[ sample ]][[ type ]] <- img_file

	if ( doCall ){
		type <- 'Called'
		img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png',  sep='')
		img_file_path <- paste( outputPath, '/', img_file, sep='' )
		png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES ); 
			par( PAR_SET )
			plot( copyNumbersCalled[ ,sample ] ); 
		dev.off()
		plotted_images[[ sample ]][[ type ]] <- img_file

		catMsg( c(" Fusing regions of sample: ", sample) )
		regions[[ sample ]] <- fuse.regions_test( cgh[, sample] )
		region_count <- nrow( data.frame( regions[[ sample ]] ) )
		catMsg( c( ' sample "', sample, '" has ', region_count, " regions" ) )
		plotted_images[[ sample ]][[ 'region_count' ]] <- region_count
	}

	## add USED read counts
	plotted_images[[ sample ]][[ 'usedReads'  ]] <- usedReads
}

if ( doCall ){
	saveRDS( regions, robjRegionPath )
	printedFiles <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath )	
}

## ------------------------
## 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"]] / 1000000, digits=2 ), 'MB' )
readCoSize <- paste( round( file.info( robjReadCoPath )[["size"]] / 1000000, digits=2 ), 'MB' )
copyNrSize <- paste( round( file.info( robjCopyNrPath )[["size"]] / 1000000, digits=2 ), 'MB' )
calledSize <- paste( round( file.info( robjCalledPath )[["size"]] / 1000000, digits=2 ), 'MB' )
regionSize <- paste( round( file.info( robjRegionPath )[["size"]] / 1000000, digits=2 ), 'MB' )
igvCopyNrSize <- paste( round( file.info( igvCopyNrPath )[["size"]] / 1000000, digits=2 ), 'MB' )
igvCalledSize <- paste( round( file.info( igvCalledPath )[["size"]] / 1000000, digits=2 ), 'MB' )

## ------------------------
## creating html output to be linked to from the middle galaxy pane
## ------------------------

#sink( file = outputHtml, type = "output" )
sink( file = htmlOutputPath, type = "output" )
		cat( "<html>\n")
		cat( "<head>\n")
			#cat( '<link rel="stylesheet" type="text/css" href="test.css" media="screen" />', "\n" )
			cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
			cat( "\t<style>\n")
				## have to include CSS into html file, because css referencing outside own dir doesn't seem to work...
				cat( paste( "\t\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
				cat( paste( "\t\t", readLines( CSS_FILE, n = -1)), sep="\n" )
				#cat( "\t\th1 {color:red;}", "\n")

			cat( "\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-success button-small pure-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"><thead><tr><th>setting</th><th>value</th></tr></thead><tbody>' )
			cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
			cat( htmlTableRow( c( "AnalysisDate", 'todo' ) ) )
			cat( htmlTableRow( c( "BinSize (kb)", binSize ) ) )
			
			for ( galaxyName in fileList[[ 'all_files' ]] ){
				sampleName <- fileList[[ galaxyName ]]
				cat( htmlTableRow( c( "InputBam", paste( galaxyName, ' (', sampleName, ')', sep='' ) ) ) )
			}
		cat( "</tbody></table></p>", "\n")
		
		## ------------------------
		## list with links to all output files
		## ------------------------
		r_code <- '<p>'
		r_code <- paste( r_code, '<code class="code">## R code to load files</code><br />', sep="\n" )
		r_code <- paste( r_code, '<code class="code">library( QDNAseq )</code><br />', sep="\n")

		cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
		cat( '<dl>', "\n" )
            cat( '<dt>', htmlLink( path=robjReadCoName, robjReadCoName ), '</dt>', "\n" ) 
            cat( '<dd>QDNAseq object with read counts per bin, ', readCoSize,'</dd>', "\n" )
            r_code <- paste( r_code, '<code class="code">readCounts <- readRDS(', robjReadCoName, ")</code><br />", sep="")

            cat( '<dt>', htmlLink( path=robjCopyNrName, robjCopyNrName ), '</dt>', "\n" ) 
            cat( '<dd>QDNAseq object with copy numbers per bin, ', copyNrSize,'</dd>', "\n" )
            r_code <- paste( r_code, '<code class="code">copyNumbersSmooth <- readRDS(', robjCopyNrName, ")</code><br />", sep="")
            
            cat( '<dt>', htmlLink( path=igvCopyNrName, igvCopyNrName ), '</dt>', "\n" ) 
            cat( '<dd>IGV formatted text file with copy numbers per bin, ', igvCopyNrSize,'</dd>', "\n" )
            
            if ( doCall ){
            	cat( '<dt>', htmlLink( path=robjCalledName, robjCalledName ), '</dt>', "\n" ) 
            	cat( '<dd>QDNAseq object with segment and call status per bin, ', calledSize,'</dd>', "\n" )
            	r_code <- paste( r_code, '<code class="code">copyNumbersCalled <- readRDS(', robjCalledName, ")</code><br />", sep="")

            	cat( '<dt>', htmlLink( path=regiOutputName, regiOutputName ), '</dt>', "\n" ) 
            	cat( '<dd>list with segmented/called regions for each sample, ', regionSize, '</dd>', "\n" )
            	r_code <- paste( r_code, '<code class="code">calledRegions <- readRDS(', regiOutputName, ")</code><br />", sep="")

            	cat( '<dt>', htmlLink( path=igvCalledName, igvCalledName ), '</dt>', "\n" ) 
            	cat( '<dd>IGV formatted text file with calls per bin , ', igvCalledSize,'</dd>', "\n" )

            }
        cat( '</dl></p>', "\n" )

        ## r_code shows example code on how to load output files
        #cat( r_code, "</p>\n", sep="\n")
        cat( '<p>See ', htmlLink( WEB_LINK, 'the bioconductor QDNAseq documentation' ), ' for more information on how to work with these files</p>', "\n", sep='' )

		## ------------------------
		## table with links to files	
		## ------------------------
		cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
		plots_html <- ''
		
		cat( '<table class="pure-table pure-table-striped"><thead><tr><th>Sample / File</th><th>CopyNr</th><th>Calls</th><th>Reads</th><th>Regions</th><th>Files</th></tr></thead><tbody>' )
			
			for ( bam_file in sampleNames ){
				
				#width <- 600; height <- 240
				width <- PLOT_WIDTH; height <- PLOT_HEIGHT
				width_t <- 100; height_t <- 40

				## add thumbnails to table with links to anchors on html page
				copy_img <- plotted_images[[ bam_file ]][[ 'CopyNumbers' ]]
				usedReads <- plotted_images[[ bam_file ]][[ 'usedReads' ]]
				usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )

				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_call_thumb <- 'NA'
				html_call_img <- ''
				html_bed <- 'NA'
				html_bedGraph <- 'NA'
				region_count <- 'NA'

				if ( doCall ){
					call_img <- plotted_images[[ bam_file ]][[ 'Called' ]]
					region_count <- plotted_images[[ bam_file ]][[ 'region_count' ]]
					html_call_thumb <- htmlLink( path=paste('#', call_img, sep=''), paste('<img src="', call_img, '" alt="', bam_file, '" width="', width_t,'" height="', height_t,'">', sep='') )

					files <- printedFiles[[ bam_file ]]
					html_bed <- htmlLink( path=files[1], 'bed' )	
					html_bedGraph <- htmlLink( path=files[2], 'bedGraph' )	
					html_call_img <- htmlLink( path=call_img, paste('<img id="', call_img,'" src="', call_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
				}

				## add info to overview table, including small thumbnails
				cat( htmlTableRow( c(bam_file, html_copy_thumb, html_call_thumb, usedReads, region_count, paste( html_bed, html_bedGraph, sep=", ")) ), "\n" )
				## now include (large) images in html page
				plots_html <- paste( plots_html, html_copy_img, "\n", html_call_img, "\n<hr \\>\n", sep='' )
			}
		cat( "</tbody></table></p>", "\n")

		## ------------------------
		## create the noise and frequency plots
		## ------------------------
		html_noise_img <- htmlLink( 
			path=noiseImgName, 
			paste('<img id="', noiseImgName,'" src="', noiseImgName, '" width="', height/2, '" height="', height/2, '" alt="NoisePlot">', sep='') 
		)
		html_freq_img  <- htmlLink( 
			path=freqImgName, 
			paste('<img id="', freqImgName,'" src="', freqImgName, '" width="', width/2, '" height="', height/2, ' alt="FrequenceyPlot">', sep='') 
		)
		extra_plots_html <- paste( 
			html_noise_img, " \n", 
			html_freq_img, "\n", sep='' 
		)
		
		## ------------------------
		## section with various output shown
		## ------------------------
		cat( '<h3 class="qdnaseq">Results: Analysis plots</h3><p>', "\n")
		cat( extra_plots_html, "\n")
		cat( '<h3 class="qdnaseq">Results: Sample plots</h3><p>', "\n")
		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="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )

			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-success button-small pure-button">the results</a> directly within galaxy</p>', "\n", sep="")
		cat( '<p>Or download a <a href="', gzipOutputName, '" class="button-success button-small pure-button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )

		cat( '<p>The zip file contains all output files, including *.rds files allowing you to load the R copyNumber object(s) and perform further detailed analysis or create your own output for further processing. You can load the rds file with <code class="code">loadRDS(file.rds)</code></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)