annotate rapsodyn/ParseBlastForUniqueMatch.pl @ 3:9332b9da7491 draft

Uploaded
author mcharles
date Thu, 11 Sep 2014 07:31:20 -0400
parents 442a7c88b886
children 3f7b0788a1c4
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
1 #!/usr/bin/perl
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
2 use strict;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
3 use warnings;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
4
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
5
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
6 my $input_variant_file = $ARGV[0];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
7 my $input_blast_file = $ARGV[1];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
8 my $window_length = $ARGV[2];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
9 my $nb_mismatch_max = $ARGV[3];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
10
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
11 my %hash_name;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
12
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
13 open(INB, $input_blast_file) or die ("Can't open $input_blast_file\n");
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
14 while (my $line =<INB>){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
15 my @fields = split (/\s+/,$line);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
16 # print $#fields,"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
17 # print $line;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
18 # exit(0);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
19 my $query_name;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
20 my $query_start;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
21 my $query_stop;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
22 my $query_aln;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
23 my $subject_aln;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
24 my $compt_mismatch_5p=0;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
25 my $compt_mismatch_3p=0;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
26
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
27 if ($#fields == 24){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
28 $query_name = $fields[0];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
29 $query_start = $fields[6];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
30 $query_stop = $fields[7];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
31 $query_aln = $fields[20];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
32 $subject_aln = $fields[21];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
33 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
34 elsif ($#fields == 4){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
35 $query_name = $fields[0];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
36 $query_start = $fields[1];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
37 $query_stop = $fields[2];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
38 $query_aln = $fields[3];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
39 $subject_aln = $fields[4];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
40 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
41 else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
42 print STDERR "Unrecongnized tabular format for blast result\nScript works with 25 column or 5 custom column(qseqid,qstart,qend,ssseq,sseq)\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
43 exit(0);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
44 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
45
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
46
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
47 my @field_name = split (/\_/,$query_name);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
48 my $name = $field_name[0]."_".$field_name[1];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
49 if ($query_name =~ /random/){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
50 $name .= "_".$field_name[2];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
51 # print $name."\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
52
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
53 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
54 if (!$hash_name{$name}){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
55 $hash_name{$name}=0;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
56 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
57 elsif ($hash_name{$name}>1){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
58 next;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
59 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
60 # if ($query_name eq "chrUnn_random_8279117_1_M1_a"){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
61 # print "READ : $query_name\n$name\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
62 # print "HASH : ",$hash_name{$name},"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
63 # # exit(0);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
64 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
65
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
66
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
67
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
68 chomp($query_aln);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
69 chomp($subject_aln);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
70
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
71 my $nb_gap_query=0;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
72
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
73 if (length($query_aln) == length($subject_aln)){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
74 if (length($query_aln)<$window_length-$nb_mismatch_max){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
75 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
76 else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
77 my @q = split(//,$query_aln);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
78 my @s = split(//,$subject_aln);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
79 for (my $i=0;$i<=$#q;$i++){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
80 my $global_idx = $query_start-1+$i-$nb_gap_query;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
81 if ($q[$i] eq "-"){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
82 if ($global_idx < $window_length){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
83 $compt_mismatch_5p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
84 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
85 elsif ($global_idx > $window_length){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
86 $compt_mismatch_3p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
87 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
88 $nb_gap_query++; #On compte les gap dans la query pour les soustraire de l'index global
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
89 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
90 else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
91 if ($q[$i] ne $s[$i]){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
92 if ($global_idx < $window_length){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
93 $compt_mismatch_5p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
94 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
95 elsif ($global_idx > $window_length){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
96 $compt_mismatch_3p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
97 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
98 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
99 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
100 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
101 $compt_mismatch_5p += $query_start-1;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
102 $compt_mismatch_3p += $window_length *2 + 1 - $query_stop;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
103
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
104 # for (my $i=0;$i<$window_length;$i++){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
105 # if ($tbl_q_aln[$i] eq "#"){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
106 # $compt_mismatch_5p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
107 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
108 # elsif ($tbl_q_aln[$i] ne $tbl_s_aln[$i]){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
109 # $compt_mismatch_5p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
110 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
111 # else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
112 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
113
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
114 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
115 # for (my $i=$window_length+1;$i<=$window_length*2;$i++){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
116 # if ($tbl_q_aln[$i] eq "#"){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
117 # $compt_mismatch_3p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
118 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
119 # elsif ($tbl_q_aln[$i] ne $tbl_s_aln[$i]){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
120 # $compt_mismatch_3p++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
121 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
122 # else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
123 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
124 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
125 if (($compt_mismatch_5p <= $nb_mismatch_max)||($compt_mismatch_3p <= $nb_mismatch_max)){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
126 $hash_name{$name}++;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
127 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
128
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
129 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
130 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
131 else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
132 print STDERR "incompatible subject and query alignement length\n $query_aln\n$subject_aln\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
133 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
134
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
135 # if ($line=~/chrCnn_random_49828229/){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
136 # print $line;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
137 # print $query_aln,"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
138 # print $subject_aln,"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
139 # print $compt_mismatch_5p,"\t",$compt_mismatch_3p,"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
140 # print $hash_name{"chrCnn_random_49828229"},"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
141 # print "\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
142 # }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
143
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
144
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
145
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
146 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
147 # exit(0);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
148
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
149 close (INB);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
150
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
151 open(INV, $input_variant_file) or die ("Can't open $input_variant_file\n");
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
152
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
153 while (my $ligne = <INV>) {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
154
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
155 my @champs = split (/\s+/,$ligne);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
156 my $header = $champs[0]."_".$champs[1];
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
157
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
158 if ($hash_name{$header}){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
159 if ($hash_name{$header}==1){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
160 print $ligne;
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
161 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
162 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
163 else {
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
164 #print STDERR "No blast result for ",$header,"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
165 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
166
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
167
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
168 }
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
169
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
170 close(INV);
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
171
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
172
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
173 # foreach my $key (sort keys %hash_name){
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
174 # print $key,"\t",$hash_name{$key},"\n";
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
175
442a7c88b886 Uploaded
mcharles
parents:
diff changeset
176 # }