|
0
|
1 #!usr/bin/perl
|
|
|
2
|
|
|
3 ### Perl modules
|
|
|
4 use warnings;
|
|
|
5 use strict;
|
|
|
6 use Getopt::Long qw(GetOptions); #Creation of script options
|
|
|
7 use Pod::Usage qw(pod2usage); #Creation of script options
|
|
|
8
|
|
|
9 #Personnal packages
|
|
|
10 use FindBin ; ## Allows you to locate the directory of original perl script
|
|
|
11 #use lib $FindBin::Bin;
|
|
|
12 use lib "$FindBin::Bin/lib";
|
|
|
13 use IonFiltration;
|
|
|
14
|
|
|
15 my ($file, $mass_file, $opt, $dataMatrix, $combined_DMVM, $repres_opt, $rt_threshold, $mass_threshold, $output_sif, $output_tabular, $correl_threshold, $intensity_threshold, $intensity_pourc); #Options to complete
|
|
|
16
|
|
|
17 ########################
|
|
|
18 ### Options and help ###
|
|
|
19 ########################
|
|
|
20
|
|
|
21 GetOptions("f=s"=>\$file, "m=s"=>\$mass_file, "o=s"=>\$opt, "d=s"=>\$dataMatrix, "v=s"=>\$combined_DMVM, "r=s"=>\$repres_opt, "rt=f"=>\$rt_threshold, "mass=f"=>\$mass_threshold, "output_sif=s"=>\$output_sif, "output_tabular=s"=>\$output_tabular, "correl=s"=>\$correl_threshold, "IT=f"=>\$intensity_threshold, "IP=f"=>\$intensity_pourc) or pod2usage(2);
|
|
|
22
|
|
|
23 ### Check required parameters :
|
|
|
24 pod2usage({-message=>q{Mandatory argument '-f' is missing}, -exitval=>1, -verbose=>0}) unless $file;
|
|
|
25 #pod2usage({-message=>q{Mandatory argument '-m' is missing}, -exitval=>1, -verbose=>0}) unless $mass_file;
|
|
|
26 pod2usage({-message=>q{Mandatory argument '-o' is missing. It correspond to the grouping method for analytical correlation groups formation.
|
|
|
27 #It should be a number (1 ; 2 or 3) :
|
|
|
28 # 1 : Don't take into acount mass information (only RT) ;
|
|
|
29 # 2 : Check that all mass differences are include in a specific list and taking into acount RT information
|
|
|
30 # 3 : Check that all mass differences are include in a specific list, ignoring RT information
|
|
|
31 #To use the tool without takinf into account mass and RT information, use option 1 and define the RT threshold to 999999999.}, -exitval=>1, -verbose=>0}) unless $opt;
|
|
|
32 pod2usage({-message=>q{Mandatory argument '-r' is missing. It correspond to the group representent choosing method for analytical correlation groups formation.
|
|
|
33 It should be one of the 3 options below :
|
|
|
34 "mass" : choose the ion with the highest mass as the representant
|
|
|
35 "intensity" : choose the ion with the highest intensity as the representant
|
|
|
36 "mixt" : choose the ion with the highest (mass^2 * intensity) as the representant
|
|
|
37 "max_intensity_max_mass" : choose tha ion witht he highest intenisty among the 5 most intense ions of the group}, -exitval=>1, -verbose=>0}) unless $repres_opt;
|
|
|
38 pod2usage({-message=>q{Mandatory argument '-d' is missing}, -exitval=>1, -verbose=>0}) unless $dataMatrix;
|
|
|
39 pod2usage({-message=>q{Mandatory argument '-v' is missing}, -exitval=>1, -verbose=>0}) unless $combined_DMVM;
|
|
|
40 #pod2usage({-message=>q{Mandatory argument '-rt' is missing}, -exitval=>1, -verbose=>0}) unless $rt_threshold;
|
|
|
41 #pod2usage({-message=>q{Mandatory argument '-mass' is missing}, -exitval=>1, -verbose=>0}) unless $mass_threshold;
|
|
|
42 pod2usage({-message=>q{Mandatory argument '-correl' is missing}, -exitval=>1, -verbose=>0}) unless $correl_threshold;
|
|
|
43 pod2usage({-message=>q{Mandatory argument '-output_tabular' is missing}, -exitval=>1, -verbose=>0}) unless $output_tabular;
|
|
|
44 pod2usage({-message=>q{Mandatory argument '-output_sif' is missing}, -exitval=>1, -verbose=>0}) unless $output_sif;
|
|
|
45
|
|
|
46
|
|
|
47 #if(($opt != 1) && ($opt != 2) && ($opt != 3)){
|
|
|
48 # print "you must indicate \"1\", \"2\" or \"3\" for the --o otpion\n";
|
|
|
49 # exit;
|
|
|
50 #}
|
|
|
51
|
|
|
52
|
|
|
53
|
|
|
54 if(($repres_opt ne "mass") && ($repres_opt ne "intensity") && ($repres_opt ne "mixt") && ($repres_opt ne "max_intensity_max_mass")){
|
|
|
55 print "you must indicate \"mass\", \"intensity\", \"mix\" or \"max_intensity_max_mass\" for the --r otpion\n";
|
|
|
56 exit;
|
|
|
57 }
|
|
|
58
|
|
|
59
|
|
|
60
|
|
|
61 #########################################################################
|
|
|
62 #### Création of a hash containing all adduits and fragments possible ###
|
|
|
63 #########################################################################
|
|
|
64
|
|
|
65 my %hmass;
|
|
|
66 if($opt != 1){
|
|
|
67 %hmass = IonFiltration::MassCollecting($mass_file);
|
|
|
68
|
|
|
69 }
|
|
|
70
|
|
|
71 my $refhmass = \%hmass;
|
|
|
72
|
|
|
73 print "Création of a hash containing all adduits and fragments possible\n";
|
|
|
74
|
|
|
75
|
|
|
76 ########################################################
|
|
|
77 ### Creation of a sif table + correlation filtration ###
|
|
|
78 ########################################################
|
|
|
79
|
|
|
80 my %hrtmz;
|
|
|
81 ($output_sif, %hrtmz) = IonFiltration::sifTableCreation($file, $output_sif, $opt, $rt_threshold, $mass_threshold, $correl_threshold, $dataMatrix, $output_tabular, $combined_DMVM, $repres_opt, $intensity_threshold, $intensity_pourc, \%hmass);
|
|
|
82 print "Creation of a sif table + correlation filtration done\n";
|
|
|
83
|
|
|
84
|
|
|
85 ######################################################
|
|
|
86 ### Analytic correlation filtrering follow options ###
|
|
|
87 ######################################################
|
|
|
88
|
|
|
89 my %hheader_file;
|
|
|
90 my %hduplicate;
|
|
|
91
|
|
|
92 my %hcorrelgroup;
|
|
|
93 my $groupct=1;
|
|
|
94
|
|
|
95 my $linenb3=0;
|
|
|
96 my %hheader_line;
|
|
|
97
|
|
|
98
|
|
|
99
|
|
|
100 open (F1, $output_sif) or die "Impossible to open $output_sif\n";
|
|
|
101
|
|
|
102 while(my $line = <F1>){
|
|
|
103 my $count=0;
|
|
|
104 chomp $line;
|
|
|
105 my @tline = split(/\t/, $line);
|
|
|
106 my $a = $tline[0];
|
|
|
107 my $b = $tline[2];
|
|
|
108
|
|
|
109 my $amass=$hrtmz{$a}{mz};
|
|
|
110 my $atemp=$hrtmz{$a}{rt};
|
|
|
111 my $bmass= $hrtmz{$b}{mz};
|
|
|
112 my $btemp=$hrtmz{$b}{rt};
|
|
|
113 print "YY : $a ==> $amass ; $b ==> $bmass\n";
|
|
|
114 my $diff = $amass-$bmass;
|
|
|
115 $diff = abs($diff);
|
|
|
116
|
|
|
117 ### Option 1: Don't take into acount mass information ###
|
|
|
118
|
|
|
119 if($opt == 1){
|
|
|
120 my $btplus = $btemp + $rt_threshold;
|
|
|
121 my $btmoins = $btemp - $rt_threshold;
|
|
|
122 if(($btmoins <= $atemp) && ($atemp <= $btplus)){
|
|
|
123 foreach my $k (keys %hcorrelgroup){
|
|
|
124 if((defined($hcorrelgroup{$k}{$a})) || (defined($hcorrelgroup{$k}{$b}))){
|
|
|
125 $hcorrelgroup{$k}{$a}=1;
|
|
|
126 $hcorrelgroup{$k}{$b}=1;
|
|
|
127 $count++;
|
|
|
128 last;
|
|
|
129 }
|
|
|
130 }
|
|
|
131 if($count == 0){
|
|
|
132 my $groupnb="group".$groupct;
|
|
|
133 $hcorrelgroup{$groupnb}{$a}=1;
|
|
|
134 $hcorrelgroup{$groupnb}{$b}=1;
|
|
|
135 $groupct ++;
|
|
|
136 }
|
|
|
137 }
|
|
|
138 }
|
|
|
139
|
|
|
140
|
|
|
141
|
|
|
142 ### Option 2: Check that all mass differences are include in a specific list taking into account RT information ###
|
|
|
143
|
|
|
144 elsif($opt == 2){
|
|
|
145
|
|
|
146 my $print = 0;
|
|
|
147 foreach my $s (keys %{$refhmass}){
|
|
|
148 foreach my $r (keys %{$refhmass->{$s}}){
|
|
|
149 my $rm = $r - $mass_threshold;
|
|
|
150 my $rp = $r + $mass_threshold;
|
|
|
151 if(($diff <= $rp) && ($diff >= $rm)){
|
|
|
152 if($print == 0){
|
|
|
153 my $btplus = $btemp + $rt_threshold;
|
|
|
154 my $btmoins = $btemp - $rt_threshold;
|
|
|
155
|
|
|
156 if(($btmoins <= $atemp) && ($atemp <= $btplus)){
|
|
|
157 foreach my $k (keys %hcorrelgroup){
|
|
|
158 if((defined($hcorrelgroup{$k}{$a})) || (defined($hcorrelgroup{$k}{$b}))){
|
|
|
159 $hcorrelgroup{$k}{$a}=1;
|
|
|
160 $hcorrelgroup{$k}{$b}=1;
|
|
|
161 $count++;
|
|
|
162 last;
|
|
|
163 }
|
|
|
164 }
|
|
|
165 if($count == 0){
|
|
|
166 my $groupnb="group".$groupct;
|
|
|
167 $hcorrelgroup{$groupnb}{$a}=1;
|
|
|
168 $hcorrelgroup{$groupnb}{$b}=1;
|
|
|
169 $groupct ++;
|
|
|
170 }
|
|
|
171 $print = 1;
|
|
|
172 }
|
|
|
173 }
|
|
|
174 }
|
|
|
175 }
|
|
|
176 }
|
|
|
177 }
|
|
|
178
|
|
|
179
|
|
|
180 ### Option 3: Check that all mass differences are include in a specific list, ignoring RT information ###
|
|
|
181
|
|
|
182 elsif($opt == 3){
|
|
|
183
|
|
|
184 my $print = 0;
|
|
|
185 foreach my $s (keys %{$refhmass}){
|
|
|
186 foreach my $r (keys %{$refhmass->{$s}}){
|
|
|
187 my $rm = $r - $mass_threshold;
|
|
|
188 my $rp = $r + $mass_threshold;
|
|
|
189 if(($diff <= $rp) && ($diff >= $rm)){
|
|
|
190 if($print == 0){
|
|
|
191
|
|
|
192 foreach my $k (keys %hcorrelgroup){
|
|
|
193 if((defined($hcorrelgroup{$k}{$a})) || (defined($hcorrelgroup{$k}{$b}))){
|
|
|
194 $hcorrelgroup{$k}{$a}=1;
|
|
|
195 $hcorrelgroup{$k}{$b}=1;
|
|
|
196 $count++;
|
|
|
197 last;
|
|
|
198 }
|
|
|
199 }
|
|
|
200 if($count == 0){
|
|
|
201 my $groupnb="group".$groupct;
|
|
|
202 $hcorrelgroup{$groupnb}{$a}=1;
|
|
|
203 $hcorrelgroup{$groupnb}{$b}=1;
|
|
|
204 $groupct ++;
|
|
|
205 }
|
|
|
206 $print = 1;
|
|
|
207 }
|
|
|
208 }
|
|
|
209 }
|
|
|
210 }
|
|
|
211 }
|
|
|
212 }
|
|
|
213 close F1;
|
|
|
214
|
|
|
215 print "Analytic correlation filtrering follow options done\n";
|
|
|
216
|
|
|
217
|
|
|
218 #############################################
|
|
|
219 ### Join groups that have been subdivided ###
|
|
|
220 #############################################
|
|
|
221
|
|
|
222 my @tdelete;
|
|
|
223
|
|
|
224 foreach my $k (keys %hcorrelgroup){
|
|
|
225 foreach my $i (keys %{$hcorrelgroup{$k}}){
|
|
|
226 foreach my $v (keys %hcorrelgroup){
|
|
|
227 my $count = 0;
|
|
|
228 if ($v ne $k){
|
|
|
229 foreach my $w (keys %{$hcorrelgroup{$v}}){
|
|
|
230 if($w eq $i){
|
|
|
231 $count = 1;
|
|
|
232 push(@tdelete, $v);
|
|
|
233 }
|
|
|
234 }
|
|
|
235 }
|
|
|
236 if($count == 1){
|
|
|
237 foreach my $w (keys %{$hcorrelgroup{$v}}){
|
|
|
238 $hcorrelgroup{$k}{$w}=$hcorrelgroup{$v}{$w};
|
|
|
239 }
|
|
|
240 delete($hcorrelgroup{$v});
|
|
|
241 }
|
|
|
242 }
|
|
|
243 }
|
|
|
244 }
|
|
|
245
|
|
|
246 foreach my $t (@tdelete){
|
|
|
247 delete($hcorrelgroup{$t});
|
|
|
248 }
|
|
|
249
|
|
|
250
|
|
|
251 ### Do it twice to see if it fix the problem of unmerge groups
|
|
|
252
|
|
|
253 foreach my $k (keys %hcorrelgroup){
|
|
|
254 foreach my $i (keys %{$hcorrelgroup{$k}}){
|
|
|
255 foreach my $v (keys %hcorrelgroup){
|
|
|
256 my $count = 0;
|
|
|
257 if ($v ne $k){
|
|
|
258 foreach my $w (keys %{$hcorrelgroup{$v}}){
|
|
|
259 if($w eq $i){
|
|
|
260 $count = 1;
|
|
|
261 push(@tdelete, $v);
|
|
|
262 }
|
|
|
263 }
|
|
|
264 }
|
|
|
265 if($count == 1){
|
|
|
266 foreach my $w (keys %{$hcorrelgroup{$v}}){
|
|
|
267 $hcorrelgroup{$k}{$w}=$hcorrelgroup{$v}{$w};
|
|
|
268 }
|
|
|
269 delete($hcorrelgroup{$v});
|
|
|
270 }
|
|
|
271 }
|
|
|
272 }
|
|
|
273 }
|
|
|
274
|
|
|
275 foreach my $t (@tdelete){
|
|
|
276 delete($hcorrelgroup{$t});
|
|
|
277 }
|
|
|
278
|
|
|
279 print "Join groups that have been subdivided done\n";
|
|
|
280
|
|
|
281 #######################################################
|
|
|
282 ### Addition of annotation information among groups ###
|
|
|
283 #######################################################
|
|
|
284
|
|
|
285 foreach my $k (keys %hcorrelgroup){
|
|
|
286 foreach my $i (keys %{$hcorrelgroup{$k}}){
|
|
|
287 foreach my $j (keys %{$hcorrelgroup{$k}}){
|
|
|
288 my $count = 0;
|
|
|
289 if ($i ne $j){
|
|
|
290
|
|
|
291 my $a = $hrtmz{$i}{mz};
|
|
|
292 my $b = $hrtmz{$j}{mz};
|
|
|
293
|
|
|
294 my $diff = $a - $b;
|
|
|
295 my $sign;
|
|
|
296 if($diff>0){
|
|
|
297 $sign="+";
|
|
|
298 }
|
|
|
299 if($diff<0){
|
|
|
300 $sign="-";
|
|
|
301 }
|
|
|
302 $diff = abs($diff);
|
|
|
303
|
|
|
304 foreach my $z (keys %{$refhmass}){
|
|
|
305
|
|
|
306 foreach my $y (keys %{$refhmass->{$z}}){
|
|
|
307 my $ym = $y - $mass_threshold;
|
|
|
308 my $yp = $y + $mass_threshold;
|
|
|
309
|
|
|
310
|
|
|
311 if(($diff <= $yp) && ($diff >= $ym)){
|
|
|
312 my $diff_list = $diff - $y;
|
|
|
313 $diff_list = abs($diff_list);
|
|
|
314 $diff_list = sprintf ("%0.6f", $diff_list);
|
|
|
315
|
|
|
316 if($hcorrelgroup{$k}{$i} eq 1){
|
|
|
317 my $val = "@".$j."|".$sign."(".$z.")(".$diff_list.")|";
|
|
|
318 $hcorrelgroup{$k}{$i}=$val;
|
|
|
319 $count ++;
|
|
|
320 }
|
|
|
321 else{
|
|
|
322 if($count == 0){
|
|
|
323 my $val = "@".$j."|".$sign."(".$z.")(".$diff_list.")|";
|
|
|
324 $hcorrelgroup{$k}{$i}.=$val;
|
|
|
325 $count ++;
|
|
|
326 }
|
|
|
327 else{
|
|
|
328 my $val = $sign."(".$z.")(".$diff_list.")|";
|
|
|
329 $hcorrelgroup{$k}{$i}.=$val;
|
|
|
330 $count ++;
|
|
|
331 }
|
|
|
332 }
|
|
|
333 }
|
|
|
334 }
|
|
|
335 }
|
|
|
336 }
|
|
|
337 }
|
|
|
338 }
|
|
|
339 }
|
|
|
340
|
|
|
341
|
|
|
342 print "Addition of annotation information among groups done\n";
|
|
|
343
|
|
|
344
|
|
|
345 ####################################################
|
|
|
346 ### Choose the representative ion for each group ###
|
|
|
347 ####################################################
|
|
|
348
|
|
|
349 my %hgrouprepres;
|
|
|
350
|
|
|
351 open(F3, $dataMatrix);
|
|
|
352
|
|
|
353 while (my $line = <F3>){
|
|
|
354 chomp $line;
|
|
|
355
|
|
|
356 my @tline = split (/\t/, $line);
|
|
|
357
|
|
|
358 foreach my $k (keys %hcorrelgroup){
|
|
|
359 foreach my $i (keys %{$hcorrelgroup{$k}}){
|
|
|
360 if($tline[0] eq $i){
|
|
|
361 $hgrouprepres{$k}{$i}{mass}=$hrtmz{$tline[0]}{mz};
|
|
|
362 my $intensity;
|
|
|
363 my $nbsubjects=0;
|
|
|
364 for(my $y=1;$y<scalar(@tline);$y++){
|
|
|
365 $intensity += $tline[$y];
|
|
|
366 $nbsubjects ++;
|
|
|
367 }
|
|
|
368 my $meanintensity = $intensity/$nbsubjects;
|
|
|
369 $hgrouprepres{$k}{$i}{intensity}=$meanintensity;
|
|
|
370 $hgrouprepres{$k}{$i}{squaredmassint}=($hgrouprepres{$k}{$i}{mass}**2)/($hgrouprepres{$k}{$i}{intensity});
|
|
|
371 }
|
|
|
372 }
|
|
|
373 }
|
|
|
374 }
|
|
|
375 close F3;
|
|
|
376
|
|
|
377 foreach my $z (keys %hgrouprepres){
|
|
|
378 my $max_intensity = 0;
|
|
|
379 my $max_int_ion = "";
|
|
|
380 my $max_mass = 0;
|
|
|
381 my $max_mass_ion = "";
|
|
|
382 my $max_squared = 0;
|
|
|
383 my $max_squared_ion = "";
|
|
|
384 foreach my $w (keys %{$hgrouprepres{$z}}){
|
|
|
385 if($hgrouprepres{$z}{$w}{intensity} > $max_intensity){
|
|
|
386 $max_intensity = $hgrouprepres{$z}{$w}{intensity};
|
|
|
387 $max_int_ion = $w;
|
|
|
388 }
|
|
|
389 if($hgrouprepres{$z}{$w}{mass} > $max_mass){
|
|
|
390 $max_mass = $hgrouprepres{$z}{$w}{mass};
|
|
|
391 $max_mass_ion = $w;
|
|
|
392 }
|
|
|
393 if($hgrouprepres{$z}{$w}{squaredmassint} > $max_squared){
|
|
|
394 $max_squared = $hgrouprepres{$z}{$w}{squaredmassint};
|
|
|
395 $max_squared_ion = $w;
|
|
|
396 }
|
|
|
397 }
|
|
|
398
|
|
|
399 my $max_int_max_mass_ion="";
|
|
|
400
|
|
|
401 if($repres_opt eq "max_intensity_max_mass"){
|
|
|
402 my %hfirst;
|
|
|
403 my $first=0;
|
|
|
404 foreach my $w (reverse sort {$hgrouprepres{$z}{$a}{intensity} <=> $hgrouprepres{$z}{$b}{intensity} } keys %{$hgrouprepres{$z}}){
|
|
|
405 $first ++;
|
|
|
406 if ($first <= 3){
|
|
|
407 $hfirst{$w} = $hgrouprepres{$z}{$w}{intensity};
|
|
|
408 }
|
|
|
409 }
|
|
|
410
|
|
|
411 my $first_2 = 0;
|
|
|
412 my $intens_max = 0;
|
|
|
413 my $mass_max = 0;
|
|
|
414
|
|
|
415 foreach my $y (reverse sort {$hfirst{$a} <=> $hfirst{$b}} keys %hfirst){
|
|
|
416
|
|
|
417 $first_2 ++;
|
|
|
418 if($first_2 == 1){
|
|
|
419 $intens_max = $hfirst{$y};
|
|
|
420 if($intensity_threshold > $intens_max){
|
|
|
421 $intensity_threshold = 0;
|
|
|
422 }
|
|
|
423 $max_int_max_mass_ion = $y;
|
|
|
424 $mass_max = $hgrouprepres{$z}{$y}{mass};
|
|
|
425 }
|
|
|
426 if($hgrouprepres{$z}{$y}{mass} > $mass_max){
|
|
|
427 if($hfirst{$y}>$intensity_threshold){
|
|
|
428 my $a = $intens_max * $intensity_pourc;
|
|
|
429 if($hfirst{$y} > $a){
|
|
|
430 $max_int_max_mass_ion = $y;
|
|
|
431 $mass_max = $hgrouprepres{$z}{$y}{mass};
|
|
|
432 }
|
|
|
433 }
|
|
|
434 }
|
|
|
435 }
|
|
|
436 }
|
|
|
437
|
|
|
438 $hgrouprepres{$z}{max_int}=$max_int_ion;
|
|
|
439 $hgrouprepres{$z}{max_mass}=$max_mass_ion;
|
|
|
440 $hgrouprepres{$z}{max_squared}=$max_squared_ion;
|
|
|
441 $hgrouprepres{$z}{max_int_max_mass}=$max_int_max_mass_ion;
|
|
|
442
|
|
|
443 }
|
|
|
444
|
|
|
445
|
|
|
446 print "Choose the representative ion for each group done\n";
|
|
|
447
|
|
|
448 #############################################################################
|
|
|
449 ### Addition of annotation information relative to the representative ion ###
|
|
|
450 #############################################################################
|
|
|
451
|
|
|
452 my %hreprescomparison;
|
|
|
453
|
|
|
454 my $representative="";
|
|
|
455
|
|
|
456 if($opt != 1){
|
|
|
457 foreach my $k (keys %hcorrelgroup){
|
|
|
458 foreach my $i (keys %{$hcorrelgroup{$k}}){
|
|
|
459
|
|
|
460 if($repres_opt eq "mass"){$representative = $hgrouprepres{$k}{max_mass}}
|
|
|
461 if($repres_opt eq "intensity"){$representative = $hgrouprepres{$k}{max_int}}
|
|
|
462 if($repres_opt eq "mixt"){$representative = $hgrouprepres{$k}{max_squared}}
|
|
|
463 if($repres_opt eq "max_intensity_max_mass"){$representative = $hgrouprepres{$k}{max_int_max_mass}}
|
|
|
464
|
|
|
465
|
|
|
466 my $count = 0;
|
|
|
467 if ($i ne $representative){
|
|
|
468
|
|
|
469 my $a = $hrtmz{$i}{mz};
|
|
|
470 my $b = $hrtmz{$representative}{mz};
|
|
|
471
|
|
|
472 my $diff = $a - $b;
|
|
|
473 my $sign;
|
|
|
474 if($diff>0){
|
|
|
475 $sign="+";
|
|
|
476 }
|
|
|
477 if($diff<0){
|
|
|
478 $sign="-";
|
|
|
479 }
|
|
|
480 $diff = abs($diff);
|
|
|
481
|
|
|
482 foreach my $z (keys %{$refhmass}){
|
|
|
483
|
|
|
484 foreach my $y (keys %{$refhmass->{$z}}){
|
|
|
485 my $ym = $y - $mass_threshold;
|
|
|
486 my $yp = $y + $mass_threshold;
|
|
|
487
|
|
|
488 if(($diff <= $yp) && ($diff >= $ym)){
|
|
|
489 my $diff_list = $diff - $y;
|
|
|
490 $diff_list = abs($diff_list);
|
|
|
491 $diff_list = sprintf ("%0.4f", $diff_list);
|
|
|
492 if($hcorrelgroup{$k}{$i} eq 1){
|
|
|
493 my $valrep = "[M ".$sign."(".$z.")]|";
|
|
|
494 $hreprescomparison{$k}{$i}{repres_diff}=$valrep;
|
|
|
495 $count ++;
|
|
|
496 }
|
|
|
497 else{
|
|
|
498 if($count == 0){
|
|
|
499 my $valrep = "[M ".$sign."(".$z.")]|";
|
|
|
500 $hreprescomparison{$k}{$i}{repres_diff}.=$valrep;
|
|
|
501 $count ++;
|
|
|
502 }
|
|
|
503 else{
|
|
|
504 my $valrep = "[M ".$sign."(".$z.")]|";
|
|
|
505 $hreprescomparison{$k}{$i}{repres_diff}.=$valrep;
|
|
|
506 $count ++;
|
|
|
507 }
|
|
|
508 }
|
|
|
509 }
|
|
|
510 }
|
|
|
511 }
|
|
|
512 }
|
|
|
513 else{
|
|
|
514 $hreprescomparison{$k}{$i}{repres_diff}="M";
|
|
|
515 }
|
|
|
516 }
|
|
|
517 }
|
|
|
518 }
|
|
|
519
|
|
|
520
|
|
|
521 print "Addition of annotation information relative to the representative ion done\n";
|
|
|
522
|
|
|
523 ##############################
|
|
|
524 ### Print in result file ! ###
|
|
|
525 ##############################
|
|
|
526
|
|
|
527 open(F4, ">$output_tabular");
|
|
|
528 open(F5, $combined_DMVM);
|
|
|
529
|
|
|
530 my $line_nb = 0;
|
|
|
531 my %hheader;
|
|
|
532 while (my $line = <F5>){
|
|
|
533 chomp $line;
|
|
|
534
|
|
|
535
|
|
|
536 my @tline = split (/\t/, $line);
|
|
|
537
|
|
|
538 if($line_nb == 0){
|
|
|
539 print F4 "$line\tACorF_groups";
|
|
|
540 if($opt == 1){
|
|
|
541 if($repres_opt eq "intensity"){print F4 "\tACorF_filter\tintensity_repres\n"}
|
|
|
542 if($repres_opt eq "mass"){print F4 "\tACorF_filter\tmass_repres\n"}
|
|
|
543 if($repres_opt eq "mixt"){print F4 "\tACorF_filter\tmass2intens_repres\n"}
|
|
|
544 if($repres_opt eq "max_intensity_max_mass"){print F4 "\tACorF_filter\tmax_intensity_max_mass_repres\n"}
|
|
|
545 }
|
|
|
546 else{
|
|
|
547 if($repres_opt eq "intensity"){print F4 "\tisotopes_adducts_fragments_[\@id|annotation(delta_annotation)]\tACorF_filter\tintensity_repres\tannotation_relative_to_representative\n"}
|
|
|
548 if($repres_opt eq "mass"){print F4 "\tisotopes_adducts_fragments_[\@id|annotation(delta_annotation)]\tACorF_filter\tmass_repres\tannotation_relative_to_representative\n"}
|
|
|
549 if($repres_opt eq "mixt"){print F4 "\tisotopes_adducts_fragments_[\@id|annotation(delta_annotation)]\tACorF_filter\tmass2intens_repres\tannotation_relative_to_representative\n"}
|
|
|
550 if($repres_opt eq "max_intensity_max_mass"){print F4 "\tisotopes_adducts_fragments_[\@id|annotation(delta_annotation)]\tACorF_filter\tmax_intensity_max_mass_repres\tannotation_relative_to_representative\n"}
|
|
|
551 }
|
|
|
552
|
|
|
553
|
|
|
554 ### Creation of a header hash
|
|
|
555 for(my $i=0; $i<scalar(@tline);$i++){
|
|
|
556 my $a = $tline[$i];
|
|
|
557 $hheader{$a}=$i;
|
|
|
558 }
|
|
|
559 }
|
|
|
560
|
|
|
561 else{
|
|
|
562 my $find = 0;
|
|
|
563 foreach my $v (keys %hcorrelgroup){
|
|
|
564 if(defined($hgrouprepres{$v}{$tline[0]})){
|
|
|
565 print F4 "$line\t$v";
|
|
|
566
|
|
|
567 if($opt != 1){
|
|
|
568 if(defined($hcorrelgroup{$v}{$tline[0]})){
|
|
|
569 print F4 "\t$hcorrelgroup{$v}{$tline[0]}\t";
|
|
|
570
|
|
|
571 }
|
|
|
572 else{
|
|
|
573 print F4 "\t";
|
|
|
574 }
|
|
|
575 }
|
|
|
576
|
|
|
577 if($repres_opt eq "intensity"){
|
|
|
578 if($tline[0] eq $hgrouprepres{$v}{max_int}){
|
|
|
579 print F4 "1\t";
|
|
|
580 }
|
|
|
581 else{
|
|
|
582 print F4 "0\t";
|
|
|
583 }
|
|
|
584 $find = 1;
|
|
|
585 }
|
|
|
586 if($repres_opt eq "mass"){
|
|
|
587 if($tline[0] eq $hgrouprepres{$v}{max_mass}){
|
|
|
588 print F4 "1\t";
|
|
|
589 }
|
|
|
590 else{
|
|
|
591 print F4 "0\t";
|
|
|
592 }
|
|
|
593 $find = 1;
|
|
|
594 }
|
|
|
595 if($repres_opt eq "mixt"){
|
|
|
596 if($tline[0] eq $hgrouprepres{$v}{max_squared}){
|
|
|
597 print F4 "1\t";
|
|
|
598 }
|
|
|
599 else{
|
|
|
600 print F4 "0\t";
|
|
|
601 }
|
|
|
602 $find = 1;
|
|
|
603 }
|
|
|
604 if($repres_opt eq "max_intensity_max_mass"){
|
|
|
605 if($tline[0] eq $hgrouprepres{$v}{max_int_max_mass}){
|
|
|
606 print F4 "1\t";
|
|
|
607 }
|
|
|
608 else{
|
|
|
609 print F4 "0\t";
|
|
|
610 }
|
|
|
611 $find = 1;
|
|
|
612 }
|
|
|
613
|
|
|
614 if($repres_opt eq "intensity"){print F4 "$hgrouprepres{$v}{max_int}\t"}
|
|
|
615 if($repres_opt eq "mass"){print F4 "$hgrouprepres{$v}{max_mass}\t"}
|
|
|
616 if($repres_opt eq "mixt"){print F4 "$hgrouprepres{$v}{max_squared}\t"}
|
|
|
617 if($repres_opt eq "max_intensity_max_mass"){print F4 "$hgrouprepres{$v}{max_int_max_mass}\t"}
|
|
|
618
|
|
|
619 if(defined($hreprescomparison{$v}{$tline[0]}{repres_diff})){
|
|
|
620 print F4 "$hreprescomparison{$v}{$tline[0]}{repres_diff}\n";
|
|
|
621 }
|
|
|
622 else{
|
|
|
623 print F4 "-\n";
|
|
|
624 }
|
|
|
625 }
|
|
|
626 }
|
|
|
627 if($find == 0){
|
|
|
628 $groupct ++;
|
|
|
629 my $group = "group".$groupct;
|
|
|
630 if($opt != 1){
|
|
|
631 print F4 "$line\t$group\t-\t-\t-\t-\n";
|
|
|
632 }
|
|
|
633 else{
|
|
|
634 print F4 "$line\t$group\t-\t-\n";
|
|
|
635 }
|
|
|
636 }
|
|
|
637 }
|
|
|
638 $line_nb ++;
|
|
|
639 }
|
|
|
640
|
|
|
641 print "Print in result file done\n";
|
|
|
642
|
|
|
643 print "All steps done\n";
|
|
|
644
|
