annotate bin/ppp.pm @ 50:372ce358ef09 draft

Uploaded
author pierre.pouchin
date Thu, 06 Sep 2018 11:25:35 -0400
parents 4bc00caa60b4
children 9645d995fb3c
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
40
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
1 package ppp;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
2
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
3 use strict;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
4 use warnings;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
5 use FindBin;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
6 use lib $FindBin::Bin;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
7 use Rcall qw ( histogram );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
8 use Math::CDF;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
9
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
10 use Exporter;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
11 our @ISA = qw( Exporter );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
12 our @EXPORT_OK = qw( &ping_pong_partners );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
13
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
14 sub ping_pong_partners
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
15 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
16 my ( $TE_fai, $sam, $dir, $max ) = @_;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
17 my ( $hashRef, $dupRef, $hasPpp ) = count_mapped ( $TE_fai, $sam );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
18 my ( %num_per_overlap_size, $overlap_number, $reverseR, $begRev, $endRev, $sensR, $begSens, $endSens, $snum, $rnum, $overlap );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
19 my ( $SP, $AP, $SN, $AN, $txt );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
20 my $flag = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
21 my @distri_overlap = (); my @overlaps_names = ();
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
22
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
23 open my $ppp_f, '>', $dir."ppp.txt" || die "cannot create ppp.txt $!\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
24 foreach my $k ( sort keys %{$hashRef} )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
25 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
26 my $v = $hashRef->{$k};
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
27 my $TE_dir = $dir.$k.'/';
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
28
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
29 %num_per_overlap_size = (); $overlap_number = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
30 $flag = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
31 for ( my $i = 0; $i <= $#{$v->[1]} ; $i++ )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
32 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
33 $reverseR = ${$v->[1]}[$i] ;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
34 $begRev = $reverseR->[0];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
35 $endRev = $begRev + length($reverseR->[1]) - 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
36
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
37 my $revR = reverse($reverseR->[1]);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
38 $revR =~ tr/atgcuATGCU/tacgaTACGA/;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
39
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
40 for ( my $j = 0; $j <= $#{$v->[0]}; $j++ )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
41 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
42 $sensR = ${$v->[0]}[$j];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
43 $begSens = $sensR->[0];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
44 $endSens = $begSens + length($sensR->[1]) - 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
45
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
46 if ( $begSens <= $endRev && $endSens > $endRev )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
47 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
48 $flag = 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
49 mkdir $TE_dir;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
50 open $txt, '>', $TE_dir.'overlap_size.txt' || die "cannot open repartition\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
51
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
52 $overlap = $endRev - $begSens + 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
53 $snum = $dupRef->{$sensR->[0].$sensR->[1].$sensR->[2].$sensR->[3]};
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
54 $rnum = $dupRef->{$reverseR->[0].$reverseR->[1].$reverseR->[2].$reverseR->[3]};
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
55
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
56 if ( $overlap == 10 )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
57 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
58 $hasPpp->{ $sensR->[0].$sensR->[1].$sensR->[2].$sensR->[3] } = 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
59 $hasPpp->{ $reverseR->[0].$reverseR->[1].$reverseR->[2].$reverseR->[3] } = 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
60 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
61 next if $overlap > $max;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
62 if ( $snum < $rnum )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
63 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
64 $num_per_overlap_size{$overlap} += $snum;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
65 $overlap_number += $snum;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
66 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
67 else
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
68 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
69 $num_per_overlap_size{$overlap} += $rnum ;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
70 $overlap_number += $rnum ;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
71 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
72 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
73 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
74 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
75 if ( $max != 0 )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
76 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
77 my @overlaps = ();
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
78 push @overlaps_names, $k;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
79 for my $i (1..$max)
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
80 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
81 $num_per_overlap_size{$i} = 0 unless exists( $num_per_overlap_size{$i} );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
82 push @overlaps, $num_per_overlap_size{$i};
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
83 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
84 push @distri_overlap, \@overlaps;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
85 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
86
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
87 if ( $flag == 1 )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
88 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
89 open $AP, '>', $TE_dir."antisensPPP.txt" || die "cannot create antisensPPP\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
90 open $AN, '>', $TE_dir."antisens.txt" || die "cannot create antisens\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
91 for ( my $i = 0; $i <= $#{$v->[1]} ; $i++ )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
92 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
93 $reverseR = ${$v->[1]}[$i] ;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
94 my $revR = reverse($reverseR->[1]);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
95 $revR =~ tr/atgcuATGCU/tacgaTACGA/;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
96 $rnum = $dupRef->{$reverseR->[0].$reverseR->[1].$reverseR->[2].$reverseR->[3]};
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
97 if ( $hasPpp->{ $reverseR->[0].$reverseR->[1].$reverseR->[2].$reverseR->[3] } == 1 )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
98 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
99 print $AP ">$reverseR->[0]|$reverseR->[2]|$reverseR->[3]|$rnum\n$revR\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
100 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
101 else
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
102 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
103 print $AN ">$reverseR->[0]|$reverseR->[2]|$reverseR->[3]|$rnum\n$revR\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
104 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
105 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
106 close $AP; close $AN;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
107
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
108 open $SP, '>', $TE_dir."sensPPP.txt" || die "cannot create sensPPP\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
109 open $SN, '>', $TE_dir."sens.txt" || die "cannot create sens\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
110 for ( my $j = 0; $j <= $#{$v->[0]}; $j++ )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
111 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
112 $sensR = ${$v->[0]}[$j];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
113 $snum = $dupRef->{$sensR->[0].$sensR->[1].$sensR->[2].$sensR->[3]};
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
114 if ( $hasPpp->{ $sensR->[0].$sensR->[1].$sensR->[2].$sensR->[3] } == 1 )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
115 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
116 print $SP ">$sensR->[0]|$sensR->[2]|$sensR->[3]|$snum\n$sensR->[1]\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
117 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
118 else
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
119 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
120 print $SN ">$sensR->[0]|$sensR->[2]|$sensR->[3]|$snum\n$sensR->[1]\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
121 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
122 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
123 close $SP; close $SN;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
124
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
125 my $histo_png = $TE_dir.'histogram.png';
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
126 histogram( \%num_per_overlap_size, $histo_png, $overlap_number );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
127 print $txt "size\tnumber\tpercentage of the total overlap number\n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
128 foreach my $k ( sort {$a <=> $b} keys %num_per_overlap_size )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
129 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
130 my $percentage = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
131 $percentage = $num_per_overlap_size{$k} * 100 / $overlap_number unless $overlap_number == 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
132 print $txt "$k\t$num_per_overlap_size{$k}\t"; printf $txt "%.2f\n",$percentage;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
133 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
134 close $txt;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
135 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
136 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
137
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
138 foreach my $tabP ( @distri_overlap )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
139 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
140 my $sum = sum($tabP);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
141 my $ten = $tabP->[9];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
142 my $mean = mean($tabP);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
143 my $std = standard_deviation($tabP, $mean);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
144 my $zsc = z_significance($ten, $mean, $std);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
145 my $name = shift @overlaps_names;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
146 my $prob = 'NA';
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
147 $prob = 1 - &Math::CDF::pnorm( $zsc ) if $zsc ne 'NA';
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
148 print $ppp_f (join ("\t", $name, $sum, $ten, $mean, $std, $zsc, $prob ),"\n" );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
149 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
150 close $ppp_f;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
151 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
152
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
153 sub count_mapped
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
154 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
155 my ( $fai, $in_file ) = @_;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
156 my ( %mapped, %dup, %has_ppp );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
157
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
158 open my $f, '<', $fai || die "cannot open $fai $! \n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
159 while(<$f>)
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
160 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
161 if ($_ =~ /(.*)\t(\d+)\n/)
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
162 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
163 $mapped{$1} = [];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
164 $mapped{$1}->[0] = []; $mapped{$1}->[1] = [];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
165 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
166 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
167 close $f;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
168
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
169 open my $infile, "samtools view $in_file |"|| die "cannot open input file $! \n";
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
170 while(<$infile>)
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
171 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
172 unless ($_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
173 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
174 my @line = split (/\t/,$_);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
175 if ($line[1] == 0)
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
176 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
177 unless ( exists ($dup{$line[3].$line[9].$line[1].$line[2]}) )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
178 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
179 push @{$mapped{$line[2]}->[0]} , [$line[3], $line[9], $line[1], $line[2]];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
180 $has_ppp {$line[3].$line[9].$line[1].$line[2]} = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
181 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
182 $dup{$line[3].$line[9].$line[1].$line[2]}+=1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
183 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
184 elsif ($line[1] == 16)
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
185 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
186 unless ( exists ($dup{$line[3].$line[9].$line[1].$line[2]}) )
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
187 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
188 push @{$mapped{$line[2]}->[1]} , [$line[3], $line[9], $line[1], $line[2]];
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
189 $has_ppp{$line[3].$line[9].$line[1].$line[2]} = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
190 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
191 $dup{$line[3].$line[9].$line[1].$line[2]}+=1
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
192 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
193 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
194 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
195 close $infile;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
196 return (\%mapped, \%dup, \%has_ppp );
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
197 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
198
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
199 sub sum
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
200 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
201 my $arrayref = shift;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
202 my $result = 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
203 foreach (@$arrayref) {$result += $_}
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
204 return $result;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
205 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
206
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
207 sub mean
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
208 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
209 my $arrayref = shift;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
210 my $result;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
211 foreach (@$arrayref) {$result += $_}
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
212 return $result / scalar(@$arrayref);
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
213 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
214
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
215 sub standard_deviation
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
216 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
217 my ($arrayref, $mean) = @_;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
218 return sqrt ( mean ( [map $_**2 , @$arrayref ]) - ($mean**2));
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
219 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
220
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
221 sub z_significance
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
222 {
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
223 my ($ten, $mean, $std) = @_;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
224 my $z = 'NA';
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
225 $z = (($ten - $mean) / $std) if $std != 0;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
226 return $z;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
227 }
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
228
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
229 1;
4bc00caa60b4 Uploaded
pierre.pouchin
parents: 1
diff changeset
230