Mercurial > repos > pmac > normalizemediannegcntrl

comparison normalizeMedianNegCntrl.pl @ 0:087ab474dfa4 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author pmac
date Wed, 01 Jun 2016 03:52:52 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:087ab474dfa4
1 ###############################################################################
2 # In summary this script normalizes each plate of input data to each of its negative
3 # controls.
4 #
5 # Args:
6 # input.data.frame contains one or more columns of plate data. The first column
7 # must contain the well annotation ie A01, A02 etc. Each subsequent column represents
8 # a plate of data where the column name is the plate name.
9 #
10 # plate.conf.data.frame is a file of the following format:
11 # well, type (syntax must be either "poscontr" or "negcontr" or "empty"), name (name of control)
12 #
13 # Returns:
14 # plate values normalized to each negative control type per plate.
15 #
16 # Author: kate gould
17 ###############################################################################
18
19 use strict;
20 use warnings;
21 use IO::Handle;
22 use File::Temp qw/ tempfile tempdir /;
23 my $tdir = tempdir( CLEANUP => 0 );
24
25 # check to make sure having correct input and output files
26 my $usage = "usage: normalizeMedianNegCntrl.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] \n";
27 die $usage unless @ARGV == 3;
28
29 #get the input arguments
30 my $plateData = $ARGV[0];
31 my $plateConfig = $ARGV[1];
32 my $normMedNeg = $ARGV[2];
33
34 #open the input files
35 open (INPUT1, "<", $plateData) || die("Could not open file $plateData \n");
36 open (INPUT2, "<", $plateConfig) || die("Could not open file $plateConfig \n");
37 open (OUTPUT1, ">", $normMedNeg) || die("Could not open file $normMedNeg \n");
38
39 #variable to store the name of the R script file
40 my $r_script;
41
42 # R script to implement the calcualtion of q-values based on multiple simultaneous tests p-values
43 # construct an R script file and save it in a temp directory
44 chdir $tdir;
45 $r_script = "normalizeMedianNegCntrl.r";
46
47 open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
48 print Rcmd "
49 #options(show.error.messages = FALSE);
50 header <- read.table(\"$plateData\", head=F, sep=\"\\t\", comment=\"\", nrows=1);
51 input.data.frame <- read.table(\"$plateData\", head=F, sep=\"\\t\", skip=1);
52 plate.conf.data.frame <- read.table(\"$plateConfig\", head=T, sep=\"\\t\", comment=\"\");
53
54 normalise <- function(x, y){if (!is.na(x)) return (x/y) else (return (NA))}
55
56 # assumed second column is type ie negcontr and third column is name of control
57 negative.controls.list <- unique(plate.conf.data.frame[which (plate.conf.data.frame[,2] == \"negcontr\"),][,3]);
58
59 normalized.data.frame <- data.frame(well=input.data.frame[,1], stringsAsFactors=FALSE);
60 column.name.list <- c();
61 for (negative.control in negative.controls.list){
62 negative.control.wells <- plate.conf.data.frame[which (plate.conf.data.frame[3] == negative.control),][,1];
63
64 for (i in 2:length(colnames(input.data.frame))){
65 column.name <- paste(colnames(input.data.frame)[i], \".normalised.to.\", negative.control, sep=\"\");
66 column.name.list <- append(column.name.list, column.name);
67
68 negative.control.values <- input.data.frame[((input.data.frame[,1] %in% negative.control.wells)&(!(is.na(input.data.frame[,i])))),][,i];
69 median.negative.value <- median(negative.control.values);
70
71 normalised.to.median.negative.control.value <- round(sapply(input.data.frame[,i], FUN=normalise, median.negative.value), 2);
72 normalized.data.frame <- cbind(normalized.data.frame, normalised.to.median.negative.control.value);
73 }
74 }
75 header <- c(header[1], column.name.list)
76 write.table(header, file=\"$normMedNeg\", quote=F, sep=\"\\t\", row.names=F, col.names=F);
77 write.table(normalized.data.frame, file=\"$normMedNeg\", quote=F, sep=\"\\t\", row.names=F, col.names=F, append=T);
78 #eof\n";
79
80 close Rcmd;
81
82 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
83
84 #close the input and output files
85 close(OUTPUT1);
86 close(INPUT1);
87 close(INPUT2);