|
78
|
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
|
|
|
10
|
|
|
11 ## ==================================================
|
|
|
12 ## Start of analysis
|
|
|
13 ## ==================================================
|
|
|
14 MAIN_NAME <- '[INFO] '
|
|
|
15 catMsg( "Starting QDNAseq-plot wrapper" )
|
|
|
16 catMsg( "Loading R libraries" )
|
|
|
17
|
|
|
18 ## supress msg to allow R to finish with non-error msg
|
|
|
19 suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
|
|
|
20
|
|
|
21 ## only one param: the tmp config file
|
|
|
22 cmdLineArgs <- commandArgs(TRUE)
|
|
|
23 config <- cmdLineArgs[1]
|
|
|
24
|
|
|
25 ## sourcing the config file will load all input params
|
|
|
26 ## many variables are imported via sourced "config"
|
|
|
27 source( config ) # outputPngPath, outputPdfPath, allOrOne, rdsFilePath
|
|
|
28 #cat( "ALL? ", allOrOne, sep='' )
|
|
|
29
|
|
|
30 ## desparate tries to make png text scale well, damn you R...!
|
|
|
31 PLOT_RES <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 6.3
|
|
|
32 PAR_SET <- list( pch=22 )
|
|
|
33 systemUser <- system("whoami",T)
|
|
|
34 qdnaseqVersion <- packageDescription( "QDNAseq" )$Version
|
|
|
35 rVersion <- R.version.string
|
|
|
36 catMsg( c("QDNAseq version: ", qdnaseqVersion) )
|
|
|
37 catMsg( c( rVersion ) )
|
|
|
38
|
|
|
39 qdnaseqObject <- readRDS( rdsFilePath )
|
|
81
|
40 chromosomesToPlot <- unlist( strsplit( chromosomesToPlotString, ",") )
|
|
|
41
|
|
|
42 #cat( "CHROM: ", chromosomesToPlotString, "\n" )
|
|
|
43 #cat( "REGION: ", regionToPlotString, "\n" )
|
|
|
44 #cat( "What to plot: ", whatToPlot, "\n" )
|
|
|
45
|
|
78
|
46 ## COPYNUMBER PLOT
|
|
|
47 sample <- SAMPLE_INDEX
|
|
|
48 png( outputPngPath, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES )
|
|
|
49 par( PAR_SET )
|
|
87
|
50 ylab_text <- "log2 read counts"
|
|
81
|
51 if ( whatToPlot == 'everything' ){
|
|
|
52 catMsg( c( "Plotting all data in object" ) )
|
|
87
|
53 plot( qdnaseqObject[ ,sample ], ylab=ylab_text )
|
|
88
|
54 abline( h=c(-4,-3,-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
|
|
81
|
55 } else if( whatToPlot == 'chromosomes' ){
|
|
|
56 catMsg( c( "Plotting subset of chromosomes" ) )
|
|
|
57 fdata <- qdnaseqObject@featureData@data
|
|
|
58 idx_region <- which( fdata$chromosome %in% chromosomesToPlot )
|
|
87
|
59 plot( qdnaseqObject[ idx_region, sample ], ylab=ylab_text )
|
|
88
|
60 abline( h=c(-4,-3,-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
|
|
81
|
61 } else if( whatToPlot == 'region' ){
|
|
|
62 regionC <- chrName
|
|
|
63 regionS <- chrStart
|
|
|
64 regionE <- chrEnd
|
|
|
65 if ( regionS > regionE ) stop("Chosen start is > end")
|
|
|
66 catMsg( c( "Plotting genomic region (chr=", regionC, " start=", regionS, " end=", regionE, ")" ) )
|
|
|
67
|
|
|
68 fdata <- qdnaseqObject@featureData@data
|
|
|
69 idx_region <- which( fdata$chromosome == regionC & fdata$start > regionS & fdata$end < regionE )
|
|
87
|
70 cat( idx_region, "\n")
|
|
81
|
71
|
|
87
|
72 plot( qdnaseqObject[ idx_region, sample ], doCalls=FALSE, ylab=ylab_text )
|
|
81
|
73 }
|
|
78
|
74 #mtext( "plotted in galaxy", 3 )
|
|
87
|
75
|
|
78
|
76 dev.off()
|
|
|
77
|
|
|
78
|
|
|
79 ## all ok
|
|
|
80 q(status=0)
|
