annotate t_test_two_samples.pl @ 0:3d24b35de20a draft

Uploaded
author devteam
date Wed, 21 Aug 2013 10:58:59 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
3d24b35de20a Uploaded
devteam
parents:
diff changeset
1 # A program to implement the non-pooled t-test for two samples where the alternative hypothesis is two-sided or one-sided.
3d24b35de20a Uploaded
devteam
parents:
diff changeset
2 # The first input file is a TABULAR format file representing the first sample and consisting of one column only.
3d24b35de20a Uploaded
devteam
parents:
diff changeset
3 # The second input file is a TABULAR format file representing the first sample nd consisting of one column only.
3d24b35de20a Uploaded
devteam
parents:
diff changeset
4 # The third input is the sidedness of the t-test: either two-sided or, one-sided with m1 less than m2 or,
3d24b35de20a Uploaded
devteam
parents:
diff changeset
5 # one-sided with m1 greater than m2.
3d24b35de20a Uploaded
devteam
parents:
diff changeset
6 # The fourth input is the equality status of the standard deviations of both populations
3d24b35de20a Uploaded
devteam
parents:
diff changeset
7 # The output file is a TXT file representing the result of the two sample t-test.
3d24b35de20a Uploaded
devteam
parents:
diff changeset
8
3d24b35de20a Uploaded
devteam
parents:
diff changeset
9 use strict;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
10 use warnings;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
11
3d24b35de20a Uploaded
devteam
parents:
diff changeset
12 #variable to handle the motif information
3d24b35de20a Uploaded
devteam
parents:
diff changeset
13 my $motif;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
14 my $motifName = "";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
15 my $motifNumber = 0;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
16 my $totalMotifsNumber = 0;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
17 my @motifNamesArray = ();
3d24b35de20a Uploaded
devteam
parents:
diff changeset
18
3d24b35de20a Uploaded
devteam
parents:
diff changeset
19 # check to make sure having correct files
3d24b35de20a Uploaded
devteam
parents:
diff changeset
20 my $usage = "usage: non_pooled_t_test_two_samples_galaxy.pl [TABULAR.in] [TABULAR.in] [testSidedness] [standardDeviationEquality] [TXT.out] \n";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
21 die $usage unless @ARGV == 5;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
22
3d24b35de20a Uploaded
devteam
parents:
diff changeset
23 #get the input arguments
3d24b35de20a Uploaded
devteam
parents:
diff changeset
24 my $firstSampleInputFile = $ARGV[0];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
25 my $secondSampleInputFile = $ARGV[1];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
26 my $testSidedness = $ARGV[2];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
27 my $standardDeviationEquality = $ARGV[3];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
28 my $outputFile = $ARGV[4];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
29
3d24b35de20a Uploaded
devteam
parents:
diff changeset
30 #open the input files
3d24b35de20a Uploaded
devteam
parents:
diff changeset
31 open (INPUT1, "<", $firstSampleInputFile) || die("Could not open file $firstSampleInputFile \n");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
32 open (INPUT2, "<", $secondSampleInputFile) || die("Could not open file $secondSampleInputFile \n");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
33 open (OUTPUT, ">", $outputFile) || die("Could not open file $outputFile \n");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
34
3d24b35de20a Uploaded
devteam
parents:
diff changeset
35
3d24b35de20a Uploaded
devteam
parents:
diff changeset
36 #variables to store the name of the R script file
3d24b35de20a Uploaded
devteam
parents:
diff changeset
37 my $r_script;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
38
3d24b35de20a Uploaded
devteam
parents:
diff changeset
39 # R script to implement the two-sample test on the motif frequencies in upstream flanking region
3d24b35de20a Uploaded
devteam
parents:
diff changeset
40 #construct an R script file and save it in the same directory where the perl file is located
3d24b35de20a Uploaded
devteam
parents:
diff changeset
41 $r_script = "non_pooled_t_test_two_samples.r";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
42
3d24b35de20a Uploaded
devteam
parents:
diff changeset
43 open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
44 print Rcmd "
3d24b35de20a Uploaded
devteam
parents:
diff changeset
45 sampleTable1 <- read.table(\"$firstSampleInputFile\", header=FALSE);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
46 sample1 <- sampleTable1[, 1];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
47
3d24b35de20a Uploaded
devteam
parents:
diff changeset
48 sampleTable2 <- read.table(\"$secondSampleInputFile\", header=FALSE);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
49 sample2 <- sampleTable2[, 1];
3d24b35de20a Uploaded
devteam
parents:
diff changeset
50
3d24b35de20a Uploaded
devteam
parents:
diff changeset
51 testSideStatus <- \"$testSidedness\";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
52 STEqualityStatus <- \"$standardDeviationEquality\";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
53
3d24b35de20a Uploaded
devteam
parents:
diff changeset
54 #open the output a text file
3d24b35de20a Uploaded
devteam
parents:
diff changeset
55 sink(file = \"$outputFile\");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
56
3d24b35de20a Uploaded
devteam
parents:
diff changeset
57 #check if the t-test is two-sided
3d24b35de20a Uploaded
devteam
parents:
diff changeset
58 if (testSideStatus == \"two-sided\"){
3d24b35de20a Uploaded
devteam
parents:
diff changeset
59
3d24b35de20a Uploaded
devteam
parents:
diff changeset
60 #check if the standard deviations are equal in both populations
3d24b35de20a Uploaded
devteam
parents:
diff changeset
61 if (STEqualityStatus == \"equal\"){
3d24b35de20a Uploaded
devteam
parents:
diff changeset
62 #two-sample t-test where standard deviations are assumed to be unequal, the test is two-sided
3d24b35de20a Uploaded
devteam
parents:
diff changeset
63 testResult <- t.test(sample1, sample2, var.equal = TRUE);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
64 } else{
3d24b35de20a Uploaded
devteam
parents:
diff changeset
65 #two-sample t-test where standard deviations are assumed to be unequal, the test is two-sided
3d24b35de20a Uploaded
devteam
parents:
diff changeset
66 testResult <- t.test(sample1, sample2, var.equal = FALSE);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
67 }
3d24b35de20a Uploaded
devteam
parents:
diff changeset
68 } else{ #the t-test is one sided
3d24b35de20a Uploaded
devteam
parents:
diff changeset
69
3d24b35de20a Uploaded
devteam
parents:
diff changeset
70 #check if the t-test is two-sided with m1 < m2
3d24b35de20a Uploaded
devteam
parents:
diff changeset
71 if (testSideStatus == \"one-sided:_m1_less_than_m2\"){
3d24b35de20a Uploaded
devteam
parents:
diff changeset
72
3d24b35de20a Uploaded
devteam
parents:
diff changeset
73 #check if the standard deviations are equal in both populations
3d24b35de20a Uploaded
devteam
parents:
diff changeset
74 if (STEqualityStatus == \"equal\"){
3d24b35de20a Uploaded
devteam
parents:
diff changeset
75 #two-sample t-test where standard deviations are assumed to be unequal, the test is one-sided: Halt: m1 < m2
3d24b35de20a Uploaded
devteam
parents:
diff changeset
76 testResult <- t.test(sample1, sample2, var.equal = TRUE, alternative = \"less\");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
77 } else{
3d24b35de20a Uploaded
devteam
parents:
diff changeset
78 #two-sample t-test where standard deviations are assumed to be unequal, the test is one-sided: Halt: m1 < m2
3d24b35de20a Uploaded
devteam
parents:
diff changeset
79 testResult <- t.test(sample1, sample2, var.equal = FALSE, alternative = \"less\");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
80 }
3d24b35de20a Uploaded
devteam
parents:
diff changeset
81 } else{ #the t-test is one-sided with m1 > m2
3d24b35de20a Uploaded
devteam
parents:
diff changeset
82 #check if the standard deviations are equal in both populations
3d24b35de20a Uploaded
devteam
parents:
diff changeset
83 if (STEqualityStatus == \"equal\"){
3d24b35de20a Uploaded
devteam
parents:
diff changeset
84 #two-sample t-test where standard deviations are assumed to be unequal, the test is one-sided: Halt: m1 < m2
3d24b35de20a Uploaded
devteam
parents:
diff changeset
85 testResult <- t.test(sample1, sample2, var.equal = TRUE, alternative = \"greater\");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
86 } else{
3d24b35de20a Uploaded
devteam
parents:
diff changeset
87 #two-sample t-test where standard deviations are assumed to be unequal, the test is one-sided: Halt: m1 < m2
3d24b35de20a Uploaded
devteam
parents:
diff changeset
88 testResult <- t.test(sample1, sample2, var.equal = FALSE, alternative = \"greater\");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
89 }
3d24b35de20a Uploaded
devteam
parents:
diff changeset
90 }
3d24b35de20a Uploaded
devteam
parents:
diff changeset
91 }
3d24b35de20a Uploaded
devteam
parents:
diff changeset
92
3d24b35de20a Uploaded
devteam
parents:
diff changeset
93 #save the output of the t-test into the output text file
3d24b35de20a Uploaded
devteam
parents:
diff changeset
94 testResult;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
95
3d24b35de20a Uploaded
devteam
parents:
diff changeset
96 #close the output text file
3d24b35de20a Uploaded
devteam
parents:
diff changeset
97 sink();
3d24b35de20a Uploaded
devteam
parents:
diff changeset
98
3d24b35de20a Uploaded
devteam
parents:
diff changeset
99 #eof" . "\n";
3d24b35de20a Uploaded
devteam
parents:
diff changeset
100
3d24b35de20a Uploaded
devteam
parents:
diff changeset
101 close Rcmd;
3d24b35de20a Uploaded
devteam
parents:
diff changeset
102
3d24b35de20a Uploaded
devteam
parents:
diff changeset
103 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
3d24b35de20a Uploaded
devteam
parents:
diff changeset
104
3d24b35de20a Uploaded
devteam
parents:
diff changeset
105 #close the input and output files
3d24b35de20a Uploaded
devteam
parents:
diff changeset
106 close(OUTPUT);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
107 close(INPUT2);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
108 close(INPUT1);
3d24b35de20a Uploaded
devteam
parents:
diff changeset
109