annotate compRegions.pl @ 6:088c679fb58a draft default tip

Uploaded
author antmarge
date Tue, 28 Mar 2017 21:53:57 -0400
parents 198cfbcae096
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
1 #!/usr/bin/perl -w
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
2
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
3 #Margaret Antonio 17.01.10 For GALAXY
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
4
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
5 use Data::Dumper;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
6 use strict;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
7 use Getopt::Long;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
8 use Statistics::Distributions;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
9 no warnings;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
10
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
11 #ASSIGN INPUTS TO VARIABLES USING FLAGS
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
12 our ($input1, $input2, $help,$out,$round,$l1, $l2);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
13
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
14 GetOptions(
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
15 'input1:s' => \$input1,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
16 'input2:s' => \$input2,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
17 'h' => \$help,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
18 'o:s' =>\$out,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
19 'r:i'=> \$round,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
20 'l1:s'=> \$l1,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
21 'l2:s'=> \$l2,
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
22 );
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
23
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
24 sub print_usage() {
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
25 print "\n####################################################################\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
26 print "\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
27 print "compWindows.pl: compare regions outputted by the slidingWindow tool\n\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
28 print "DESCRIPTION: Takes two slidingWindows.csv files and compares them by\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
29 print "calculating the difference in mean fitness, the pval for each gene.\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
30 print "Example: compare control vs antibiotic, where both from same strain (genome).\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
31 print "Note: For different strains/genomes, use compStrains for a gene comparison.pl\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
32
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
33 print "\nUSAGE:\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
34 print "perl compGenes.pl -d inputs/ -l cond1,cond2 -r 2 -o comp-cond1cond2_date.csv\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
35
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
36 print "\nREQUIRED:\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
37 print " -d\tDirectory containing two slidingWindow.csv files (output of\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
38 print " \tslidingWindow tool)\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
39 print " \tOR\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
40 print " \tIn the command line (without a flag), input the name(s) of\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
41 print " \ttwo files.\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
42
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
43 print "\nOPTIONAL:\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
44 print " -h\tPrints usage and exits program\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
45 print " -o\tOutput file for comparison data. Default: label1label2.csv\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
46 print " -r\tRound final output numbers to this number of decimals\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
47 print " -s\tColumn number to sort by. Default: difference of means\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
48 print " -l\tLabels for input files. Default: filenames\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
49 print " \tTwo strings, comma separated (i.e. -l expt1,expt2).\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
50 print " \tOrder should match file order.\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
51
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
52 print " \n~~~~Always check that file paths are correctly specified~~~~\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
53 print "\n##################################################################\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
54 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
55 if ($help){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
56 print_usage();
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
57 exit;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
58 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
59
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
60
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
61
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
62 #THE @files ARRAY WILL CONTAIN INPUT FILE NAMES
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
63 my @files=($input1,$input2);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
64
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
65 #GET LABELS: USE (-l) OR USE FILNEAMES AS LABELS FOR COLUMNS IN OUTPUT FILE
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
66
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
67 my @labels=($l1,$l2);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
68
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
69 #INDICES 0:start 1:end 2:mutants 3:insertions 4:TA_sites 5:ratio 6:p-value
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
70 #7:average 8:variance 9:stdev 10:stderr 11:genes 12:fit-mean
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
71
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
72 #CHECK IF REQ. VARIABLES WERE DEFINED USING FLAGS. IF NOT THEN USE DEFAULT VALUES
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
73
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
74 if (!$out) {$out="comp.".$labels[0].$labels[1].".csv"}
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
75 if (!$round){$round='%.4f'}
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
76
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
77 #READ INPUT FILES AND PUT GENE INFO FROM EACH LINE INTO HASH %all WHERE KEY=GENE_ID AND
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
78 #VALUE=GENE INFO (MEAN FITNESS, # INSERTIONS, ETC.)
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
79
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
80 my %all;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
81 my @header;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
82 for (my $i=0; $i<2; $i++){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
83 print "File #",$i+1,"\t";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
84 my $file=$files[$i];
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
85 print $file,"\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
86
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
87 open(DATA, '<', $file) or (print "Could not open '$file'\n" and print_usage() and exit);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
88
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
89 #EXTRACT THE HEADER (COLUMN NAMES) OF THE FILE AND KEEP FOR OUTPUT HEADER
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
90 #APPEND FILE NAME OR GIVEN LABEL (-l) TO COLUMN NAME SO ITS DIFFERENT FROM OTHER INPUT FILE
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
91
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
92 my $head=<DATA>;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
93 my @cols=split(',',$head);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
94 for (my $j=0;$j<scalar @cols;$j++){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
95 my $colname=$cols[$j];
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
96 chomp($colname);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
97 $cols[$j]=$colname.'-'.$labels[$i];
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
98 print $cols[$j],"\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
99 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
100 push (@header,@cols);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
101
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
102 while (my $entry = <DATA>) {
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
103 chomp $entry;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
104 my @line=split(",",$entry);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
105 if (!$line[11]){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
106 $line[11]="NA";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
107 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
108 my ($start,$end,$mutants,$insertions,$TAsites,$ratio,$pval,$avg,$var,$sd,$se,$genes)=@line;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
109 $avg=sprintf("%.3f",$avg);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
110
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
111 #if (!$line[6]){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
112 # $line[6]=0;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
113 #}
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
114 #@line=@line[0,1,2,3,4,5,6];
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
115
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
116
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
117 #PUT GENE AND INFO INTO THE HASH FOR EXISTING KEY (1ST FILE) OR CREATE NEW KEY (2ND FILE)
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
118 my $key=$start;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
119 if(!exists $all{$key}){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
120 my @info;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
121 push (@info,@line);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
122 $all{$key}=\@info;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
123 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
124 #Otherwise the window existed for the prior file, and now need to calcualte extra stats
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
125 else{
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
126 my @info=@{$all{$key}};
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
127 my ($start2,$end2,$mutants2,$insertions2,$TAsites2,$ratio2,$pval2,$avg2,$var2,$sd2,$se2,$genes2)=@info;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
128 my $diff=sprintf("$round",($avg-$avg2));
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
129 my $df=$insertions+$insertions2-2;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
130 my $tdist;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
131 my $pval;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
132
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
133 # CHECK TO MAKE SURE ALL VARIABLES IN TDIST,PVAL CALCULATIONS ARE NUMBERS AND NO
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
134 # ZEROS (0) IN THE DENOMINATOR
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
135
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
136 if ($se eq "NA" or $se2 eq "NA" or $sd eq "NA" or $sd2 eq "NA" or $insertions==0 or $insertions2==0 or $sd==0 or $sd2==0 or $df<=0){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
137 ($tdist,$pval)=(" "," ");
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
138 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
139 else{
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
140 $tdist=sqrt((($diff)/(sqrt((($sd**2)/$insertions)+(($sd2**2)/$insertions2))))**2);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
141 $pval=Statistics::Distributions::tprob($df,$tdist);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
142 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
143 push (@info,@line,$diff,$df,$tdist,$pval);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
144 $all{$key}=\@info;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
145 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
146 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
147 close DATA;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
148 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
149
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
150 #READ ALL HASH CONTENTS INTO 2D ARRAY FOR EASY SORTING AND PRINTING TO OUT FILE
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
151 my @unsorted;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
152
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
153 foreach my $entry (keys %all) {
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
154 my @info=@{$all{$entry}};
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
155 my @temp;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
156 push (@temp,@info);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
157 push (@unsorted,\@temp);
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
158 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
159
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
160 #SORT GENES BY PVALUE OR FITNESS DEPENDING ON WHICH FLAG WAS USED IN COMMANDLINE
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
161
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
162 my $sortkey=26; #default: sort by difference of means
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
163
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
164 my @sorted = sort { $b->[$sortkey] <=> $a->[$sortkey] } @unsorted;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
165
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
166 #ADD NEW FIELDS TO HEADER (COLUMN NAMES)
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
167 my $field="Mean".$labels[0].'.'.$labels[1];
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
168 push (@header,$field,"DOF","TDIST","PVALUE");
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
169
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
170 #PRINT TO OUT FILE
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
171 open OUT, '>',"$out";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
172 print OUT (join(',',@header),"\n");
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
173 foreach (@sorted){
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
174 my @woo=@{$_};
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
175 print OUT join(',',@woo),"\n";
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
176 }
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
177 close OUT;
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
178
198cfbcae096 Uploaded
antmarge
parents:
diff changeset
179