Mercurial > repos > antmarge > compstrains

comparison compStrains.pl @ 12:e627f5820272 draft

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