|
2
|
1 #!/usr/bin/env Rscript
|
|
|
2 ## IPL selection script by Peter Waltman
|
|
|
3 ## August 21, 2011
|
|
|
4 ## License under Creative Commons Attribution 3.0 Unported (CC BY 3.0)
|
|
|
5 ##
|
|
|
6 #usage, options and doc goes here
|
|
|
7 argspec <- c("ipl.feature.selection.R takes a set of results from Paradigm, and filters for features that are
|
|
|
8 active, inactive or modulated above a given IPL threshold over a sufficient percentage of samples.
|
|
|
9
|
|
|
10 Usage:
|
|
|
11 ipl.feature.selection.R -d <data.file>
|
|
|
12 Optional:
|
|
|
13 -o <output.name>
|
|
|
14 -g <genes-only> ## to set if only returning genes (default is all features)
|
|
|
15 -f <filter.type> ## filter.type must be either 'modulated', 'active'or 'inactive' (default is modulated)
|
|
|
16 -t <threshold> ## the threshold to use for the filter (default is 0.25)
|
|
|
17 -p <perc.pass> ## the percentage of samples that must pass the filter (default is 0.33)
|
|
|
18 -v <verbose> ## to set verbose on
|
|
|
19
|
|
|
20 \n\n")
|
|
|
21 args <- commandArgs(TRUE)
|
|
|
22 if ( length( args ) == 1 && args =="--help") {
|
|
|
23 write(argspec, stderr())
|
|
|
24 q();
|
|
|
25 }
|
|
|
26
|
|
|
27 lib.load.quiet <- function( package ) {
|
|
|
28 package <- as.character(substitute(package))
|
|
|
29 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
30 }
|
|
|
31 lib.load.quiet(getopt)
|
|
|
32
|
|
|
33 spec <- matrix( c( "data.fname", "d", 1, "character",
|
|
|
34 "output.name", "o", 2, "character",
|
|
|
35 "genes.only", "g", 0, "logical",
|
|
|
36 "filter.type", "f", 2, "character", ## must be either 'active', 'inactive' or 'modulated'
|
|
|
37 "threshold", "t", 2, "numeric",
|
|
|
38 "empirical.fname", "e", 2, "character",
|
|
|
39 "perc.pass", "p", 2, "numeric",
|
|
|
40 "verbose", "v", 0, "logical", ## to set verbose on
|
|
|
41 "help", "h", 0, "logical"
|
|
|
42 ),
|
|
|
43 nc=4,
|
|
|
44 byrow=TRUE
|
|
|
45 )
|
|
|
46
|
|
|
47 opt <- getopt( spec=spec )
|
|
|
48 ##save.image( "~/work.local/tmp/ipl.feature.sel.dbg.rda" )
|
|
|
49 #set some reasonable defaults for the options that are needed,
|
|
|
50 #but were not specified.
|
|
|
51 if ( is.null(opt$verbose ) ) { opt$verbose = FALSE }
|
|
|
52 if ( is.null(opt$genes.only ) ) {
|
|
|
53 opt$genes.only <- FALSE
|
|
|
54 }
|
|
|
55
|
|
|
56 if ( is.null(opt$filter.type ) ) { opt$filter.type = 'modulated' }
|
|
|
57 if ( is.null( opt$threshold ) ) { opt$threshold=0.25 }
|
|
|
58 if ( is.null( opt$perc.pass ) ) { opt$perc.pass=1/3 }
|
|
|
59 if ( opt$perc.pass < 0 ) {
|
|
|
60 stop( "please specify a positive number for the percentage of samples that pass the filter (if applicable)" )
|
|
|
61 }
|
|
|
62 ## now set filter.type, threshold & perc.pass if an empirical result has been passed in
|
|
|
63 if ( ! is.null( opt$empirical.fname ) ) {
|
|
|
64
|
|
|
65 if ( ! file.exists( opt$empirical.fname ) ) stop( "can't file empirical result file:", opt$empirical.fname, "\n" )
|
|
|
66 ## assume this is an RData file
|
|
|
67 emp.fname.contents <- load( opt$empirical.fname )
|
|
|
68 if ( ! "opt.thresh" %in% emp.fname.contents ) stop( "no optimal threshold value found in RData file passed in\n" )
|
|
|
69 opt$threshold <- opt.thresh
|
|
|
70
|
|
|
71 if ( ! "filter.type" %in% emp.fname.contents ) stop( "no filter type value found in RData file passed in\n" )
|
|
|
72 opt$filter.type <- filter.type
|
|
|
73
|
|
|
74 if ( ! "perc.pass" %in% emp.fname.contents ) stop( "no percentage passing value found in RData file passed in\n" )
|
|
|
75 opt$perc.pass <- perc.pass
|
|
|
76 }
|
|
|
77 if ( ! opt$filter.type %in% c( 'active', 'inactive', 'modulated' ) ) stop( 'invalid filter.type specified:', opt$filter.type, "\n" )
|
|
|
78 if ( is.null( opt$output.name ) ) {
|
|
|
79 opt$output.name <- file.path( getwd(),
|
|
|
80 paste( opt$filter.type, basename( opt$data.fname ), sep="." ) )
|
|
|
81 }
|
|
|
82
|
|
|
83
|
|
|
84
|
|
|
85 data <- as.matrix( read.delim( opt$data.fname, header=T, row.names=1 , check.names=FALSE ) )
|
|
|
86 if ( opt$genes.only ) {
|
|
|
87 genes <- rownames( data )
|
|
|
88 genes <- genes[ ! grepl( "abstract|complex|family", genes ) ]
|
|
|
89 data <- data[ genes, ]
|
|
|
90 }
|
|
|
91
|
|
|
92
|
|
|
93 count.samps.threshold <- function( data,
|
|
|
94 threshold,
|
|
|
95 comparator ## must be one of lte, lt, gt, gte
|
|
|
96 ) {
|
|
|
97 filter.vect <- rep( TRUE, nrow( data ) ) ## set an initial val
|
|
|
98 if ( comparator == "lt" ) {
|
|
|
99 return( apply( data,
|
|
|
100 1,
|
|
|
101 function(x) sum( x < threshold, na.rm=T ) ) )
|
|
|
102 }
|
|
|
103 if ( comparator == "lte" ) {
|
|
|
104 return( apply( data,
|
|
|
105 1,
|
|
|
106 function(x) sum( x <= threshold, na.rm=T ) ) )
|
|
|
107 }
|
|
|
108 if ( comparator == "gte" ) {
|
|
|
109 return( apply( data,
|
|
|
110 1,
|
|
|
111 function(x) sum( x >= threshold, na.rm=T ) ) )
|
|
|
112 }
|
|
|
113 if ( comparator == "gt" ) {
|
|
|
114 return( apply( data,
|
|
|
115 1,
|
|
|
116 function(x) sum( x > threshold, na.rm=T ) ) )
|
|
|
117 }
|
|
|
118 if ( comparator == "bothe" ) {
|
|
|
119 return( apply( data,
|
|
|
120 1,
|
|
|
121 function(x) sum( abs(x) >= threshold, na.rm=T ) ) )
|
|
|
122 }
|
|
|
123 if ( comparator == "both" ) {
|
|
|
124 return( apply( data,
|
|
|
125 1,
|
|
|
126 function(x) sum( abs(x) > threshold, na.rm=T ) ) )
|
|
|
127 }
|
|
|
128 }
|
|
|
129
|
|
|
130
|
|
|
131
|
|
|
132
|
|
|
133 if ( opt$filter.type=="active" ) {
|
|
|
134 ## this is an implementation of the activity filter that was used in the original PARADIGM paper
|
|
|
135 filter.vect <- count.samps.threshold( data, opt$threshold, "gt" )
|
|
|
136 } else {
|
|
|
137 if ( opt$filter.type=="inactive" ) {
|
|
|
138 filter.vect <- count.samps.threshold( data, -opt$threshold, "lt" )
|
|
|
139 } else {
|
|
|
140 if ( opt$filter.type=="modulated" ) {
|
|
|
141 filter.vect <- count.samps.threshold( data, opt$threshold, "both" )
|
|
|
142 } else {
|
|
|
143 stop( "invalid filter.type specified: ", opt$filter.type )
|
|
|
144 }
|
|
|
145 }
|
|
|
146 }
|
|
|
147
|
|
|
148 if ( opt$perc.pass <1 ) {
|
|
|
149 filter.vect <- filter.vect > floor( ncol( data ) * opt$perc.pass )
|
|
|
150 } else {
|
|
|
151 filter.vect <- filter.vect >= opt$perc.pass
|
|
|
152 }
|
|
|
153 data <- data[ filter.vect, ]
|
|
|
154
|
|
|
155 write.table( data, opt$output.name, sep="\t", row.names=TRUE, col.names=NA, quote=FALSE )
|
|
|
156
|
