Mercurial > repos > brasset_jensen > srnapipe

annotate bin/ppp.pm @ 61:9185ca0a7b43 draft

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