|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 # Compute PCA metric for ROI files.
|
|
|
4 #
|
|
|
5 # USAGE: roi_pca.Rscript -i IN.roi -o OUT.roi
|
|
|
6
|
|
|
7 # TODO(holtgrew): Bernd, can you give the English names of the metrics below in the comments?
|
|
|
8
|
|
|
9 # Author: Bernd Jagla <bernd.jagla@pasteur.fr>
|
|
|
10 # Author: Manuel Holtgrewe <manuel.holtgrewe@fu-berlin.de>
|
|
|
11
|
|
|
12 library('getopt') # parsing args
|
|
|
13 library('ngsroi') # ROI I/O
|
|
|
14 library('boot') # k3.linear()
|
|
|
15
|
|
|
16 # ----------------------------------------------------------------------------
|
|
|
17 # Command Line Parsing
|
|
|
18 # ----------------------------------------------------------------------------
|
|
|
19
|
|
|
20 # Parse command line options.
|
|
|
21 spec = matrix(c(
|
|
|
22 "help", "h", 0, "logical",
|
|
|
23 "infile", "i", 1, "character",
|
|
|
24 "outfile", "o", 1, "character"
|
|
|
25 ), byrow=TRUE, ncol=4);
|
|
|
26 opt = getopt(spec);
|
|
|
27
|
|
|
28 # If help was asked for print a friendly message and exit with a non-zero
|
|
|
29 # error code.
|
|
|
30 if (!is.null(opt$help) )
|
|
|
31 {
|
|
|
32 cat(getopt(spec, usage=TRUE));
|
|
|
33 q(status=1);
|
|
|
34 }
|
|
|
35
|
|
|
36 # Check arguments.
|
|
|
37 if (is.null(opt$infile))
|
|
|
38 {
|
|
|
39 write("-i/--infile missing.", stderr());
|
|
|
40 q(status=1)
|
|
|
41 }
|
|
|
42 if (is.null(opt$out))
|
|
|
43 {
|
|
|
44 write("-o/--outfile missing.", stderr());
|
|
|
45 q(status=1)
|
|
|
46 }
|
|
|
47
|
|
|
48 # ----------------------------------------------------------------------------
|
|
|
49 # Compute metrics.
|
|
|
50 # ----------------------------------------------------------------------------
|
|
|
51
|
|
|
52 # Read input files.x
|
|
|
53 roi = readROI(opt$infile)
|
|
|
54
|
|
|
55 # Compute PCA.
|
|
|
56 data=roi[,8:length(roi)-1]
|
|
|
57 pca=princomp(data, na.action=na.exclude)
|
|
|
58 prj=predict(pca, data)
|
|
|
59 roi$pca1=prj[,1]
|
|
|
60 roi$pca2=prj[,2]
|
|
|
61
|
|
|
62 # Write output.
|
|
|
63 writeROI(roi, opt$outfile)
|
|
|
64
|
|
|
65 # Quit with success code.
|
|
|
66 q(status=0) |
