|
17
|
1 #!/usr/bin/perl -w
|
|
|
2
|
|
|
3 #Margaret Antonio 16.01.13
|
|
|
4
|
|
|
5 #DESCRIPTION: Takes two aggregate.pl outputs and compares them using mean difference, pval for each
|
|
|
6 #gene. Can compare, for example, 19F in glucose and TIGR4 in glucose.
|
|
|
7 #DIFFERENT GENOMES (ie. diff. strains).
|
|
|
8 #Requires CONVERSION FILE
|
|
|
9
|
|
|
10 #USAGE: perl compStrains.pl -c <conversion.csv> <options>
|
|
|
11 #[<aggregateFile1.csv aggregateFile2.csv> OR -indir <indir/>]
|
|
|
12
|
|
|
13 use Data::Dumper;
|
|
|
14 use strict;
|
|
|
15 use Getopt::Long;
|
|
|
16 use warnings;
|
|
|
17 use File::Path;
|
|
|
18 use File::Basename;
|
|
|
19 use Statistics::Distributions;
|
|
|
20
|
|
|
21 #ASSIGN INPUTS TO VARIABLES USING FLAGS
|
|
|
22 our ($input1,$input2,$out,$sortkey,$round,$l1,$l2,$cfile);
|
|
|
23 GetOptions(
|
|
|
24 'o:s' =>\$out,
|
|
|
25 's:i' => \$sortkey,
|
|
|
26 'r:i'=> \$round,
|
|
|
27 'l1:s'=> \$l1,
|
|
|
28 'l2:s'=> \$l2,
|
|
|
29 'c:s'=> \$cfile,
|
|
|
30 'input1:s'=>\$input1,
|
|
|
31 'input2:s'=>\$input2
|
|
|
32 );
|
|
|
33
|
|
|
34 sub print_usage() {
|
|
|
35 print "\n";
|
|
|
36 print "\n##################################################################\n";
|
|
|
37 print "compStrains.pl: compare genes from a tn-seq experiment\n";
|
|
|
38 print "\tfor two DIFFERENT strains/genomes using aggregate files\n";
|
|
|
39
|
|
|
40 print "\nDESCRIPTION: Takes two aggregate.pl outputs and compares them by\n";
|
|
|
41 print "calculating the difference in mean fitness.\n";
|
|
|
42
|
|
|
43 print "Example: two strains tested under same condition.\n";
|
|
|
44 print "Note: For same strains (genomes), use compGenes.pl\n";
|
|
|
45
|
|
|
46 print "\nUSAGE:\n";
|
|
|
47 print "perl compStrains.pl -c conversion.csv -d inputs/\n";
|
|
|
48
|
|
|
49 print "\nREQUIRED:\n";
|
|
|
50 print " -d\tDirectory containing all input files (files from\n";
|
|
|
51 print " \taggregate fitness script)\n";
|
|
|
52 print " \tOR\n";
|
|
|
53 print " \tIn the command line (without a flag), input the name(s) of\n";
|
|
|
54 print " \ttwo files containing aggregate gene fitness values. \n\n";
|
|
|
55 print " -c\tConversion file: two columns with homologs for both organisms\n";
|
|
|
56
|
|
|
57 print "\nOPTIONAL:\n";
|
|
|
58 print " -h\tPrints usage and exits program\n";
|
|
|
59 print " -o\tOutput file for comparison data. Default: label1label2.csv\n";
|
|
|
60 print " -s\tSort output by this index of the file (indices begin at 0).\n";
|
|
|
61 print " \tDefault: by mean\n";
|
|
|
62 print " -r\tRound final output numbers to this number of decimals\n";
|
|
|
63 print " -l\tLabels for input files. Default: filenames\n";
|
|
|
64 print " \tTwo strings, comma separated (i.e. -l expt1,expt2).\n";
|
|
|
65 print " \tOrder should match file order.\n";
|
|
|
66 print " \n~~~~Always check that file paths are correctly specified~~~~\n";
|
|
|
67 print "\n##################################################################\n";
|
|
|
68 }
|
|
|
69
|
|
|
70 #THE @files ARRAY WILL CONTAIN INPUT FILE NAMES
|
|
|
71 my @files=($input1,$input2);
|
|
|
72
|
|
|
73 #GET LABELS:
|
|
|
74 my @labels = ($l1,$l2);
|
|
|
75
|
|
|
76
|
|
|
77 #CHECK IF REQ. VARIABLES WERE DEFINED USING FLAGS. IF NOT THEN USE DEFAULT VALUES
|
|
|
78
|
|
|
79 if (!$out) {$out="comp.".$labels[0].$labels[1].".csv"}
|
|
|
80 if (!$round){$round='%.4f'}
|
|
|
81
|
|
|
82 #OPEN INPUTTED AGGREGATE GENE FILES AND STORE THEIR CONTENTS INTO TWO HASHES
|
|
|
83 #FILE1 GOES INTO HASH %ONE AND FILE2 GOES INTO HASH %TWO.
|
|
|
84
|
|
|
85 #FILE1 OPENING ---> %one WHERE KEY:VALUE IS GENE_ID:(GENE_ID,INSERTIONS,MEAN,ETC.)
|
|
|
86 my @header;
|
|
|
87 my %one;
|
|
|
88
|
|
|
89 open (F1,'<',$files[0]);
|
|
|
90
|
|
|
91 #STORE COLUMN NAMES (FIRST LINE OF FILE1) FOR HEADER AND APPEND LABELS
|
|
|
92 my $head=<F1>; #the header in the file
|
|
|
93 my @cols=split(',',$head);
|
|
|
94 @cols=@cols[0,1,2,3,4,5,6]; #get rid of blank columns
|
|
|
95 for (my $j=0;$j<scalar @cols;$j++){
|
|
|
96 $cols[$j]=$cols[$j].'-'.$labels[0]; #mark each column name with file it comes from
|
|
|
97 }
|
|
|
98 push (@header,@cols);
|
|
|
99
|
|
|
100 while (my $line=<F1>){
|
|
|
101 chomp $line;
|
|
|
102 my @info=split(",",$line);
|
|
|
103 #Only keep the first 7 columns (Ones about blanks aren't needed for comparisons)
|
|
|
104 @info=@info[0,1,2,3,4,5,6];
|
|
|
105 #Sometimes genes that don't have a gene name can't be blank, so fill with NA
|
|
|
106 if (!$info[5]){
|
|
|
107 $info[5]="NA";
|
|
|
108 }
|
|
|
109 #If there are no insertions in the column "total", then make it =0 rather than blank
|
|
|
110 if (!$info[6]){
|
|
|
111 $info[6]=0;
|
|
|
112 }
|
|
|
113 $one{$info[0]}=\@info;
|
|
|
114 }
|
|
|
115 close F1;
|
|
|
116
|
|
|
117 #FILE2 OPENING ---> %two WHERE KEY:VALUE IS GENE_ID:(GENE_ID,INSERTIONS,MEAN,ETC.)
|
|
|
118
|
|
|
119 my %two;
|
|
|
120 open (F2,'<',$files[1]);
|
|
|
121
|
|
|
122 #STORE COLUMN NAMES (FIRST LINE OF FILE2) FOR HEADER AND APPEND LABELS
|
|
|
123 $head=<F2>; #the header in the file
|
|
|
124 @cols=split(',',$head);
|
|
|
125 @cols=@cols[0,1,2,3,4,5,6]; #get rid of blank columns
|
|
|
126 for (my $j=0;$j<scalar @cols;$j++){
|
|
|
127 $cols[$j]=$cols[$j].'-'.$labels[1]; #mark each column name with file it comes from
|
|
|
128 }
|
|
|
129 push (@header,@cols);
|
|
|
130
|
|
|
131 while (my $line=<F2>){
|
|
|
132 chomp $line;
|
|
|
133 my @info=split(",",$line);
|
|
|
134 @info=@info[0,1,2,3,4,5,6];
|
|
|
135 if (!$info[5]){
|
|
|
136 $info[5]="NA";
|
|
|
137 }
|
|
|
138 if (!$info[6]){
|
|
|
139 $info[6]=0;
|
|
|
140 }
|
|
|
141 $two{$info[0]}=\@info;
|
|
|
142 }
|
|
|
143 close F2;
|
|
|
144
|
|
|
145
|
|
|
146 #READ CONVERSION FILE INTO ARRAY.
|
|
|
147 #Conversion file must have strain 1 for file 1 in column 1 (index 0) and
|
|
|
148 #strain 2 for file 2 in column 2 (index 1)
|
|
|
149 #conversion file must be tab delimited with no NA fields
|
|
|
150 #If homologs (exist then take info from hashes (%one and %two) by referring to gene_id in KEY
|
|
|
151
|
|
|
152 my @all; #store all homologs in this hash
|
|
|
153 open (CONV,'<',$cfile);
|
|
|
154 while (my $line=<CONV>){
|
|
|
155 chomp $line;
|
|
|
156 my @genes=split("\t",$line); #Array @genes will contain two genes (SP_0000,SPT_0000)
|
|
|
157 if (scalar @genes==2 and $genes[0] ne "" and $genes[1] ne ""){
|
|
|
158 my @info;
|
|
|
159 my @oneArray=@{$one{$genes[0]}};
|
|
|
160 my @twoArray=@{$two{$genes[1]}};
|
|
|
161 push (@info,@oneArray,@twoArray);
|
|
|
162 my $diff=sprintf("$round",($info[1]-$info[8]));
|
|
|
163 my $total1=$info[6];
|
|
|
164 my $total2=$info[13];
|
|
|
165 my $sd1=$info[3];
|
|
|
166 my $se1=$info[4];
|
|
|
167 my $sd2=$info[10];
|
|
|
168 my $se2=$info[11];
|
|
|
169 my $df=$total1+$total2-2;
|
|
|
170 my $tdist;
|
|
|
171 my $pval;
|
|
|
172 #TDIST, PVAL calculations with fail if standard dev, error, or counts are not real numbers
|
|
|
173 #or if 0 ends up in denominator
|
|
|
174 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){
|
|
|
175 ($tdist,$pval)=("NA","NA");
|
|
|
176 }
|
|
|
177 else{
|
|
|
178 $tdist=sqrt((($diff)/(sqrt((($sd1**2)/$total1)+(($sd2**2)/$total2))))**2);
|
|
|
179 $pval=Statistics::Distributions::tprob($df,$tdist);
|
|
|
180 }
|
|
|
181 push (@info,$diff,$df,$tdist,$pval);
|
|
|
182 push (@all,\@info);
|
|
|
183 }
|
|
|
184 }
|
|
|
185 close CONV;
|
|
|
186
|
|
|
187 #SORT THE HOMOLOGS BY THE SORTKEY OR BY DEFAULT DIFFERENCE IN MEAN FITNESSES
|
|
|
188 if (!$sortkey){
|
|
|
189 $sortkey=14; #for mean difference
|
|
|
190 }
|
|
|
191 my @sorted = sort { $b->[$sortkey] <=> $a->[$sortkey] } @all;
|
|
|
192
|
|
|
193 #FINISH THE HEADER BY ADDING COLUMN NAMES FOR MEAN-DIFF, DOF, TDIST, AND PVALUE
|
|
|
194 my $field="MeanDiff(".$labels[0].'.'.$labels[1].")";
|
|
|
195 push (@header,$field,"DOF","TDIST","PVALUE");
|
|
|
196
|
|
|
197 #PRINT MATCHED HOMOLOG INFORMATION INTO A SINGLE OUTPUT FILE
|
|
|
198 open OUT, '>',"$out";
|
|
|
199 print OUT (join(',',@header),"\n");
|
|
|
200 foreach (@sorted){
|
|
|
201 my @woo=@{$_};
|
|
|
202 print OUT join(',',@woo),"\n";
|
|
|
203 }
|
|
|
204
|
|
|
205 close OUT;
|
|
|
206
|
|
|
207
|
