annotate rapsodyn/rapsosnp_stats4x.pl @ 2:761fecc07fa9 draft

Uploaded
author mcharles
date Thu, 11 Sep 2014 03:10:47 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
1 #!/usr/bin/perl
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
2 use strict;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
3 use warnings;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
4
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
5 my $read1_row = $ARGV[0];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
6 my $read2_row = $ARGV[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
7
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
8 my $read1_trimmed_part1 = $ARGV[2];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
9 my $read1_trimmed_part2 = $ARGV[3];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
10 my $read1_trimmed_part3 = $ARGV[4];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
11 my $read1_trimmed_part4 = $ARGV[5];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
12 my $read2_trimmed_part1 = $ARGV[6];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
13 my $read2_trimmed_part2 = $ARGV[7];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
14 my $read2_trimmed_part3 = $ARGV[8];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
15 my $read2_trimmed_part4 = $ARGV[9];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
16
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
17 my $sam_row_part1 = $ARGV[10];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
18 my $sam_row_part2 = $ARGV[11];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
19 my $sam_row_part3 = $ARGV[12];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
20 my $sam_row_part4 = $ARGV[13];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
21 my $sam_filtered_part1 = $ARGV[14];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
22 my $sam_filtered_part2 = $ARGV[15];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
23 my $sam_filtered_part3 = $ARGV[16];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
24 my $sam_filtered_part4 = $ARGV[17];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
25
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
26 my $mpileup_variant = $ARGV[18];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
27
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
28 my $list_filtered = $ARGV[19];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
29
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
30 my $blast_filtered_part1 = $ARGV[20];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
31 my $blast_filtered_part2 = $ARGV[21];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
32 my $blast_filtered_part3 = $ARGV[22];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
33 my $blast_filtered_part4 = $ARGV[23];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
34
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
35 my $snp_selected = $ARGV[24];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
36
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
37
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
38 open(INR1R, $read1_row) or die ("Can't open $read1_row\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
39 my $nbread=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
40 my $nbbase =0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
41 while (my $line1=<INR1R>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
42 my $line2 = <INR1R>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
43 my $line3 = <INR1R>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
44 my $line4 = <INR1R>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
45 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
46 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
47 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
48 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
49 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
50 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
51 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
52 print "Row Reads 1\t\tNumber of reads : ",$nbread,"\tnumber of bases : ",$nbbase,"\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
53 close (INR1R);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
54
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
55
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
56
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
57
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
58 open(INR2R, $read2_row) or die ("Can't open $read2_row\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
59 $nbread=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
60 $nbbase =0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
61 while (my $line1=<INR2R>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
62 my $line2 = <INR2R>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
63 my $line3 = <INR2R>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
64 my $line4 = <INR2R>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
65 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
66 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
67 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
68 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
69 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
70 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
71 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
72 print "Row Reads 2\t\tNumber of reads : ",$nbread,"\tnumber of bases : ",$nbbase,"\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
73 close (INR2R);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
74
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
75
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
76
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
77
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
78
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
79 open(INR1TP1, $read1_trimmed_part1) or die ("Can't open $read1_trimmed_part1\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
80 $nbread=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
81 $nbbase =0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
82 while (my $line1=<INR1TP1>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
83 my $line2 = <INR1TP1>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
84 my $line3 = <INR1TP1>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
85 my $line4 = <INR1TP1>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
86 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
87 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
88 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
89 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
90 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
91 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
92 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
93 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
94 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
95 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
96 close (INR1TP1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
97 open(INR1TP2, $read1_trimmed_part2) or die ("Can't open $read1_trimmed_part2\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
98 while (my $line1=<INR1TP2>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
99 my $line2 = <INR1TP2>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
100 my $line3 = <INR1TP2>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
101 my $line4 = <INR1TP2>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
102 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
103 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
104 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
105 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
106 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
107 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
108 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
109 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
110 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
111 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
112 close (INR1TP2);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
113 open(INR1TP3, $read1_trimmed_part3) or die ("Can't open $read1_trimmed_part3\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
114 while (my $line1=<INR1TP3>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
115 my $line2 = <INR1TP3>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
116 my $line3 = <INR1TP3>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
117 my $line4 = <INR1TP3>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
118 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
119 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
120 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
121 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
122 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
123 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
124 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
125 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
126 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
127 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
128 close (INR1TP3);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
129 open(INR1TP4, $read1_trimmed_part4) or die ("Can't open $read1_trimmed_part4\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
130 while (my $line1=<INR1TP4>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
131 my $line2 = <INR1TP4>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
132 my $line3 = <INR1TP4>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
133 my $line4 = <INR1TP4>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
134 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
135 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
136 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
137 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
138 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
139 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
140 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
141 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
142 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
143 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
144 close (INR1TP4);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
145 print "Trimmed Reads 1\t\tNumber of reads : ",$nbread,"\tnumber of bases : ",$nbbase,"\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
146
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
147
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
148
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
149
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
150 open(INR2TP1, $read2_trimmed_part1) or die ("Can't open $read2_trimmed_part1\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
151 $nbread=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
152 $nbbase =0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
153 while (my $line1=<INR2TP1>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
154 my $line2 = <INR2TP1>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
155 my $line3 = <INR2TP1>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
156 my $line4 = <INR2TP1>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
157 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
158 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
159 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
160 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
161 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
162 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
163 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
164 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
165 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
166 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
167 close (INR2TP2);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
168 open(INR2TP2, $read2_trimmed_part2) or die ("Can't open $read2_trimmed_part2\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
169 while (my $line1=<INR2TP2>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
170 my $line2 = <INR2TP2>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
171 my $line3 = <INR2TP2>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
172 my $line4 = <INR2TP2>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
173 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
174 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
175 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
176 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
177 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
178 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
179 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
180 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
181 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
182 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
183 close (INR2TP2);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
184 open(INR2TP3, $read2_trimmed_part3) or die ("Can't open $read2_trimmed_part3\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
185 while (my $line1=<INR2TP3>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
186 my $line2 = <INR2TP3>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
187 my $line3 = <INR2TP3>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
188 my $line4 = <INR2TP3>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
189 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
190 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
191 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
192 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
193 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
194 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
195 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
196 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
197 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
198 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
199 close (INR2TP3);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
200 open(INR2TP4, $read2_trimmed_part4) or die ("Can't open $read2_trimmed_part4\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
201 while (my $line1=<INR2TP4>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
202 my $line2 = <INR2TP4>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
203 my $line3 = <INR2TP4>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
204 my $line4 = <INR2TP4>;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
205 if ($line1 =~ /^@/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
206 $nbread++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
207 if ($line2=~/([ATGCNX]+)/i){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
208 $nbbase += length($1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
209 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
210 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
211 print STDERR "$line1\n$line2\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
212 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
213 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
214 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
215 close (INR2TP4);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
216 print "Trimmed Reads 2\t\tNumber of reads : ",$nbread,"\tnumber of bases : ",$nbbase,"\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
217
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
218
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
219
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
220
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
221 print "\nSAM row\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
222 open(SAMP1, $sam_row_part1) or die ("Can't open $sam_row_part1\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
223 my %bitscore;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
224 while (my $line=<SAMP1>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
225 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
226 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
227 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
228 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
229 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
230 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
231 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
232 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
233 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
234 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
235 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
236 close (SAMP1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
237 open(SAMP2, $sam_row_part2) or die ("Can't open $sam_row_part2\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
238 while (my $line=<SAMP2>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
239 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
240 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
241 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
242 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
243 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
244 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
245 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
246 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
247 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
248 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
249 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
250 close (SAMP2);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
251 open(SAMP3, $sam_row_part3) or die ("Can't open $sam_row_part3\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
252 while (my $line=<SAMP3>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
253 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
254 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
255 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
256 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
257 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
258 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
259 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
260 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
261 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
262 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
263 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
264 close (SAMP3);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
265 open(SAMP4, $sam_row_part4) or die ("Can't open $sam_row_part4\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
266 while (my $line=<SAMP4>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
267 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
268 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
269 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
270 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
271 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
272 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
273 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
274 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
275 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
276 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
277 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
278 close (SAMP4);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
279 print "bitscore\t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
280 foreach my $key (sort {$bitscore{$b} <=> $bitscore{$a}} keys %bitscore) {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
281 print $key,"\t*\t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
282 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
283 print "\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
284 print " number \t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
285 foreach my $key (sort {$bitscore{$b} <=> $bitscore{$a}} keys %bitscore) {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
286 print $bitscore{$key},"\t*\t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
287 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
288 print "\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
289
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
290
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
291
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
292
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
293 print "\nSAM filtered\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
294 open(SAMFP1, $sam_filtered_part1) or die ("Can't open $sam_filtered_part1\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
295 undef %bitscore;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
296 while (my $line=<SAMFP1>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
297 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
298 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
299 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
300 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
301 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
302 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
303 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
304 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
305 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
306 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
307 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
308 close (SAMFP1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
309 open(SAMFP2, $sam_filtered_part2) or die ("Can't open $sam_filtered_part2\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
310 while (my $line=<SAMFP2>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
311 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
312 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
313 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
314 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
315 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
316 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
317 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
318 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
319 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
320 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
321 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
322 close (SAMFP2);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
323 open(SAMFP3, $sam_filtered_part3) or die ("Can't open $sam_filtered_part3\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
324 while (my $line=<SAMFP3>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
325 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
326 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
327 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
328 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
329 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
330 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
331 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
332 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
333 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
334 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
335 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
336 close (SAMFP3);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
337 open(SAMFP4, $sam_filtered_part4) or die ("Can't open $sam_filtered_part4\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
338 while (my $line=<SAMFP4>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
339 if (($line !~ /^\@SQ/)&&($line !~ /^\@PG/)){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
340 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
341 my $bit = $fields[1];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
342 if ($bitscore{$bit}){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
343 $bitscore{$bit}++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
344 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
345 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
346 $bitscore{$bit}=1;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
347 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
348 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
349 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
350 close (SAMFP4);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
351 print "bitscore\t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
352 foreach my $key (sort {$bitscore{$b} <=> $bitscore{$a}} keys %bitscore) {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
353 print $key,"\t*\t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
354 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
355 print "\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
356 print " number \t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
357 foreach my $key (sort {$bitscore{$b} <=> $bitscore{$a}} keys %bitscore) {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
358 print $bitscore{$key},"\t*\t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
359 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
360 print "\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
361
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
362
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
363
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
364
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
365 print "\nMPILEUP variant\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
366 open(MPV, $mpileup_variant) or die ("Can't open $mpileup_variant\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
367 my $nbvariant=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
368 while (my $line=<MPV>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
369 my @fields = split(/\s+/,$line);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
370 if ($#fields >= 4){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
371 my $match = $fields[4];
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
372 $match =~ s/\$//g; #the read start at this position
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
373 $match =~ s/\^.//g; #the read end at this position followed by quality char
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
374 if ($match =~/[ACGTNacgtn]+/){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
375 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
376 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
377 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
378 else {
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
379 #print STDERR "Erreur : $line\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
380 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
381 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
382 print "Variant detected :\t$nbvariant\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
383 close (MPV);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
384
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
385
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
386
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
387
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
388
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
389
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
390 print "\nMPILEUP filtered without dubious position\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
391 open(LF, $list_filtered) or die ("Can't open $list_filtered\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
392 $nbvariant=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
393 while (my $line=<LF>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
394 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
395 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
396 print "Variant selected :\t$nbvariant\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
397 close (LF);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
398
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
399
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
400
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
401
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
402
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
403 print "\nMPILEUP filtered without dubious position and BLAST\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
404 open(BFP1, $blast_filtered_part1) or die ("Can't open $blast_filtered_part1\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
405 $nbvariant=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
406 while (my $line=<BFP1>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
407 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
408 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
409 close (BFP1);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
410 open(BFP2, $blast_filtered_part2) or die ("Can't open $blast_filtered_part2\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
411 while (my $line=<BFP2>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
412 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
413 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
414 close (BFP2);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
415 open(BFP3, $blast_filtered_part3) or die ("Can't open $blast_filtered_part3\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
416 while (my $line=<BFP3>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
417 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
418 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
419 close (BFP3);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
420 open(BFP4, $blast_filtered_part4) or die ("Can't open $blast_filtered_part4\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
421 while (my $line=<BFP4>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
422 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
423 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
424 close (BFP4);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
425 print "Variant selected :\t$nbvariant\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
426
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
427
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
428
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
429
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
430
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
431 print "\nSNP selected after mpileup filtering : \t";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
432 open(SNP, $snp_selected) or die ("Can't open $snp_selected\n");
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
433 $nbvariant=0;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
434 while (my $line=<SNP>){
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
435 $nbvariant++;
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
436 }
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
437
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
438 print "$nbvariant\n";
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
439 close (SNP);
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
440
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
441
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
442
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
443
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
444
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
445
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
446
761fecc07fa9 Uploaded
mcharles
parents:
diff changeset
447