0
|
1 #!/usr/bin/Rscript
|
|
2
|
|
3 #-----------------------------------#
|
|
4 # Author: Maude #
|
|
5 # Script: estimateSign_Galaxy.r #
|
|
6 # Last update: 22/07/15 #
|
|
7 #-----------------------------------#
|
|
8
|
|
9 #########################################################################################################################################
|
|
10 # Estimate the number of signatures for NMF #
|
|
11 #########################################################################################################################################
|
|
12
|
|
13 #-------------------------------------------------------------------------------
|
|
14 # Load library for recovering the arguments
|
|
15 #-------------------------------------------------------------------------------
|
|
16 suppressMessages(suppressWarnings(require("getopt")))
|
|
17
|
|
18
|
|
19 #-------------------------------------------------------------------------------
|
|
20 # Recover the arguments
|
|
21 #-------------------------------------------------------------------------------
|
|
22 spec = matrix(c(
|
|
23 "input" , "i", 1, "character",
|
|
24 "stop", "stop", 1, "numeric",
|
|
25 "cpu", "cpu", 1, "integer",
|
|
26 "output", "o", 1, "character",
|
|
27 "help", "h", 0, "logical"
|
|
28 ),
|
|
29 byrow=TRUE, ncol=4
|
|
30 )
|
|
31
|
|
32 opt = getopt(spec);
|
|
33
|
|
34 # No argument is pass to the command line
|
|
35 if(length(opt) == 1)
|
|
36 {
|
|
37 cat(paste("Usage:\n estimateSign_Galaxy.r --input <matrix> --stop <maxNbSign> --cpu <cpu> --output <output_filename.png>\n",sep=""))
|
|
38 q(status=1)
|
|
39 }
|
|
40
|
|
41 # Help was asked for.
|
|
42 if ( !is.null(opt$help) )
|
|
43 {
|
|
44 # print a friendly message and exit with a non-zero error code
|
|
45 cat(paste("Usage:\n estimateSign_Galaxy.r --input <matrix> --stop <maxNbSign> --cpu <cpu> --output <output_filename.png>\n",sep=""))
|
|
46 q(status=1)
|
|
47 }
|
|
48
|
|
49
|
|
50
|
|
51 #-------------------------------------------------------------------------------
|
|
52 # Load library
|
|
53 #-------------------------------------------------------------------------------
|
|
54 suppressMessages(suppressWarnings(library(NMF)))
|
|
55
|
|
56
|
|
57 ###############################################################################
|
|
58 # Load the functions #
|
|
59 ###############################################################################
|
|
60
|
|
61 #-------------------------------------------------------------------------------
|
|
62 # Check the file doesn't have lines equal to zero
|
|
63 #-------------------------------------------------------------------------------
|
|
64 CheckFile <- function(rowsum, dataFrame, x)
|
|
65 {
|
|
66 if(rowsum == 0)
|
|
67 {
|
|
68 write("\n\nERROR: There is not enough mutations for running NMF!!!", stderr())
|
|
69 write(paste0("Input matrix contains at least one null row ", rownames(dataFrame)[x], "\n\n"), stderr())
|
|
70 stop()
|
|
71 }
|
|
72 }
|
|
73
|
|
74
|
|
75 ###############################################################################
|
|
76 # Check file #
|
|
77 ###############################################################################
|
|
78
|
|
79 # The input musn't contains lines equal to zero !!!
|
|
80 matrixNMF <- read.table(opt$input, header=T)
|
|
81 # suppresses the return of sapply function
|
|
82 invisible( sapply(1:nrow(matrixNMF), function(x) { CheckFile(rowSums(matrixNMF)[x], matrixNMF, x) } ) )
|
|
83
|
|
84
|
|
85
|
|
86 ###############################################################################
|
|
87 # Estimate the number of signatures #
|
|
88 ###############################################################################
|
|
89 # Estimate the number of signatures with our data
|
|
90 nbCPU <- paste0("vP", opt$cpu)
|
|
91 nbSign <- 2:opt$stop # The minum number of signatures can't be lower than 2
|
|
92
|
|
93 estim_r <- nmf(matrixNMF, method="brunet", nbSign, nrun=50, .opt=nbCPU)
|
|
94
|
|
95 # Shuffle original data
|
|
96 v_random <- randomize(matrixNMF)
|
|
97 # Estimate quality measures from the shuffled data
|
|
98 estim_r_random <- nmf(v_random, method="brunet", nbSign, nrun=50, .opt=nbCPU)
|
|
99
|
|
100 # Plot the estimation for our data and the random ones
|
|
101 graphics.off()
|
|
102 options(bitmapType='cairo')
|
|
103 png(opt$output, width=3000, height=2000, res=300)
|
|
104 plot(estim_r, estim_r_random)
|
|
105 invisible( dev.off() )
|