|
1
|
1 #!/usr/bin/perl -w
|
|
|
2
|
|
|
3 #Margaret Antonio 17.01.06 without essentiality. For original, unmodified aggregate.pl output (where unique insertions are not included).
|
|
|
4
|
|
|
5 use Getopt::Long;
|
|
|
6 use Statistics::Distributions;
|
|
|
7 use strict;
|
|
|
8 use autodie;
|
|
|
9 no warnings;
|
|
|
10
|
|
|
11 #no warnings;
|
|
|
12
|
|
|
13 #ASSIGN INPUTS TO VARIABLES USING FLAGS
|
|
|
14 our ($input1, $input2, $h, $out, $sortkey, $round, $l1, $l2, $outfile);
|
|
|
15
|
|
|
16 GetOptions(
|
|
|
17 'input1:s' => \$input1,
|
|
|
18 'input2:s' => \$input2,
|
|
|
19 'h' => \$h,
|
|
|
20 'o:s' =>\$outfile,
|
|
|
21 'r:i'=> \$round,
|
|
|
22
|
|
|
23 );
|
|
|
24
|
|
|
25 sub print_usage() {
|
|
|
26 print "\n####################################################################\n";
|
|
|
27
|
|
|
28 print "compGenes: compare genes for an organism under different conditions\n\n";
|
|
|
29 print "DESCRIPTION: Takes two aggregate.pl outputs and compares them by\n";
|
|
|
30 print "calculating the difference in mean fitness for each gene.\n";
|
|
|
31 print "Example: compare organism in presence of control vs antibiotic.\n";
|
|
|
32 print "Note: For different strains/genomes, use compStrains.pl\n";
|
|
|
33
|
|
|
34 print "\nUSAGE:";
|
|
|
35 print "perl compGenes.pl -d inputs/ \n\n";
|
|
|
36
|
|
|
37 print "REQUIRED:\n";
|
|
|
38 print " -d\tDirectory containing all input files (files from\n";
|
|
|
39 print " \taggregate script)\n";
|
|
|
40 print " \tOR\n";
|
|
|
41 print " \tIn the command line (without a flag), input the name(s) of\n";
|
|
|
42 print " \ttwo files containing aggregate gene fitness values. \n\n";
|
|
|
43
|
|
|
44 print "OPTIONAL:\n";
|
|
|
45 print " -h\tPrints usage and exits program\n";
|
|
|
46 print " -o\tOutput file for comparison data. Default: label1label2.csv\n";
|
|
|
47 print " -r\tRound final output numbers to this number of decimals\n";
|
|
|
48 print " -l1\tLabels for input files. Default: filenames\n";
|
|
|
49 print " \tTwo strings, comma separated (i.e. -l expt1,expt2).\n";
|
|
|
50 print " \tOrder should match file order.\n";
|
|
|
51
|
|
|
52 print " \n~~~~Always check that file paths are correctly specified~~~~\n";
|
|
|
53 print "\n##################################################################\n";
|
|
|
54
|
|
|
55 }
|
|
|
56
|
|
|
57 # Check if help needed or if improper inputs
|
|
|
58
|
|
|
59 if ($h){
|
|
|
60 print_usage();
|
|
|
61 exit;
|
|
|
62 }
|
|
|
63
|
|
|
64 if (! $outfile){
|
|
|
65 $outfile="geneComp.csv"
|
|
|
66 }
|
|
|
67 $round = '%.4f';
|
|
|
68
|
|
|
69 my @files=($input1,$input2);
|
|
|
70
|
|
|
71 #GET LABELS: USE (-l) OR USE FILNEAMES AS LABELS FOR COLUMNS IN OUTPUT FILE
|
|
|
72
|
|
|
73 my @labels=($l1,$l2);
|
|
|
74
|
|
|
75 #CHECK IF REQ. VARIABLES WERE DEFINED USING FLAGS. IF NOT THEN USE DEFAULT VALUES
|
|
|
76
|
|
|
77 #$out="comp.".$labels[0].$labels[1].".csv";
|
|
|
78
|
|
|
79 #READ INPUT FILES AND PUT GENE INFO FROM EACH LINE INTO HASH %all WHERE KEY=GENE_ID AND
|
|
|
80 #VALUE=GENE INFO (MEAN FITNESS, # INSERTIONS, ETC.)
|
|
|
81
|
|
|
82 my %all;
|
|
|
83 my @header;
|
|
|
84
|
|
|
85 #OPEN TWO COMPARISON FILES, ONE PER LOOP
|
|
|
86 for (my $i=0; $i<2; $i++){
|
|
|
87 print "File #",$i+1,"\t";
|
|
|
88 my $file=$files[$i];
|
|
|
89 print $file,"\n";
|
|
|
90
|
|
|
91 open(DATA, '<', $file) or die "Could not open '$file'\n";
|
|
|
92
|
|
|
93 #EXTRACT THE HEADER (COLUMN NAMES) OF THE FILE AND KEEP FOR OUTPUT HEADER
|
|
|
94 #APPEND FILE NAME OR GIVEN LABEL (-l) TO COLUMN NAME SO ITS DIFFERENT FROM OTHER INPUT FILE
|
|
|
95
|
|
|
96 my $head=<DATA>;
|
|
|
97 my @temp = split("\n",$head);
|
|
|
98 $head = $temp[0];
|
|
|
99
|
|
|
100 my @cols=split(',',$head);
|
|
|
101 @cols = @cols[0,1,2,3,4,5];
|
|
|
102 for (my $j=0;$j<scalar @cols;$j++){
|
|
|
103 $cols[$j]=$cols[$j].'-'.$labels[$i];
|
|
|
104 }
|
|
|
105 push (@header,@cols);
|
|
|
106 while (my $entry = <DATA>) {
|
|
|
107 chomp $entry;
|
|
|
108 my @line=split(",",$entry);
|
|
|
109 if (!$line[5]){
|
|
|
110 $line[5]="NA";
|
|
|
111 }
|
|
|
112 if (!$line[6]){
|
|
|
113 $line[6]=0;
|
|
|
114 }
|
|
|
115 @line=@line[0,1,2,3,4,5];
|
|
|
116 my $gene=$line[0];
|
|
|
117 chomp($gene);
|
|
|
118
|
|
|
119 #PUT GENE AND INFO INTO THE HASH FOR EXISTING KEY (1ST FILE) OR CREATE NEW KEY (2ND FILE)
|
|
|
120
|
|
|
121 if (!exists $all{$gene}){
|
|
|
122 my @info;
|
|
|
123 push (@info,@line);
|
|
|
124 $all{$gene}=\@info;
|
|
|
125 }
|
|
|
126 else{
|
|
|
127 my @info=@{$all{$gene}};
|
|
|
128 push (@info,@line);
|
|
|
129 my $diff=sprintf("$round",($info[1]-$info[7]));
|
|
|
130 my $total1=$info[2];
|
|
|
131 my $total2=$info[8];
|
|
|
132 my $sd1=$info[3];
|
|
|
133 my $se1=$info[4];
|
|
|
134 my $sd2=$info[9];
|
|
|
135 my $se2=$info[10];
|
|
|
136 my $df=$total1+$total2-2;
|
|
|
137 my $tdist;
|
|
|
138 my $pval;
|
|
|
139
|
|
|
140 # CHECK TO MAKE SURE ALL VARIABLES IN TDIST,PVAL CALCULATIONS ARE NUMBERS AND NO
|
|
|
141 # ZEROS (0) IN THE DENOMINATOR
|
|
|
142
|
|
|
143 if ($se1 eq "X" or $se2 eq "X" or $sd1 eq "X" or $sd2 eq "X" or $total1==0 or $total2==0 or $sd1==0 or $sd2==0 or $df<=0){
|
|
|
144 ($tdist,$pval)=("NA","NA");
|
|
|
145 }
|
|
|
146 else{
|
|
|
147 $tdist=sqrt((($diff)/(sqrt((($sd1**2)/$total1)+(($sd2**2)/$total2))))**2);
|
|
|
148
|
|
|
149 $pval=Statistics::Distributions::tprob($df,$tdist);
|
|
|
150 }
|
|
|
151 push (@info,$diff,$df,$tdist,$pval);
|
|
|
152 $all{$gene}=\@info;
|
|
|
153 }
|
|
|
154 }
|
|
|
155 close DATA;
|
|
|
156 }
|
|
|
157
|
|
|
158 #READ ALL HASH CONTENTS INTO 2D ARRAY FOR EASY SORTING AND PRINTING TO OUT FILE
|
|
|
159 my @unsorted;
|
|
|
160
|
|
|
161 foreach my $entry (keys %all) {
|
|
|
162 my @info=@{$all{$entry}};
|
|
|
163 my @temp;
|
|
|
164 push (@temp,@info);
|
|
|
165 push (@unsorted,\@temp);
|
|
|
166 }
|
|
|
167
|
|
|
168 #SORT GENES BY PVALUE OR FITNESS DEPENDING ON WHICH FLAG WAS USED IN COMMANDLINE
|
|
|
169
|
|
|
170 $sortkey=15; #default: sort by difference of means
|
|
|
171
|
|
|
172 my @sorted = sort { $a->[$sortkey] <=> $b->[$sortkey] } @unsorted;
|
|
|
173
|
|
|
174 #ADD NEW FIELDS TO HEADER (COLUMN NAMES)
|
|
|
175 my $field="Mean".$labels[0].'.'.$labels[1];
|
|
|
176 push (@header,$field,"DOF","TDIST","PVALUE");
|
|
|
177
|
|
|
178 #PRINT TO OUT FILE
|
|
|
179 open OUT, '>',"$outfile";
|
|
|
180 print OUT (join(',',@header),"\n");
|
|
|
181 foreach (@sorted){
|
|
|
182 my @woo=@{$_};
|
|
|
183 print OUT join(',',@woo),"\n";
|
|
|
184 }
|
|
|
185 close OUT;
|
|
|
186
|
|
|
187
|
