|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 # Compute metrics (all except PCA) for ROI files.
|
|
|
4 #
|
|
|
5 # USAGE: roi_metrics.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 some basic statistics.
|
|
|
56 roi$min = as.numeric(lapply(roi$counts, min))
|
|
|
57 roi$median = as.numeric(lapply(roi$counts, median))
|
|
|
58 roi$mean = as.numeric(lapply(roi$counts, mean))
|
|
|
59 roi$quantile75 = as.numeric(lapply(roi$counts, quantile, probs=0.75))
|
|
|
60 roi$quantile95 = as.numeric(lapply(roi$counts, quantile, probs=0.95))
|
|
|
61
|
|
|
62 # Compute area under curve metric.
|
|
|
63 aoc <- function(v){
|
|
|
64 m = max(v)
|
|
|
65 return (sum(v/m)/length(v))
|
|
|
66 }
|
|
|
67 roi$aoc = as.numeric(lapply(roi$counts, aoc))
|
|
|
68
|
|
|
69 # Compute xreaXX metric.
|
|
|
70 xreaXX <- function(v){
|
|
|
71 vo=order(v)
|
|
|
72 mx = max(v)
|
|
|
73 mi = min(v)
|
|
|
74 if((mx-mi)==0) return(0)
|
|
|
75 vn=(v-mi)/(mx-mi)
|
|
|
76 x=c(1:length(v))/length(v)
|
|
|
77 vno=vn[vo]
|
|
|
78 return(sum(x-vno)/length(v))
|
|
|
79 }
|
|
|
80 roi$xreaXX = as.numeric(lapply(roi$counts, xreaXX))
|
|
|
81
|
|
|
82 # Compute r3linear metric.
|
|
|
83 roi$r3linear = as.numeric(lapply(roi$counts, k3.linear))
|
|
|
84
|
|
|
85 # Compute distMax3p metric.
|
|
|
86 distMax3p <- function(v)
|
|
|
87 {
|
|
|
88 return(length(v)-which.max(rev(v))+1)
|
|
|
89 }
|
|
|
90 roi$distMax3p = as.numeric(lapply(roi$count, distMax3p))
|
|
|
91
|
|
|
92 # Compute distMax5p metric.
|
|
|
93 distMax5p <- function(v)
|
|
|
94 {
|
|
|
95 return (which.max(v))
|
|
|
96 }
|
|
|
97 roi$distMax5p = as.numeric(lapply(roi$count, distMax5p))
|
|
|
98
|
|
|
99 # Write output.
|
|
|
100 writeROI(roi, opt$outfile)
|
|
|
101
|
|
|
102 # Quit with success code.
|
|
|
103 q(status=0) |
