annotate precursors.pl @ 51:3202911efdae draft

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