Mercurial > repos > stef > qdnaseq

changeset 28:40ae95ad9d8c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author stef
date Wed, 18 Jun 2014 05:29:59 -0400
parents 47c3ecce5544
children 4c4e4d88779c
files QDNAseq.R
diffstat 1 files changed, 35 insertions(+), 37 deletions(-) [+]
line wrap: on
line diff
--- a/QDNAseq.R	Wed Jun 18 05:29:50 2014 -0400
+++ b/QDNAseq.R	Wed Jun 18 05:29:59 2014 -0400
@@ -7,7 +7,7 @@
 	cat( MAIN_NAME, paste( msg, collapse="" ), "\n", sep='')
 }
 ## --------------------
-## return the directory this script exists
+## return the location of this script
 ## --------------------
 getScriptPath <- function(){
     cmd.args <- commandArgs()
@@ -18,7 +18,7 @@
     return(script.dir)
 }
 ## --------------------
-## Some html helper functions
+## Some html creation functions
 ## --------------------
 htmlTableRow <- function( string_array=c() ){
 	td_cells <- ''
@@ -31,7 +31,7 @@
 	return( paste( '<a href="', path, '">', desc, "</a>", sep='') )
 }
 ## --------------------
-## constructs a list with bam file info
+## constructs a list with input bam file info
 ## --------------------
 makeBamFileList <- function( paths, names ){	
 	tmp <- list()
@@ -52,7 +52,7 @@
 }
 
 ## --------------------
-## took a function from Matias/MarkVdWuel? for extracting the regions by call
+## 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...')
@@ -150,24 +150,24 @@
 	outFiles
 }
 
-printIgvFile <- function( dat, filename ){
-	
-	if ( 'calls' %in% assayDataElementNames(dat) ) {
-	 	#output <- paste(output, '-called.igv', sep="")
-	  	cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
-	  	df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), calls(dat), check.names=FALSE, stringsAsFactors=FALSE)
-	}else {
-	  #output <- paste(output, '-normalized.igv', sep="")
-	  cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
-	  df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), round(copynumber(dat), digits=2), check.names=FALSE, stringsAsFactors=FALSE)
-	}
-
-	df$chromosome[df$chromosome == '23'] <- 'X'
-	df$chromosome[df$chromosome == '24'] <- 'Y'
-	df$chromosome[df$chromosome == '25'] <- 'MT'
-	#return( df )
-	write.table( df, file=filename, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )	
-}
+#printIgvFile <- function( dat, filename ){
+#	
+#	if ( 'calls' %in% assayDataElementNames(dat) ) {
+#	 	#output <- paste(output, '-called.igv', sep="")
+#	  	cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
+#	  	df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), calls(dat), check.names=FALSE, stringsAsFactors=FALSE)
+#	}else {
+#	  #output <- paste(output, '-normalized.igv', sep="")
+#	  cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
+#	  df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), round(copynumber(dat), digits=2), check.names=FALSE, stringsAsFactors=FALSE)
+#	}
+#
+#	df$chromosome[df$chromosome == '23'] <- 'X'
+#	df$chromosome[df$chromosome == '24'] <- 'Y'
+#	df$chromosome[df$chromosome == '25'] <- 'MT'
+#	#return( df )
+#	write.table( df, file=filename, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )	
+#}
 
 
 ## ==================================================
@@ -211,18 +211,18 @@
 ## ------------------------
 
 ## setup bam filelist for easy retrieval later
-catMsg( "Setting up input bam list" )
-cat( bamsPaths, "\n" )
-catMsg( "Namews" )
-cat( bamsNames, "\n" )
+#catMsg( "Setting up input bam list" )
+#cat( bamsPaths, "\n" )
+#catMsg( "Namews" )
+#cat( bamsNames, "\n" )
 fileList    <- makeBamFileList( bamsPaths, bamsNames )
 bamCount    <- length( fileList[[ 'all_paths' ]] )
 
 ## help msg still needs work!
-if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){ 
-	cat( paste( MAIN_NAME, "Usage: ", params_help, sep=''), "\n" )
-	quit( save = 'no', status=0, runLast=F )
-}
+#if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){ 
+#	cat( paste( MAIN_NAME, "Usage: ", params_help, sep=''), "\n" )
+#	quit( save = 'no', status=0, runLast=F )
+#}
 
 if ( !file.exists( outputPath) ){
 	dir.create( outputPath )
@@ -280,12 +280,12 @@
 copyNumbers           <- correctBins( readCountsFiltered )
 copyNumbersNormalized <- normalizeBins( copyNumbers )
 copyNumbersSmooth     <- smoothOutlierBins( copyNumbersNormalized )
+sampleNames           <- readCountsFiltered@phenoData@data$name
 
 ## save objects to output dir
 saveRDS( readCounts, robjReadCoPath );
-saveRDS( copyNumbersSmooth, robjCopyNrPath );
-#printIgvFile( copyNumbersSmooth, igvCopyNrName )
-#exportBins(copyNumbersSmooth, file=igvCopyNrName, format="igv")
+saveRDS( copyNumbersSmooth, robjCopyNrPath );s
+exportBins(copyNumbersSmooth, file=igvCopyNrName, format="igv")
 
 ## also save objects for galaxy history output if requested
 if ( doOutputReadcountsRds ){
@@ -305,12 +305,10 @@
 	if ( doOutputCallsRds ){
 		saveRDS( copyNumbersCalled, calledSegmentsDatasetFile );
 	}
-	#df <- getIGVdf( copyNumbersCalled )
-	#printIgvFile( copyNumbersCalled, igvCalledName )
-	#write.table( df, file=igvCalledName, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )
+	exportBins( copyNumbersCalled, file=igvCalledName, format="igv")
 }
 
-sampleNames <- readCountsFiltered@phenoData@data$name
+
 
 ## ------------------------
 ## create output files