|
17
|
1 #!/usr/bin/perl -w
|
|
|
2 #Filename:
|
|
|
3 #Author: Tian Dongmei
|
|
|
4 #Email: tiandm@big.ac.cn
|
|
|
5 #Date: 2013/7/19
|
|
|
6 #Modified:
|
|
|
7 #Description:
|
|
|
8 my $version=1.00;
|
|
|
9
|
|
|
10 use strict;
|
|
|
11 use Getopt::Long;
|
|
|
12 use RNA;
|
|
|
13
|
|
|
14 my %opts;
|
|
|
15 GetOptions(\%opts,"map=s","g=s","d:i","f:i","o=s","e:f","s=s","h");
|
|
|
16 if (!(defined $opts{map} and defined $opts{g} and defined $opts{o} and defined $opts{s} ) || defined $opts{h}) { #necessary arguments
|
|
|
17 &usage;
|
|
|
18 }
|
|
|
19
|
|
|
20 my $filein=$opts{'map'};
|
|
|
21 my $faout=$opts{'o'};
|
|
|
22 my $strout=$opts{'s'};
|
|
|
23 my $genome= $opts{'g'};
|
|
|
24
|
|
|
25 my $maxd=defined $opts{'d'} ? $opts{'d'} : 200;
|
|
|
26 my $flank=defined $opts{'f'}? $opts{'f'} : 10;
|
|
|
27
|
|
|
28 my $MAX_ENERGY=-18;
|
|
|
29 if (defined $opts{'e'}) {$MAX_ENERGY=$opts{'e'};}
|
|
|
30 my $MAX_UNPAIR=5;
|
|
|
31 my $MIN_PAIR=15;
|
|
|
32 my $MAX_SIZEDIFF=4;
|
|
|
33 my $MAX_BULGE=2;
|
|
|
34 my $ASYMMETRY=5;
|
|
|
35 my $MIN_UNPAIR=0;
|
|
|
36 my $MIN_SPACE=5;
|
|
|
37 my $MAX_SPACE=$maxd;
|
|
|
38 my $FLANK=$flank;
|
|
|
39
|
|
|
40 ######### load in genome sequences start ########
|
|
|
41 my %genome;
|
|
|
42 my %lng;
|
|
|
43 my $name;
|
|
|
44 open IN,"<$genome";
|
|
|
45 while (my $aline=<IN>) {
|
|
|
46 chomp $aline;
|
|
|
47 next if($aline=~/^\#/);
|
|
|
48 if ($aline=~/^>(\S+)/) {
|
|
|
49 $name=$1;
|
|
|
50 next;
|
|
|
51 }
|
|
|
52 $genome{$name} .=$aline;
|
|
|
53 }
|
|
|
54 close IN;
|
|
|
55 foreach my $key (keys %genome) {
|
|
|
56 $lng{$key}=length($genome{$key});
|
|
|
57 }
|
|
|
58 ####### load in genome sequences end ##########
|
|
|
59
|
|
|
60 my %breaks; ### reads number bigger than 3
|
|
|
61 open IN,"<$filein"; #input file
|
|
|
62 while (my $aline=<IN>) {
|
|
|
63 chomp $aline;
|
|
|
64 my @tmp=split/\t/,$aline;
|
|
|
65 $tmp[0]=~/_x(\d+)$/;
|
|
|
66 my $no=$1;
|
|
|
67 next if($no<3);
|
|
|
68 #my $trand=&find_strand($tmp[9]);
|
|
|
69 #my @pos=split/\.\./,$tmp[5];
|
|
|
70 my $end=$tmp[3]+length($tmp[4])-1;
|
|
|
71 if($tmp[1] eq "-"){$tmp[4]=revcom($tmp[4]);}
|
|
|
72 push @{$breaks{$tmp[2]}{$tmp[1]}},[$tmp[3],$end,$no,$tmp[4]]; ### 0 base
|
|
|
73 }
|
|
|
74 close IN;
|
|
|
75
|
|
|
76 my %cites; ### peaks
|
|
|
77 foreach my $chr (keys %breaks) {
|
|
|
78 foreach my $strand (keys %{$breaks{$chr}}) {
|
|
|
79 my @array=@{$breaks{$chr}{$strand}};
|
|
|
80 @array=sort{$a->[0]<=>$b->[0]} @array;
|
|
|
81 for (my $i=0;$i<@array;$i++) {
|
|
|
82 my $start=$array[$i][0];my $end=$array[$i][1];
|
|
|
83 my @subarray=();
|
|
|
84 push @subarray,$array[$i];
|
|
|
85
|
|
|
86 for (my $j=$i+1;$j<@array;$j++) {
|
|
|
87 if ($start<$array[$j][1] && $end>$array[$j][0]) {
|
|
|
88 push @subarray,$array[$j];
|
|
|
89 ($start,$end)=&newpos($start,$end,$array[$j][0],$array[$j][1]);
|
|
|
90 }
|
|
|
91 else{
|
|
|
92 $i=$j;
|
|
|
93 &find_cites(\@subarray,$chr,$strand);
|
|
|
94 last;
|
|
|
95 }
|
|
|
96 }
|
|
|
97 }
|
|
|
98 }
|
|
|
99 }
|
|
|
100
|
|
|
101 open FA,">$faout"; #output file
|
|
|
102 open STR,">$strout";
|
|
|
103 foreach my $chr (keys %cites) {
|
|
|
104 foreach my $strand (keys %{$cites{$chr}}) {
|
|
|
105
|
|
|
106 my @array2=@{$cites{$chr}{$strand}};
|
|
|
107 @array2=sort{$a->[0]<=>$b->[0]} @array2;
|
|
|
108 &excise(\@array2,$chr,$strand);
|
|
|
109 }
|
|
|
110 }
|
|
|
111 close FA;
|
|
|
112 close STR;
|
|
|
113 sub oneCiteDn{
|
|
|
114 my ($array,$a,$chr,$strand)=@_;
|
|
|
115
|
|
|
116 my $ss=$$array[$a][0]-$flank;
|
|
|
117 $ss=0 if($ss<0);
|
|
|
118 my $ee=$$array[$a][1]+$maxd+$flank;
|
|
|
119 $ee=$lng{$chr} if($ee>$lng{$chr});
|
|
|
120
|
|
|
121 my $seq=substr($genome{$chr},$ss,$ee-$ss+1);
|
|
|
122 if($strand eq "-"){$seq=revcom($seq);}
|
|
|
123
|
|
|
124 my $val=&ffw1($seq,$$array[$a][3],$chr,$strand,$ss,$ee);
|
|
|
125 return $val;
|
|
|
126 }
|
|
|
127 sub oneCiteUp{
|
|
|
128 my ($array,$a,$chr,$strand)=@_;
|
|
|
129
|
|
|
130 my $ss=$$array[$a][0]-$maxd-$flank;
|
|
|
131 $ss=0 if($ss<0);
|
|
|
132 my $ee=$$array[$a][1]+$flank;
|
|
|
133 $ee=$lng{$chr} if($ee>$lng{$chr});
|
|
|
134
|
|
|
135 my $seq=substr($genome{$chr},$ss,$ee-$ss+1);
|
|
|
136 if($strand eq "-"){$seq=revcom($seq);}
|
|
|
137
|
|
|
138 my $val=&ffw1($seq,$$array[$a][3],$chr,$strand,$ss,$ee);
|
|
|
139 return $val;
|
|
|
140
|
|
|
141 }
|
|
|
142
|
|
|
143 sub twoCites{
|
|
|
144 my ($array,$a,$b,$chr,$strand)=@_;
|
|
|
145
|
|
|
146 my $ss=$$array[$a][0]-$flank;
|
|
|
147 $ss=0 if($ss<0);
|
|
|
148 my $ee=$$array[$b][1]+$flank;
|
|
|
149 $ee=$lng{$chr} if($ee>$lng{$chr});
|
|
|
150
|
|
|
151 my $seq=substr($genome{$chr},$ss,$ee-$ss+1);
|
|
|
152 if($strand eq "-"){$seq=revcom($seq);}
|
|
|
153
|
|
|
154 # my( $str,$mfe)=RNA::fold($seq);
|
|
|
155 # return 0 if($mfe>$MAX_ENERGY); ### minimum mfe
|
|
|
156 my $val=&ffw2($seq,$$array[$a][3],$$array[$b][3],$chr,$strand,$ss,$ee);
|
|
|
157
|
|
|
158 return $val;
|
|
|
159
|
|
|
160 }
|
|
|
161 sub excise{
|
|
|
162 my ($cluster,$chr,$strand)=@_;
|
|
|
163
|
|
|
164 my $last_pos=0;
|
|
|
165 for (my $i=0;$i<@{$cluster};$i++) {
|
|
|
166 next if($$cluster[$i][0]<$last_pos);
|
|
|
167 my $ok=0;
|
|
|
168 for (my $j=$i+1;$j<@{$cluster} ;$j++) {
|
|
|
169 if($$cluster[$j][0]-$$cluster[$i][1]>$maxd){
|
|
|
170 $i=$j;
|
|
|
171 last;
|
|
|
172 }else{
|
|
|
173 $ok=&twoCites($cluster,$i,$j,$chr,$strand);
|
|
|
174 if($ok){ $last_pos=$$cluster[$j][1]+$flank; $i=$j; last;}
|
|
|
175 }
|
|
|
176 }
|
|
|
177 next if($ok);
|
|
|
178
|
|
|
179 $ok=&oneCiteDn($cluster,$i,$chr,$strand);
|
|
|
180 if($ok){$last_pos=$$cluster[$i][1]+$maxd+$flank; next;}
|
|
|
181 $ok=&oneCiteUp($cluster,$i,$chr,$strand);
|
|
|
182 if($ok){$last_pos=$$cluster[$i][1]+$flank;next;}
|
|
|
183 }
|
|
|
184
|
|
|
185
|
|
|
186 }
|
|
|
187
|
|
|
188 sub ffw2{
|
|
|
189 my ($seq,$tag1,$tag2,$chr,$strand,$ss,$ee)=@_;
|
|
|
190
|
|
|
191 my $N_count=$seq=~tr/N//; ## precursor sequence has not more than 5 Ns
|
|
|
192 if ($N_count > 5) {
|
|
|
193 return 0;
|
|
|
194 }
|
|
|
195
|
|
|
196 my $seq_length=length $seq;
|
|
|
197 # position tag1 and tag2
|
|
|
198 my $tag1_beg=index($seq,$tag1,0)+1;
|
|
|
199 if ($tag1_beg < 1) {
|
|
|
200 warn "[ffw2] coordinate error.\n";
|
|
|
201 # $fold->{reason}="coordinate error";
|
|
|
202 return 0;
|
|
|
203 }
|
|
|
204 my $tag2_beg=index($seq,$tag2,0)+1;
|
|
|
205 if ($tag2_beg < 1) {
|
|
|
206 warn "[ffw2] coordinate error.\n";
|
|
|
207 # $fold->{reason}="coordinate error";
|
|
|
208 return 0;
|
|
|
209 }
|
|
|
210 if ($tag2_beg < $tag1_beg) {
|
|
|
211 # swap tag1 and tag2
|
|
|
212 ($tag1,$tag2)=($tag2,$tag1);
|
|
|
213 ($tag1_beg,$tag2_beg)=($tag2_beg,$tag1_beg);
|
|
|
214 }
|
|
|
215 my $tag1_end=$tag1_beg+length($tag1)-1;
|
|
|
216 my $tag2_end=$tag2_beg+length($tag2)-1;
|
|
|
217 # re-clipping
|
|
|
218 my $beg=$tag1_beg-$FLANK; $beg=1 if $beg < 1;
|
|
|
219 my $end=$tag2_end+$FLANK; $end=$seq_length if $end > $seq_length;
|
|
|
220 $seq=substr($seq,$beg-1,$end-$beg+1);
|
|
|
221 $seq_length=length $seq;
|
|
|
222 # re-reposition
|
|
|
223 $tag1_beg=index($seq,$tag1,0)+1;
|
|
|
224 if ($tag1_beg < 1) {
|
|
|
225 warn "[ffw2] coordinate error.\n";
|
|
|
226 # $fold->{reason}="coordinate error";
|
|
|
227 return 0;
|
|
|
228 }
|
|
|
229
|
|
|
230 $tag2_beg=index($seq,$tag2,0)+1;
|
|
|
231 if ($tag2_beg < 1) {
|
|
|
232 warn "[ffw2] coordinate error.\n";
|
|
|
233 # $fold->{reason}="coordinate error";
|
|
|
234 return 0;
|
|
|
235 }
|
|
|
236 $tag1_end=$tag1_beg+length($tag1)-1;
|
|
|
237 $tag2_end=$tag2_beg+length($tag2)-1;
|
|
|
238
|
|
|
239 # fold
|
|
|
240 my ($struct,$mfe)=RNA::fold($seq);
|
|
|
241 $mfe=sprintf "%.2f", $mfe;
|
|
|
242 if ($mfe > $MAX_ENERGY) {return 0;}
|
|
|
243
|
|
|
244 # tag1
|
|
|
245 my $tag1_length=$tag1_end-$tag1_beg+1;
|
|
|
246 my $tag1_struct=substr($struct,$tag1_beg-1,$tag1_length);
|
|
|
247 my $tag1_arm=which_arm($tag1_struct);
|
|
|
248 my $tag1_unpair=$tag1_struct=~tr/.//;
|
|
|
249 my $tag1_pair=$tag1_length-$tag1_unpair;
|
|
|
250 my $tag1_max_bulge=biggest_bulge($tag1_struct);
|
|
|
251 if ($tag1_arm ne "5p") { return 0;} # tag not in stem
|
|
|
252 # if ($tag1_unpair > $MAX_UNPAIR) {$fold->{reason}="unpair=$tag1_unpair ($MAX_UNPAIR)"; return $pass}
|
|
|
253 if ($tag1_pair < $MIN_PAIR) {return 0;}
|
|
|
254 if ($tag1_max_bulge > $MAX_BULGE) {return 0;}
|
|
|
255
|
|
|
256 # tag2
|
|
|
257 my $tag2_length=$tag2_end-$tag2_beg+1;
|
|
|
258 my $tag2_struct=substr($struct,$tag2_beg-1,$tag2_length);
|
|
|
259 my $tag2_arm=which_arm($tag2_struct);
|
|
|
260 my $tag2_unpair=$tag2_struct=~tr/.//;
|
|
|
261 my $tag2_pair=$tag2_length-$tag2_unpair;
|
|
|
262 my $tag2_max_bulge=biggest_bulge($tag2_struct);
|
|
|
263 if ($tag2_arm ne "3p") {return 0;} # star not in stem
|
|
|
264 # if ($tag2_unpair > $MAX_UNPAIR) {$fold->{reason}="unpair=$tag2_unpair ($MAX_UNPAIR)"; return $pass}
|
|
|
265 if ($tag2_pair < $MIN_PAIR) {return 0;}
|
|
|
266 if ($tag2_max_bulge > $MAX_BULGE) {return 0;}
|
|
|
267
|
|
|
268 # space size between miR and miR*
|
|
|
269 my $space=$tag2_beg-$tag1_end-1;
|
|
|
270 if ($space < $MIN_SPACE) {return 0;}
|
|
|
271 if ($space > $MAX_SPACE) {return 0;}
|
|
|
272
|
|
|
273 # size diff of miR and miR*
|
|
|
274 my $size_diff=abs($tag1_length-$tag2_length);
|
|
|
275 if ($size_diff > $MAX_SIZEDIFF) {return 0;}
|
|
|
276
|
|
|
277 # build base pairing table
|
|
|
278 my %pairtable;
|
|
|
279 &parse_struct($struct,\%pairtable); # coords count from 1
|
|
|
280
|
|
|
281 my $asy1=get_asy(\%pairtable,$tag1_beg,$tag1_end);
|
|
|
282 my $asy2=get_asy(\%pairtable,$tag2_beg,$tag2_end);
|
|
|
283 my $asy=($asy1 < $asy2) ? $asy1 : $asy2;
|
|
|
284 if ($asy > $ASYMMETRY) {return 0}
|
|
|
285
|
|
|
286 # duplex fold, determine whether two matures like a miR/miR* ike duplex
|
|
|
287 my ($like_mir_duplex1,$duplex_pair,$overhang1,$overhang2)=likeMirDuplex1($tag1,$tag2);
|
|
|
288 # parse hairpin, determine whether two matures form miR/miR* duplex in hairpin context
|
|
|
289 my ($like_mir_duplex2,$duplex_pair2,$overhang_b,$overhang_t)=likeMirDuplex2(\%pairtable,$tag1_beg,$tag1_end,$tag2_beg,$tag2_end);
|
|
|
290 if ($like_mir_duplex1==0 && $like_mir_duplex2==0) {
|
|
|
291 return 0;
|
|
|
292 }
|
|
|
293
|
|
|
294 print FA ">$chr:$strand:$ss..$ee\n$seq\n";
|
|
|
295 print STR ">$chr:$strand:$ss..$ee\n$seq\n$struct\t($mfe)\n";
|
|
|
296
|
|
|
297 return 1;
|
|
|
298 }
|
|
|
299
|
|
|
300 sub ffw1{
|
|
|
301 my ($seq,$tag,$chr,$strand,$ss,$ee)=@_;
|
|
|
302 my $pass=0;
|
|
|
303
|
|
|
304 my $N_count=$seq=~tr/N//;
|
|
|
305 if ($N_count > 5) {
|
|
|
306 return 0;
|
|
|
307 }
|
|
|
308
|
|
|
309 my $seq_length=length $seq;
|
|
|
310 my $tag_length=length $tag;
|
|
|
311
|
|
|
312 # position
|
|
|
313 my $tag_beg=index($seq,$tag,0)+1;
|
|
|
314 if ($tag_beg < 1) {
|
|
|
315 warn "[ffw1] coordinate error.\n";
|
|
|
316 return $pass;
|
|
|
317 }
|
|
|
318 my $tag_end=$tag_beg+length($tag)-1;
|
|
|
319
|
|
|
320
|
|
|
321 # define candidate precursor by hybrid short arm to long arm, not solid enough
|
|
|
322 my($beg,$end)=define_precursor($seq,$tag);
|
|
|
323 if (not defined $beg) {
|
|
|
324 return $pass;
|
|
|
325 }
|
|
|
326 if (not defined $end) {
|
|
|
327 return $pass;
|
|
|
328 }
|
|
|
329 $seq=substr($seq,$beg-1,$end-$beg+1);
|
|
|
330 $seq_length=length $seq;
|
|
|
331
|
|
|
332
|
|
|
333 # fold
|
|
|
334 my ($struct,$mfe)=RNA::fold($seq);
|
|
|
335 $mfe=sprintf "%.2f",$mfe;
|
|
|
336 if ($mfe > $MAX_ENERGY) {
|
|
|
337 $pass=0;
|
|
|
338 return $pass;
|
|
|
339 }
|
|
|
340
|
|
|
341 # reposition
|
|
|
342 $tag_beg=index($seq,$tag,0)+1;
|
|
|
343 if ($tag_beg < 1) {
|
|
|
344 warn "[ffw1] coordinate error.\n";
|
|
|
345 return 0;
|
|
|
346 }
|
|
|
347 $tag_end=$tag_beg+length($tag)-1;
|
|
|
348
|
|
|
349 my $tag_struct=substr($struct,$tag_beg-1,$tag_length);
|
|
|
350 my $tag_arm=which_arm($tag_struct);
|
|
|
351 my $tag_unpair=$tag_struct=~tr/.//;
|
|
|
352 my $tag_pair=$tag_length-$tag_unpair;
|
|
|
353 my $tag_max_bulge=biggest_bulge($tag_struct);
|
|
|
354 if ($tag_arm eq "-") { return $pass;}
|
|
|
355 # if ($tag_unpair > $MAX_UNPAIR) {$fold->{reason}="unpair=$tag_unpair ($MAX_UNPAIR)"; return $pass}
|
|
|
356 if ($tag_pair < $MIN_PAIR) { return $pass;}
|
|
|
357 if ($tag_max_bulge > $MAX_BULGE) {return $pass;}
|
|
|
358
|
|
|
359 # build base pairing table
|
|
|
360 my %pairtable;
|
|
|
361 &parse_struct($struct,\%pairtable); # coords count from 1
|
|
|
362
|
|
|
363 # get star
|
|
|
364 my ($star_beg,$star_end)=get_star(\%pairtable,$tag_beg,$tag_end);
|
|
|
365 my $star=substr($seq,$star_beg-1,$star_end-$star_beg+1);
|
|
|
366 my $star_length=$star_end-$star_beg+1;
|
|
|
367 my $star_struct=substr($struct,$star_beg-1,$star_end-$star_beg+1);
|
|
|
368 my $star_arm=which_arm($star_struct);
|
|
|
369 my $star_unpair=$star_struct=~tr/.//;
|
|
|
370 my $star_pair=$star_length-$star_unpair;
|
|
|
371 my $star_max_bulge=biggest_bulge($star_struct);
|
|
|
372 if ($star_arm eq "-") { return $pass;}
|
|
|
373 # if ($star_unpair > $MAX_UNPAIR) {$fold->{reason}="unpair=$star_unpair ($MAX_UNPAIR)"; return $pass}
|
|
|
374 if ($star_pair < $MIN_PAIR) {return $pass;}
|
|
|
375 if ($star_max_bulge > $MAX_BULGE) {return $pass;}
|
|
|
376
|
|
|
377 if ($tag_arm eq $star_arm) {return $pass;}
|
|
|
378
|
|
|
379 # space size between miR and miR*
|
|
|
380 my $space;
|
|
|
381 if ($tag_beg < $star_beg) {
|
|
|
382 $space=$star_beg-$tag_end-1;
|
|
|
383 }
|
|
|
384 else {
|
|
|
385 $space=$tag_beg-$star_end-1;
|
|
|
386 }
|
|
|
387 if ($space < $MIN_SPACE) { return $pass;}
|
|
|
388 if ($space > $MAX_SPACE) { return $pass;}
|
|
|
389
|
|
|
390 # size diff
|
|
|
391 my $size_diff=abs($tag_length-$star_length);
|
|
|
392 if ($size_diff > $MAX_SIZEDIFF) { return $pass;}
|
|
|
393
|
|
|
394 # asymmetry
|
|
|
395 my $asy=get_asy(\%pairtable,$tag_beg,$tag_end);
|
|
|
396 if ($asy > $ASYMMETRY) {return $pass;}
|
|
|
397
|
|
|
398 $pass=1;
|
|
|
399 print FA ">$chr:$strand:$ss..$ee\n$seq\n";
|
|
|
400 print STR ">$chr:$strand:$ss..$ee\n$seq\n$struct\t($mfe)\n";
|
|
|
401 return $pass;
|
|
|
402
|
|
|
403 }
|
|
|
404 sub get_star {
|
|
|
405 my($table,$beg,$end)=@_;
|
|
|
406
|
|
|
407 my ($s1,$e1,$s2,$e2); # s1 pair to s2, e1 pair to e2
|
|
|
408 foreach my $i ($beg..$end) {
|
|
|
409 if (defined $table->{$i}) {
|
|
|
410 my $j=$table->{$i};
|
|
|
411 $s1=$i;
|
|
|
412 $s2=$j;
|
|
|
413 last;
|
|
|
414 }
|
|
|
415 }
|
|
|
416 foreach my $i (reverse ($beg..$end)) {
|
|
|
417 if (defined $table->{$i}) {
|
|
|
418 my $j=$table->{$i};
|
|
|
419 $e1=$i;
|
|
|
420 $e2=$j;
|
|
|
421 last;
|
|
|
422 }
|
|
|
423 }
|
|
|
424 # print "$s1,$e1 $s2,$e2\n";
|
|
|
425
|
|
|
426 # correct terminus
|
|
|
427 my $off1=$s1-$beg;
|
|
|
428 my $off2=$end-$e1;
|
|
|
429 $s2+=$off1;
|
|
|
430 $s2+=2; # 081009
|
|
|
431 $e2-=$off2; $e2=1 if $e2 < 1;
|
|
|
432 $e2+=2; $e2=1 if $e2 < 1; # 081009
|
|
|
433 ($s2,$e2)=($e2,$s2) if ($s2 > $e2);
|
|
|
434 return ($s2,$e2);
|
|
|
435 }
|
|
|
436
|
|
|
437 sub define_precursor {
|
|
|
438 my $seq=shift;
|
|
|
439 my $tag=shift;
|
|
|
440
|
|
|
441 my $seq_length=length $seq;
|
|
|
442 my $tag_length=length $tag;
|
|
|
443 my $tag_beg=index($seq,$tag,0)+1;
|
|
|
444 my $tag_end=$tag_beg+$tag_length-1;
|
|
|
445
|
|
|
446 # split the candidate region into short arm and long arm
|
|
|
447 my $tag_arm;
|
|
|
448 my ($larm,$larm_beg,$larm_end);
|
|
|
449 my ($sarm,$sarm_beg,$sarm_end);
|
|
|
450 if ($tag_beg-1 < $seq_length-$tag_end) { # on 5' arm
|
|
|
451 $sarm=substr($seq,0,$tag_end);
|
|
|
452 $larm=substr($seq,$tag_end);
|
|
|
453 $sarm_beg=1;
|
|
|
454 $sarm_end=$tag_end;
|
|
|
455 $larm_beg=$tag_end+1;
|
|
|
456 $larm_end=$seq_length;
|
|
|
457 $tag_arm="5p";
|
|
|
458 }
|
|
|
459 else {
|
|
|
460 $larm=substr($seq,0,$tag_beg-1); # on 3' arm
|
|
|
461 $sarm=substr($seq,$tag_beg-1);
|
|
|
462 $larm_beg=1;
|
|
|
463 $larm_end=$tag_beg-1;
|
|
|
464 $sarm_beg=$tag_beg;
|
|
|
465 $sarm_end=$seq_length;
|
|
|
466 $tag_arm="3p";
|
|
|
467 }
|
|
|
468
|
|
|
469 # print "$sarm_beg,$sarm_end $sarm\n";
|
|
|
470 # print "$larm_beg,$larm_end $larm\n";
|
|
|
471
|
|
|
472 # clipping short arm
|
|
|
473 if ($tag_arm eq "5p") {
|
|
|
474 $sarm_beg=$tag_beg-$flank; $sarm_beg=1 if $sarm_beg < 1;
|
|
|
475 $sarm=substr($seq,$sarm_beg-1,$sarm_end-$sarm_beg+1);
|
|
|
476 }
|
|
|
477 else {
|
|
|
478 $sarm_end=$tag_end+$flank; $sarm_end=$seq_length if $sarm_end > $seq_length;
|
|
|
479 $sarm=substr($seq,$sarm_beg-1,$sarm_end-$sarm_beg+1);
|
|
|
480 }
|
|
|
481 # print "$sarm_beg,$sarm_end $sarm\n";
|
|
|
482 # print "$larm_beg,$larm_end $larm\n";
|
|
|
483
|
|
|
484 # define the precursor by hybriding short arm to long arm
|
|
43
|
485 =cut #modify in 2014-10-28
|
|
17
|
486 my $duplex=RNA::duplexfold($sarm,$larm);
|
|
|
487 my $struct=$duplex->{structure};
|
|
|
488 my $energy=sprintf "%.2f", $duplex->{energy};
|
|
|
489 my ($str1,$str2)=split(/&/,$struct);
|
|
|
490 my $pair=$str1=~tr/(//;
|
|
|
491 # print "pair=$pair\n";
|
|
|
492 my $beg1=$duplex->{i}+1-length($str1);
|
|
|
493 my $end1=$duplex->{i};
|
|
|
494 my $beg2=$duplex->{j};
|
|
|
495 my $end2=$duplex->{j}+length($str2)-1;
|
|
43
|
496 =cut
|
|
|
497 ###### new codes begin
|
|
|
498 my $duplex=`perl -e 'print "$sarm\n$larm"' | RNAduplex`;
|
|
|
499 #(.(.(((.....(((.&))))))...).). 1,16 : 1,13 (-7.20)
|
|
|
500 my @tmpduplex=split/\s+/,$duplex;
|
|
|
501 my $struct=$tmpduplex[0];
|
|
|
502 $tmpduplex[-1]=~s/[(|)]//g;
|
|
|
503 my $energy=$tmpduplex[-1];
|
|
|
504 my ($str1,$str2)=split(/&/,$struct);
|
|
|
505 my $pair=$str1=~tr/(//;
|
|
|
506 my ($beg1,$end1)=split/,/,$tmpduplex[1];
|
|
|
507 my ($beg2,$end2)=split/,/,$tmpduplex[3];
|
|
|
508 ######## new codes end
|
|
|
509
|
|
17
|
510 # print "$beg1:$end1 $beg2:$end2\n";
|
|
|
511 # transform coordinates
|
|
|
512 $beg1=$beg1+$sarm_beg-1;
|
|
|
513 $end1=$end1+$sarm_beg-1;
|
|
|
514 $beg2=$beg2+$larm_beg-1;
|
|
|
515 $end2=$end2+$larm_beg-1;
|
|
|
516 # print "$beg1:$end1 $beg2:$end2\n";
|
|
|
517
|
|
|
518 my $off5p=$beg1-$sarm_beg;
|
|
|
519 my $off3p=$sarm_end-$end1;
|
|
|
520 $beg2-=$off3p; $beg2=1 if $beg2 < 1;
|
|
|
521 $end2+=$off5p; $end2=$seq_length if $end2 > $seq_length;
|
|
|
522
|
|
|
523 # print "$beg1:$end1 $beg2:$end2\n";
|
|
|
524
|
|
|
525 my $beg=$sarm_beg < $beg2 ? $sarm_beg : $beg2;
|
|
|
526 my $end=$sarm_end > $end2 ? $sarm_end : $end2;
|
|
|
527
|
|
|
528 return if $pair < $MIN_PAIR;
|
|
|
529 # print "$beg,$end\n";
|
|
|
530 return ($beg,$end);
|
|
|
531 }
|
|
|
532
|
|
|
533
|
|
|
534 # duplex fold, judge whether two short seqs like a miRNA/miRNA* duplex
|
|
|
535 sub likeMirDuplex1 {
|
|
|
536 my $seq1=shift;
|
|
|
537 my $seq2=shift;
|
|
|
538 my $like_mir_duplex=1;
|
|
|
539
|
|
|
540 my $length1=length $seq1;
|
|
|
541 my $length2=length $seq2;
|
|
43
|
542 =cut
|
|
17
|
543 my $duplex=RNA::duplexfold($seq1, $seq2);
|
|
|
544 my $duplex_struct=$duplex->{structure};
|
|
|
545 my $duplex_energy=sprintf "%.2f", $duplex->{energy};
|
|
|
546 my ($str1,$str2)=split(/&/,$duplex_struct);
|
|
|
547 my $beg1=$duplex->{i}+1-length($str1);
|
|
|
548 my $end1=$duplex->{i};
|
|
|
549 my $beg2=$duplex->{j};
|
|
|
550 my $end2=$duplex->{j}+length($str2)-1;
|
|
43
|
551 =cut
|
|
|
552 my $duplex=`perl -e 'print "$seq1\n$seq2"' | RNAduplex`;
|
|
|
553 #(.(.(((.....(((.&))))))...).). 1,16 : 1,13 (-7.20)
|
|
|
554 my @tmpduplex=split/\s+/,$duplex;
|
|
|
555 my $duplex_struct=$tmpduplex[0];
|
|
|
556 $tmpduplex[-1]=~s/[(|)]//g;
|
|
|
557 my $duplex_energy=$tmpduplex[-1];
|
|
|
558 my ($str1,$str2)=split(/&/,$duplex_struct);
|
|
|
559 #my $pair=$str1=~tr/(//;
|
|
|
560 my ($beg1,$end1)=split/,/,$tmpduplex[1];
|
|
|
561 my ($beg2,$end2)=split/,/,$tmpduplex[3];
|
|
|
562
|
|
17
|
563 # revise beg1, end1, beg2, end2
|
|
|
564 $str1=~/^(\.*)/;
|
|
|
565 $beg1+=length($1);
|
|
|
566 $str1=~/(\.*)$/;
|
|
|
567 $end1-=length($1);
|
|
|
568 $str2=~/^(\.*)/;
|
|
|
569 $beg2+=length($1);
|
|
|
570 $str2=~/(\.*)$/;
|
|
|
571 $end2-=length($1);
|
|
|
572
|
|
|
573 my $pair_num=$str1=~tr/(//;
|
|
|
574 my $overhang1=($length2-$end2)-($beg1-1); # 3' overhang at hairpin bottom
|
|
|
575 my $overhang2=($length1-$end1)-($beg2-1); # 3' overhang at hairpin neck
|
|
|
576 # print $pair_num,"\n";
|
|
|
577 # print $overhang1,"\n";
|
|
|
578 # print $overhang2,"\n";
|
|
|
579 if ($pair_num < 13) {
|
|
|
580 $like_mir_duplex=0;
|
|
|
581 }
|
|
|
582 if ($overhang1 < 0 || $overhang2 < 0 ) {
|
|
|
583 $like_mir_duplex=0;
|
|
|
584 }
|
|
|
585 if ($overhang1 > 4 || $overhang2 > 4) {
|
|
|
586 $like_mir_duplex=0;
|
|
|
587 }
|
|
|
588 return ($like_mir_duplex,$pair_num,$overhang1,$overhang2);
|
|
|
589 }
|
|
|
590
|
|
|
591 # judge whether two matures form miR/miR* duplex, in hairpin context
|
|
|
592 sub likeMirDuplex2 {
|
|
|
593 my ($table,$beg1,$end1,$beg2,$end2)=@_;
|
|
|
594 my $like_mir_duplex=1;
|
|
|
595
|
|
|
596 # s1 e1
|
|
|
597 # 5 ----------------------------3
|
|
|
598 # | | |||| ||| |
|
|
|
599 #3 -------------------------------5
|
|
|
600 # e2 s2
|
|
|
601
|
|
|
602 my $pair_num=0;
|
|
|
603 my $overhang1=0;
|
|
|
604 my $overhang2=0;
|
|
|
605 my ($s1,$e1,$s2,$e2);
|
|
|
606 foreach my $i ($beg1..$end1) {
|
|
|
607 if (defined $table->{$i}) {
|
|
|
608 my $j=$table->{$i};
|
|
|
609 if ($j <= $end2 && $j >= $beg2) {
|
|
|
610 $s1=$i;
|
|
|
611 $e2=$j;
|
|
|
612 last;
|
|
|
613 }
|
|
|
614 }
|
|
|
615 }
|
|
|
616 foreach my $i (reverse ($beg1..$end1)) {
|
|
|
617 if (defined $table->{$i}) {
|
|
|
618 my $j=$table->{$i};
|
|
|
619 if ($j <= $end2 && $j >= $beg2) {
|
|
|
620 $e1=$i;
|
|
|
621 $s2=$j;
|
|
|
622 last;
|
|
|
623 }
|
|
|
624 }
|
|
|
625 }
|
|
|
626
|
|
|
627 # print "$beg1,$end1 $s1,$e1\n";
|
|
|
628 # print "$beg2,$end2 $s2,$e2\n";
|
|
|
629
|
|
|
630 foreach my $i ($beg1..$end1) {
|
|
|
631 if (defined $table->{$i}) {
|
|
|
632 my $j=$table->{$i};
|
|
|
633 if ($j <= $end2 && $j >= $beg2) {
|
|
|
634 ++$pair_num;
|
|
|
635 }
|
|
|
636 }
|
|
|
637 }
|
|
|
638 if (defined $s1 && defined $e2) {
|
|
|
639 $overhang1=($end2-$e2)-($s1-$beg1);
|
|
|
640 }
|
|
|
641 if (defined $e1 && defined $s2) {
|
|
|
642 $overhang2=($end1-$e1)-($s2-$beg2);
|
|
|
643 }
|
|
|
644
|
|
|
645 if ($pair_num < 13) {
|
|
|
646 $like_mir_duplex=0;
|
|
|
647 }
|
|
|
648 if ($overhang1 < 0 && $overhang2 < 0) {
|
|
|
649 $like_mir_duplex=0;
|
|
|
650 }
|
|
|
651 return ($like_mir_duplex,$pair_num,$overhang1,$overhang2);
|
|
|
652 }
|
|
|
653 sub parse_struct {
|
|
|
654 my $struct=shift;
|
|
|
655 my $table=shift;
|
|
|
656
|
|
|
657 my @t=split('',$struct);
|
|
|
658 my @lbs; # left brackets
|
|
|
659 foreach my $k (0..$#t) {
|
|
|
660 if ($t[$k] eq "(") {
|
|
|
661 push @lbs, $k+1;
|
|
|
662 }
|
|
|
663 elsif ($t[$k] eq ")") {
|
|
|
664 my $lb=pop @lbs;
|
|
|
665 my $rb=$k+1;
|
|
|
666 $table->{$lb}=$rb;
|
|
|
667 $table->{$rb}=$lb;
|
|
|
668 }
|
|
|
669 }
|
|
|
670 if (@lbs) {
|
|
|
671 warn "unbalanced RNA struct.\n";
|
|
|
672 }
|
|
|
673 }
|
|
|
674 sub which_arm {
|
|
|
675 my $substruct=shift;
|
|
|
676 my $arm;
|
|
|
677 if ($substruct=~/\(/ && $substruct=~/\)/) {
|
|
|
678 $arm="-";
|
|
|
679 }
|
|
|
680 elsif ($substruct=~/\(/) {
|
|
|
681 $arm="5p";
|
|
|
682 }
|
|
|
683 else {
|
|
|
684 $arm="3p";
|
|
|
685 }
|
|
|
686 return $arm;
|
|
|
687 }
|
|
|
688 sub biggest_bulge {
|
|
|
689 my $struct=shift;
|
|
|
690 my $bulge_size=0;
|
|
|
691 my $max_bulge=0;
|
|
|
692 while ($struct=~/(\.+)/g) {
|
|
|
693 $bulge_size=length $1;
|
|
|
694 if ($bulge_size > $max_bulge) {
|
|
|
695 $max_bulge=$bulge_size;
|
|
|
696 }
|
|
|
697 }
|
|
|
698 return $max_bulge;
|
|
|
699 }
|
|
|
700 sub get_asy {
|
|
|
701 my($table,$a1,$a2)=@_;
|
|
|
702 my ($pre_i,$pre_j);
|
|
|
703 my $asymmetry=0;
|
|
|
704 foreach my $i ($a1..$a2) {
|
|
|
705 if (defined $table->{$i}) {
|
|
|
706 my $j=$table->{$i};
|
|
|
707 if (defined $pre_i && defined $pre_j) {
|
|
|
708 my $diff=($i-$pre_i)+($j-$pre_j);
|
|
|
709 $asymmetry += abs($diff);
|
|
|
710 }
|
|
|
711 $pre_i=$i;
|
|
|
712 $pre_j=$j;
|
|
|
713 }
|
|
|
714 }
|
|
|
715 return $asymmetry;
|
|
|
716 }
|
|
|
717
|
|
|
718 sub peaks{
|
|
|
719 my @cluster=@{$_[0]};
|
|
|
720
|
|
|
721 return if(@cluster<1);
|
|
|
722
|
|
|
723 my $max=0; my $index=-1;
|
|
|
724 for (my $i=0;$i<@cluster;$i++) {
|
|
|
725 if($cluster[$i][2]>$max){
|
|
|
726 $max=$cluster[$i][2];
|
|
|
727 $index=$i;
|
|
|
728 }
|
|
|
729 }
|
|
|
730 # &excise(\@cluster,$index,$_[1],$_[2]);
|
|
|
731 return($index);
|
|
|
732 }
|
|
|
733
|
|
|
734 sub find_cites{
|
|
|
735 my @tmp=@{$_[0]};
|
|
|
736 my $i=&peaks(\@tmp);
|
|
|
737
|
|
|
738 my $start=$tmp[$i][0];
|
|
|
739 my $total=0; my $node5=0;
|
|
|
740 for (my $j=0;$j<@tmp ;$j++) {
|
|
|
741 $total+=$tmp[$j][2];
|
|
|
742 $node5 +=$tmp[$j][2] if($tmp[$j][0]-$start<=2 && $tmp[$j][0]-$start>=-2);
|
|
|
743 }
|
|
|
744 push @{$cites{$_[1]}{$_[2]}},$tmp[$i] if($node5/$total>0.80 && $tmp[$i][2]/$node5>0.5);
|
|
|
745 }
|
|
|
746
|
|
|
747 sub newpos{
|
|
|
748 my ($a,$b,$c,$d)=@_;
|
|
|
749 my $s= $a>$c ? $c : $a;
|
|
|
750 my $e=$b>$d ? $b : $d;
|
|
|
751 return($s,$e);
|
|
|
752 }
|
|
|
753
|
|
|
754 sub rev{
|
|
|
755
|
|
|
756 my($sequence)=@_;
|
|
|
757
|
|
|
758 my $rev=reverse $sequence;
|
|
|
759
|
|
|
760 return $rev;
|
|
|
761 }
|
|
|
762
|
|
|
763 sub com{
|
|
|
764
|
|
|
765 my($sequence)=@_;
|
|
|
766
|
|
|
767 $sequence=~tr/acgtuACGTU/TGCAATGCAA/;
|
|
|
768
|
|
|
769 return $sequence;
|
|
|
770 }
|
|
|
771
|
|
|
772 sub revcom{
|
|
|
773
|
|
|
774 my($sequence)=@_;
|
|
|
775
|
|
|
776 my $revcom=rev(com($sequence));
|
|
|
777
|
|
|
778 return $revcom;
|
|
|
779 }
|
|
|
780
|
|
|
781 sub find_strand{
|
|
|
782
|
|
|
783 #A subroutine to find the strand, parsing different blast formats
|
|
|
784 my($other)=@_;
|
|
|
785
|
|
|
786 my $strand="+";
|
|
|
787
|
|
|
788 if($other=~/-/){
|
|
|
789 $strand="-";
|
|
|
790 }
|
|
|
791
|
|
|
792 if($other=~/minus/i){
|
|
|
793 $strand="-";
|
|
|
794 }
|
|
|
795
|
|
|
796 return($strand);
|
|
|
797 }
|
|
|
798 sub usage{
|
|
|
799 print <<"USAGE";
|
|
|
800 Version $version
|
|
|
801 Usage:
|
|
|
802 $0 -map -g -d -f -o -s -e
|
|
|
803 options:
|
|
|
804 -map input file# align result # bst. format
|
|
|
805 -g input file # genome sequence fasta format
|
|
|
806 -d <int> Maximal space between miRNA and miRNA* (200)
|
|
|
807 -f <int> Flank sequence length of miRNA precursor (10)
|
|
|
808 -o output file# percursor fasta file
|
|
|
809 -s output file# precursor structure file
|
|
|
810 -e <folat> Maximal free energy allowed for a miRNA precursor (-18 kcal/mol)
|
|
|
811
|
|
|
812 -h help
|
|
|
813 USAGE
|
|
|
814 exit(1);
|
|
|
815 }
|
|
|
816
|
