|
12
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 use warnings;
|
|
|
4 use strict;
|
|
|
5 use Getopt::Std;
|
|
44
|
6 use RNA;
|
|
12
|
7
|
|
|
8
|
|
|
9 ################################# MIRDEEP #################################################
|
|
|
10
|
|
|
11 ################################## USAGE ##################################################
|
|
|
12
|
|
|
13
|
|
|
14 my $usage=
|
|
|
15 "$0 file_signature file_structure temp_out_directory
|
|
|
16
|
|
|
17 This is the core algorithm of miRDeep. It takes as input a file in blastparsed format with
|
|
|
18 information on the positions of reads aligned to potential precursor sequences (signature).
|
|
|
19 It also takes as input an RNAfold output file, giving information on the sequence, structure
|
|
|
20 and mimimum free energy of the potential precursor sequences.
|
|
|
21
|
|
|
22 Extra arguments can be given. -s specifies a fastafile containing the known mature miRNA
|
|
|
23 sequences that should be considered for conservation purposes. -t prints out the potential
|
|
|
24 precursor sequences that do _not_ exceed the cut-off (default prints out the sequences that
|
|
|
25 exceeds the cut-off). -u gives limited output, that is only the ids of the potential precursors
|
|
|
26 that exceed the cut-off. -v varies the cut-off. -x is a sensitive option for Sanger sequences
|
|
|
27 obtained through conventional cloning. -z consider the number of base pairings in the lower
|
|
|
28 stems (this option is not well tested).
|
|
|
29
|
|
|
30 -h print this usage
|
|
|
31 -s fasta file with known miRNAs
|
|
|
32 #-o temp directory ,maked befor running the program.
|
|
|
33 -t print filtered
|
|
|
34 -u limited output (only ids)
|
|
|
35 -v cut-off (default 1)
|
|
|
36 -x sensitive option for Sanger sequences
|
|
|
37 -y use Randfold
|
|
|
38 -z consider Drosha processing
|
|
|
39 ";
|
|
|
40
|
|
|
41
|
|
|
42
|
|
|
43
|
|
|
44
|
|
|
45 ############################################################################################
|
|
|
46
|
|
|
47 ################################### INPUT ##################################################
|
|
|
48
|
|
|
49
|
|
|
50 #signature file in blast_parsed format
|
|
|
51 my $file_blast_parsed=shift or die $usage;
|
|
|
52
|
|
|
53 #structure file outputted from RNAfold
|
|
|
54 my $file_struct=shift or die $usage;
|
|
|
55
|
|
|
56 my $tmpdir=shift or die $usage;
|
|
|
57 #options
|
|
|
58 my %options=();
|
|
|
59 getopts("hs:tuv:xyz",\%options);
|
|
|
60
|
|
|
61
|
|
|
62
|
|
|
63
|
|
|
64
|
|
|
65
|
|
|
66 #############################################################################################
|
|
|
67
|
|
|
68 ############################# GLOBAL VARIABLES ##############################################
|
|
|
69
|
|
|
70
|
|
|
71 #parameters
|
|
|
72 my $nucleus_lng=11;
|
|
|
73
|
|
|
74 my $score_star=3.9;
|
|
|
75 my $score_star_not=-1.3;
|
|
|
76 my $score_nucleus=7.63;
|
|
|
77 my $score_nucleus_not=-1.17;
|
|
|
78 my $score_randfold=1.37;
|
|
|
79 my $score_randfold_not=-3.624;
|
|
|
80 my $score_intercept=0.3;
|
|
|
81 my @scores_stem=(-3.1,-2.3,-2.2,-1.6,-1.5,0.1,0.6,0.8,0.9,0.9,0);
|
|
|
82 my $score_min=1;
|
|
|
83 if($options{v}){$score_min=$options{v};}
|
|
|
84 if($options{x}){$score_min=-5;}
|
|
|
85
|
|
|
86 my $e=2.718281828;
|
|
|
87
|
|
|
88 #hashes
|
|
|
89 my %hash_desc;
|
|
|
90 my %hash_seq;
|
|
|
91 my %hash_struct;
|
|
|
92 my %hash_mfe;
|
|
|
93 my %hash_nuclei;
|
|
|
94 my %hash_mirs;
|
|
|
95 my %hash_query;
|
|
|
96 my %hash_comp;
|
|
|
97 my %hash_bp;
|
|
|
98
|
|
|
99 #other variables
|
|
|
100 my $subject_old;
|
|
|
101 my $message_filter;
|
|
|
102 my $message_score;
|
|
|
103 my $lines;
|
|
|
104 my $out_of_bound;
|
|
|
105
|
|
|
106
|
|
|
107
|
|
|
108 ##############################################################################################
|
|
|
109
|
|
|
110 ################################ MAIN ######################################################
|
|
|
111
|
|
|
112
|
|
|
113 #print help if that option is used
|
|
|
114 if($options{h}){die $usage;}
|
|
|
115 unless ($tmpdir=~/\/$/) {$tmpdir .="/";}
|
|
|
116 if(!(-s $tmpdir)){mkdir $tmpdir;}
|
|
|
117 $tmpdir .="TMP_DIR/";
|
|
|
118 mkdir $tmpdir;
|
|
|
119
|
|
|
120 #parse structure file outputted from RNAfold
|
|
|
121 parse_file_struct($file_struct);
|
|
|
122
|
|
|
123 #if conservation is scored, the fasta file of known miRNA sequences is parsed
|
|
|
124 if($options{s}){create_hash_nuclei($options{s})};
|
|
|
125
|
|
|
126 #parse signature file in blast_parsed format and resolve each potential precursor
|
|
|
127 parse_file_blast_parsed($file_blast_parsed);
|
|
44
|
128 #`rm -rf $tmpdir`;
|
|
12
|
129 exit;
|
|
|
130
|
|
|
131
|
|
|
132
|
|
|
133
|
|
|
134 ##############################################################################################
|
|
|
135
|
|
|
136 ############################## SUBROUTINES ###################################################
|
|
|
137
|
|
|
138
|
|
|
139
|
|
|
140 sub parse_file_blast_parsed{
|
|
|
141
|
|
|
142 # read through the signature blastparsed file, fills up a hash with information on queries
|
|
|
143 # (deep sequences) mapping to the current subject (potential precursor) and resolve each
|
|
|
144 # potential precursor in turn
|
|
|
145
|
|
|
146 my $file_blast_parsed=shift;
|
|
|
147
|
|
|
148 open (FILE_BLAST_PARSED, "<$file_blast_parsed") or die "can not open $file_blast_parsed\n";
|
|
|
149 while (my $line=<FILE_BLAST_PARSED>){
|
|
|
150 if($line=~/^(\S+)\s+(\S+)\s+(\d+)\.+(\d+)\s+(\S+)\s+(\S+)\s+(\d+)\.+(\d+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(.+)$/){
|
|
|
151 my $query=$1;
|
|
|
152 my $query_lng=$2;
|
|
|
153 my $query_beg=$3;
|
|
|
154 my $query_end=$4;
|
|
|
155 my $subject=$5;
|
|
|
156 my $subject_lng=$6;
|
|
|
157 my $subject_beg=$7;
|
|
|
158 my $subject_end=$8;
|
|
|
159 my $e_value=$9;
|
|
|
160 my $pid=$10;
|
|
|
161 my $bitscore=$11;
|
|
|
162 my $other=$12;
|
|
|
163
|
|
|
164 #if the new line concerns a new subject (potential precursor) then the old subject must be resolved
|
|
|
165 if($subject_old and $subject_old ne $subject){
|
|
|
166 resolve_potential_precursor();
|
|
|
167 }
|
|
|
168
|
|
|
169 #resolve the strand
|
|
|
170 my $strand=find_strand($other);
|
|
|
171
|
|
|
172 #resolve the number of reads that the deep sequence represents
|
|
|
173 my $freq=find_freq($query);
|
|
|
174
|
|
|
175 #read information of the query (deep sequence) into hash
|
|
|
176 $hash_query{$query}{"subject_beg"}=$subject_beg;
|
|
|
177 $hash_query{$query}{"subject_end"}=$subject_end;
|
|
|
178 $hash_query{$query}{"strand"}=$strand;
|
|
|
179 $hash_query{$query}{"freq"}=$freq;
|
|
|
180
|
|
|
181 #save the signature information
|
|
|
182 $lines.=$line;
|
|
|
183
|
|
|
184 $subject_old=$subject;
|
|
|
185 }
|
|
|
186 }
|
|
|
187 resolve_potential_precursor();
|
|
|
188 }
|
|
|
189
|
|
|
190 sub resolve_potential_precursor{
|
|
|
191
|
|
|
192 # dissects the potential precursor in parts by filling hashes, and tests if it passes the
|
|
|
193 # initial filter and the scoring filter
|
|
|
194
|
|
|
195 # binary variable whether the potential precursor is still viable
|
|
|
196 my $ret=1;
|
|
|
197 #print STDERR ">$subject_old\n";
|
|
|
198
|
|
|
199 fill_structure();
|
|
|
200 #print STDERR "\%hash_bp",scalar keys %hash_bp,"\n";
|
|
|
201 fill_pri();
|
|
|
202 #print STDERR "\%hash_comp",scalar keys %hash_comp,"\n";
|
|
|
203
|
|
|
204 fill_mature();
|
|
|
205 #print STDERR "\%hash_comp",scalar keys %hash_comp,"\n";
|
|
|
206
|
|
|
207 fill_star();
|
|
|
208 #print STDERR "\%hash_comp",scalar keys %hash_comp,"\n";
|
|
|
209
|
|
|
210 fill_loop();
|
|
|
211 #print STDERR "\%hash_comp",scalar keys %hash_comp,"\n";
|
|
|
212
|
|
|
213 fill_lower_flanks();
|
|
|
214 #print STDERR "\%hash_comp",scalar keys %hash_comp,"\n";
|
|
|
215
|
|
|
216 # do_test_assemble();
|
|
|
217
|
|
|
218 # this is the actual classification
|
|
|
219 unless(pass_filtering_initial() and pass_threshold_score()){$ret=0;}
|
|
|
220
|
|
|
221 print_results($ret);
|
|
|
222
|
|
|
223 reset_variables();
|
|
|
224
|
|
|
225 return;
|
|
|
226
|
|
|
227 }
|
|
|
228
|
|
|
229
|
|
|
230
|
|
|
231 sub print_results{
|
|
|
232
|
|
|
233 my $ret=shift;
|
|
|
234
|
|
|
235 # print out if the precursor is accepted and accepted precursors should be printed out
|
|
|
236 # or if the potential precursor is discarded and discarded potential precursors should
|
|
|
237 # be printed out
|
|
|
238
|
|
|
239 if((!$options{t} and $ret) or ($options{t} and !$ret)){
|
|
|
240 #full output
|
|
|
241 unless($options{u}){
|
|
|
242 if($message_filter){print $message_filter;}
|
|
|
243 if($message_score){print $message_score;}
|
|
|
244 print_hash_comp();
|
|
|
245 print $lines,"\n\n";
|
|
|
246 return;
|
|
|
247 }
|
|
|
248 #limited output (only ids)
|
|
|
249 my $id=$hash_comp{"pri_id"};
|
|
|
250 print "$id\n";
|
|
|
251 }
|
|
|
252 }
|
|
|
253
|
|
|
254
|
|
|
255
|
|
|
256
|
|
|
257
|
|
|
258
|
|
|
259
|
|
|
260 sub pass_threshold_score{
|
|
|
261
|
|
|
262 # this is the scoring
|
|
|
263
|
|
|
264 #minimum free energy of the potential precursor
|
|
|
265 # my $score_mfe=score_mfe($hash_comp{"pri_mfe"});
|
|
|
266 my $score_mfe=score_mfe($hash_comp{"pri_mfe"},$hash_comp{"pri_end"});
|
|
|
267
|
|
|
268 #count of reads that map in accordance with Dicer processing
|
|
|
269 my $score_freq=score_freq($hash_comp{"freq"});
|
|
|
270 #print STDERR "score_mfe: $score_mfe\nscore_freq: $score_freq\n";
|
|
|
271
|
|
|
272 #basic score
|
|
|
273 my $score=$score_mfe+$score_freq;
|
|
|
274
|
|
|
275 #scoring of conserved nucleus/seed (optional)
|
|
|
276 if($options{s}){
|
|
|
277
|
|
|
278 #if the nucleus is conserved
|
|
|
279 if(test_nucleus_conservation()){
|
|
|
280
|
|
|
281 #nucleus from position 2-8
|
|
|
282 my $nucleus=substr($hash_comp{"mature_seq"},1,$nucleus_lng);
|
|
|
283
|
|
|
284 #resolve DNA/RNA ambiguities
|
|
|
285 $nucleus=~tr/[T]/[U]/;
|
|
|
286
|
|
|
287 #print score contribution
|
|
|
288 score_s("score_nucleus\t$score_nucleus");
|
|
|
289
|
|
|
290 #print the ids of known miRNAs with same nucleus
|
|
|
291 score_s("$hash_mirs{$nucleus}");
|
|
|
292 #print STDERR "score_nucleus\t$score_nucleus\n";
|
|
|
293
|
|
|
294 #add to score
|
|
|
295 $score+=$score_nucleus;
|
|
|
296
|
|
|
297 #if the nucleus is not conserved
|
|
|
298 }else{
|
|
|
299 #print (negative) score contribution
|
|
|
300 score_s("score_nucleus\t$score_nucleus_not");
|
|
|
301
|
|
|
302 #add (negative) score contribution
|
|
|
303 $score+=$score_nucleus_not;
|
|
|
304 }
|
|
|
305 }
|
|
|
306
|
|
|
307 #if the majority of potential star reads fall as expected from Dicer processing
|
|
|
308 if($hash_comp{"star_read"}){
|
|
|
309 score_s("score_star\t$score_star");
|
|
|
310 #print STDERR "score_star\t$score_star\n";
|
|
|
311 $score+=$score_star;
|
|
|
312 }else{
|
|
|
313 score_s("score_star\t$score_star_not");
|
|
|
314 #print STDERR "score_star_not\t$score_star_not\n";
|
|
|
315 $score+=$score_star_not;
|
|
|
316 }
|
|
|
317
|
|
|
318 #score lower stems for potential for Drosha recognition (highly optional)
|
|
|
319 if($options{z}){
|
|
|
320 my $stem_bp=$hash_comp{"stem_bp"};
|
|
|
321 my $score_stem=$scores_stem[$stem_bp];
|
|
|
322 $score+=$score_stem;
|
|
|
323 score_s("score_stem\t$score_stem");
|
|
|
324 }
|
|
|
325
|
|
|
326 #print STDERR "score_intercept\t$score_intercept\n";
|
|
|
327
|
|
|
328 $score+=$score_intercept;
|
|
|
329
|
|
|
330 #score for randfold (optional)
|
|
|
331 if($options{y}){
|
|
|
332
|
|
|
333 # only calculate randfold value if it can make the difference between the potential precursor
|
|
|
334 # being accepted or discarded
|
|
|
335 if($score+$score_randfold>=$score_min and $score+$score_randfold_not<=$score_min){
|
|
|
336
|
|
|
337 #randfold value<0.05
|
|
|
338 if(test_randfold()){$score+=$score_randfold;score_s("score_randfold\t$score_randfold");}
|
|
|
339
|
|
|
340 #randfold value>0.05
|
|
|
341 else{$score+=$score_randfold_not;score_s("score_randfold\t$score_randfold_not");}
|
|
|
342 }
|
|
|
343 }
|
|
|
344
|
|
|
345 #round off values to one decimal
|
|
|
346 my $round_mfe=round($score_mfe*10)/10;
|
|
|
347 my $round_freq=round($score_freq*10)/10;
|
|
|
348 my $round=round($score*10)/10;
|
|
|
349
|
|
|
350 #print scores
|
|
|
351 score_s("score_mfe\t$round_mfe\nscore_freq\t$round_freq\nscore\t$round");
|
|
|
352
|
|
|
353 #return 1 if the potential precursor is accepted, return 0 if discarded
|
|
|
354 unless($score>=$score_min){return 0;}
|
|
|
355 return 1;
|
|
|
356 }
|
|
|
357
|
|
|
358 sub test_randfold{
|
|
|
359
|
|
|
360 #print sequence to temporary file, test randfold value, return 1 or 0
|
|
|
361
|
|
|
362 # print_file("pri_seq.fa",">pri_seq\n".$hash_comp{"pri_seq"});
|
|
44
|
363 #my $tmpfile=$tmpdir.$hash_comp{"pri_id"};
|
|
|
364 #open(FILE, ">$tmpfile");
|
|
|
365 #print FILE ">pri_seq\n",$hash_comp{"pri_seq"};
|
|
|
366 #close FILE;
|
|
12
|
367
|
|
|
368 # my $p_value=`randfold -s $tmpfile 999 | cut -f 3`;
|
|
44
|
369 #my $p1=`randfold -s $tmpfile 999 | cut -f 3`;
|
|
|
370 #my $p2=`randfold -s $tmpfile 999 | cut -f 3`;
|
|
|
371 my $p1=&randfold_pvalue($hash_comp{"pri_seq"},999);
|
|
|
372 my $p2=&randfold_pvalue($hash_comp{"pri_seq"},999);
|
|
12
|
373 my $p_value=($p1+$p2)/2;
|
|
|
374 wait;
|
|
|
375 # system "rm $tmpfile";
|
|
|
376
|
|
|
377 if($p_value<=0.05){return 1;}
|
|
|
378
|
|
|
379 return 0;
|
|
|
380 }
|
|
|
381
|
|
44
|
382 sub randfold_pvalue{
|
|
|
383 my $cpt_sup = 0;
|
|
|
384 my $cpt_inf = 0;
|
|
|
385 my $cpt_ega = 1;
|
|
|
386
|
|
|
387 my ($seq,$number_of_randomizations)=@_;
|
|
|
388 my $str =$seq;
|
|
|
389 my $mfe = RNA::fold($seq,$str);
|
|
12
|
390
|
|
44
|
391 for (my $i=0;$i<$number_of_randomizations;$i++) {
|
|
|
392 $seq = shuffle_sequence_dinucleotide($seq);
|
|
|
393 $str = $seq;
|
|
|
394
|
|
|
395 my $rand_mfe = RNA::fold($str,$str);
|
|
|
396
|
|
|
397 if ($rand_mfe < $mfe) {
|
|
|
398 $cpt_inf++;
|
|
|
399 }
|
|
|
400 if ($rand_mfe == $mfe) {
|
|
|
401 $cpt_ega++;
|
|
|
402 }
|
|
|
403 if ($rand_mfe > $mfe) {
|
|
|
404 $cpt_sup++;
|
|
|
405 }
|
|
|
406 }
|
|
|
407
|
|
|
408 my $proba = ($cpt_ega + $cpt_inf) / ($number_of_randomizations + 1);
|
|
12
|
409
|
|
44
|
410 #print "$name\t$mfe\t$proba\n";
|
|
|
411 return $proba;
|
|
|
412 }
|
|
|
413
|
|
|
414 sub shuffle_sequence_dinucleotide {
|
|
|
415
|
|
|
416 my ($str) = @_;
|
|
12
|
417
|
|
44
|
418 # upper case and convert to ATGC
|
|
|
419 $str = uc($str);
|
|
|
420 $str =~ s/U/T/g;
|
|
|
421
|
|
|
422 my @nuc = ('A','T','G','C');
|
|
|
423 my $count_swap = 0;
|
|
|
424 # set maximum number of permutations
|
|
|
425 my $stop = length($str) * 10;
|
|
|
426
|
|
|
427 while($count_swap < $stop) {
|
|
|
428
|
|
|
429 my @pos;
|
|
|
430
|
|
|
431 # look start and end letters
|
|
|
432 my $firstnuc = $nuc[int(rand 4)];
|
|
|
433 my $thirdnuc = $nuc[int(rand 4)];
|
|
12
|
434
|
|
44
|
435 # get positions for matching nucleotides
|
|
|
436 for (my $i=0;$i<(length($str)-2);$i++) {
|
|
|
437 if ((substr($str,$i,1) eq $firstnuc) && (substr($str,$i+2,1) eq $thirdnuc)) {
|
|
|
438 push (@pos,($i+1));
|
|
|
439 $i++;
|
|
|
440 }
|
|
|
441 }
|
|
|
442
|
|
|
443 # swap at random trinucleotides
|
|
|
444 my $max = scalar(@pos);
|
|
|
445 for (my $i=0;$i<$max;$i++) {
|
|
|
446 my $swap = int(rand($max));
|
|
|
447 if ((abs($pos[$swap] - $pos[$i]) >= 3) && (substr($str,$pos[$i],1) ne substr($str,$pos[$swap],1))) {
|
|
|
448 $count_swap++;
|
|
|
449 my $w1 = substr($str,$pos[$i],1);
|
|
|
450 my $w2 = substr($str,$pos[$swap],1);
|
|
|
451 substr($str,$pos[$i],1,$w2);
|
|
|
452 substr($str,$pos[$swap],1,$w1);
|
|
|
453 }
|
|
|
454 }
|
|
|
455 }
|
|
|
456 return($str);
|
|
|
457 }
|
|
12
|
458
|
|
|
459 sub test_nucleus_conservation{
|
|
|
460
|
|
|
461 #test if nucleus is identical to nucleus from known miRNA, return 1 or 0
|
|
|
462
|
|
|
463 my $nucleus=substr($hash_comp{"mature_seq"},1,$nucleus_lng);
|
|
|
464 $nucleus=~tr/[T]/[U]/;
|
|
|
465 if($hash_nuclei{$nucleus}){return 1;}
|
|
|
466
|
|
|
467 return 0;
|
|
|
468 }
|
|
|
469
|
|
|
470
|
|
|
471
|
|
|
472 sub pass_filtering_initial{
|
|
|
473
|
|
|
474 #test if the structure forms a plausible hairpin
|
|
|
475 unless(pass_filtering_structure()){filter_p("structure problem"); return 0;}
|
|
|
476
|
|
|
477 #test if >90% of reads map to the hairpin in consistence with Dicer processing
|
|
|
478 unless(pass_filtering_signature()){filter_p("signature problem");return 0;}
|
|
|
479
|
|
|
480 return 1;
|
|
|
481
|
|
|
482 }
|
|
|
483
|
|
|
484
|
|
|
485 sub pass_filtering_signature{
|
|
|
486
|
|
|
487 #number of reads that map in consistence with Dicer processing
|
|
|
488 my $consistent=0;
|
|
|
489
|
|
|
490 #number of reads that map inconsistent with Dicer processing
|
|
|
491 my $inconsistent=0;
|
|
|
492
|
|
|
493 # number of potential star reads map in good consistence with Drosha/Dicer processing
|
|
|
494 # (3' overhangs relative to mature product)
|
|
|
495 my $star_perfect=0;
|
|
|
496
|
|
|
497 # number of potential star reads that do not map in good consistence with 3' overhang
|
|
|
498 my $star_fuzzy=0;
|
|
|
499
|
|
|
500
|
|
|
501 #sort queries (deep sequences) by their position on the hairpin
|
|
|
502 my @queries=sort {$hash_query{$a}{"subject_beg"} <=> $hash_query{$b}{"subject_beg"}} keys %hash_query;
|
|
|
503
|
|
|
504 foreach my $query(@queries){
|
|
|
505
|
|
|
506 #number of reads that the deep sequence represents
|
|
|
507 unless(defined($hash_query{$query}{"freq"})){next;}
|
|
|
508 my $query_freq=$hash_query{$query}{"freq"};
|
|
|
509
|
|
|
510 #test which Dicer product (if any) the deep sequence corresponds to
|
|
|
511 my $product=test_query($query);
|
|
|
512
|
|
|
513 #if the deep sequence corresponds to a Dicer product, add to the 'consistent' variable
|
|
|
514 if($product){$consistent+=$query_freq;}
|
|
|
515
|
|
|
516 #if the deep sequence do not correspond to a Dicer product, add to the 'inconsistent' variable
|
|
|
517 else{$inconsistent+=$query_freq;}
|
|
|
518
|
|
|
519 #test a potential star sequence has good 3' overhang
|
|
|
520 if($product eq "star"){
|
|
|
521 if(test_star($query)){$star_perfect+=$query_freq;}
|
|
|
522 else{$star_fuzzy+=$query_freq;}
|
|
|
523 }
|
|
|
524 }
|
|
|
525
|
|
|
526 # if the majority of potential star sequences map in good accordance with 3' overhang
|
|
|
527 # score for the presence of star evidence
|
|
|
528 if($star_perfect>$star_fuzzy){$hash_comp{"star_read"}=1;}
|
|
|
529
|
|
|
530 #total number of reads mapping to the hairpin
|
|
|
531 my $freq=$consistent+$inconsistent;
|
|
|
532 $hash_comp{"freq"}=$freq;
|
|
|
533 unless($freq>0){filter_s("read frequency too low"); return 0;}
|
|
|
534
|
|
|
535 #unless >90% of the reads map in consistence with Dicer processing, the hairpin is discarded
|
|
|
536 my $inconsistent_fraction=$inconsistent/($inconsistent+$consistent);
|
|
|
537 unless($inconsistent_fraction<=0.1){filter_p("inconsistent\t$inconsistent\nconsistent\t$consistent"); return 0;}
|
|
|
538
|
|
|
539 #the hairpin is retained
|
|
|
540 return 1;
|
|
|
541 }
|
|
|
542
|
|
|
543 sub test_star{
|
|
|
544
|
|
|
545 #test if a deep sequence maps in good consistence with 3' overhang
|
|
|
546
|
|
|
547 my $query=shift;
|
|
|
548
|
|
|
549 #5' begin and 3' end positions
|
|
|
550 my $beg=$hash_query{$query}{"subject_beg"};
|
|
|
551 my $end=$hash_query{$query}{"subject_end"};
|
|
|
552
|
|
|
553 #the difference between observed and expected begin positions must be 0 or 1
|
|
|
554 my $offset=$beg-$hash_comp{"star_beg"};
|
|
|
555 if($offset==0 or $offset==1 or $offset==-1){return 1;}
|
|
|
556
|
|
|
557 return 0;
|
|
|
558 }
|
|
|
559
|
|
|
560
|
|
|
561
|
|
|
562 sub test_query{
|
|
|
563
|
|
|
564 #test if deep sequence maps in consistence with Dicer processing
|
|
|
565
|
|
|
566 my $query=shift;
|
|
|
567
|
|
|
568 #begin, end, strand and read count
|
|
|
569 my $beg=$hash_query{$query}{"subject_beg"};
|
|
|
570 my $end=$hash_query{$query}{"subject_end"};
|
|
|
571 my $strand=$hash_query{$query}{"strand"};
|
|
|
572 my $freq=$hash_query{$query}{"freq"};
|
|
|
573
|
|
|
574 #should not be on the minus strand (although this has in fact anecdotally been observed for known miRNAs)
|
|
|
575 if($strand eq '-'){return 0;}
|
|
|
576
|
|
|
577 #the deep sequence is allowed to stretch 2 nt beyond the expected 5' end
|
|
|
578 my $fuzz_beg=2;
|
|
|
579 #the deep sequence is allowed to stretch 5 nt beyond the expected 3' end
|
|
|
580 my $fuzz_end=2;
|
|
|
581
|
|
|
582 #if in accordance with Dicer processing, return the type of Dicer product
|
|
|
583 if(contained($beg,$end,$hash_comp{"mature_beg"}-$fuzz_beg,$hash_comp{"mature_end"}+$fuzz_end)){return "mature";}
|
|
|
584 if(contained($beg,$end,$hash_comp{"star_beg"}-$fuzz_beg,$hash_comp{"star_end"}+$fuzz_end)){return "star";}
|
|
|
585 if(contained($beg,$end,$hash_comp{"loop_beg"}-$fuzz_beg,$hash_comp{"loop_end"}+$fuzz_end)){return "loop";}
|
|
|
586
|
|
|
587 #if not in accordance, return 0
|
|
|
588 return 0;
|
|
|
589 }
|
|
|
590
|
|
|
591
|
|
|
592 sub pass_filtering_structure{
|
|
|
593
|
|
|
594 #The potential precursor must form a hairpin with miRNA precursor-like characteristics
|
|
|
595
|
|
|
596 #return value
|
|
|
597 my $ret=1;
|
|
|
598
|
|
|
599 #potential mature, star, loop and lower flank parts must be identifiable
|
|
|
600 unless(test_components()){return 0;}
|
|
|
601
|
|
|
602 #no bifurcations
|
|
|
603 unless(no_bifurcations_precursor()){$ret=0;}
|
|
|
604
|
|
|
605 #minimum 14 base pairings in duplex
|
|
|
606 unless(bp_duplex()>=15){$ret=0;filter_s("too few pairings in duplex");}
|
|
|
607
|
|
|
608 #not more than 6 nt difference between mature and star length
|
|
|
609 unless(-6<diff_lng() and diff_lng()<6){$ret=0; filter_s("too big difference between mature and star length") }
|
|
|
610
|
|
|
611 return $ret;
|
|
|
612 }
|
|
|
613
|
|
|
614
|
|
|
615
|
|
|
616
|
|
|
617
|
|
|
618
|
|
|
619 sub test_components{
|
|
|
620
|
|
|
621 #tests whether potential mature, star, loop and lower flank parts are identifiable
|
|
|
622
|
|
|
623 unless($hash_comp{"mature_struct"}){
|
|
|
624 filter_s("no mature");
|
|
|
625 # print STDERR "no mature\n";
|
|
|
626 return 0;
|
|
|
627 }
|
|
|
628
|
|
|
629 unless($hash_comp{"star_struct"}){
|
|
|
630 filter_s("no star");
|
|
|
631 # print STDERR "no star\n";
|
|
|
632 return 0;
|
|
|
633 }
|
|
|
634
|
|
|
635 unless($hash_comp{"loop_struct"}){
|
|
|
636 filter_s("no loop");
|
|
|
637 # print STDERR "no loop\n";
|
|
|
638 return 0;
|
|
|
639 }
|
|
|
640
|
|
|
641 unless($hash_comp{"flank_first_struct"}){
|
|
|
642 filter_s("no flanks");
|
|
|
643 #print STDERR "no flanks_first_struct\n";
|
|
|
644 return 0;
|
|
|
645 }
|
|
|
646
|
|
|
647 unless($hash_comp{"flank_second_struct"}){
|
|
|
648 filter_s("no flanks");
|
|
|
649 # print STDERR "no flanks_second_struct\n";
|
|
|
650 return 0;
|
|
|
651 }
|
|
|
652 return 1;
|
|
|
653 }
|
|
|
654
|
|
|
655
|
|
|
656
|
|
|
657
|
|
|
658
|
|
|
659 sub no_bifurcations_precursor{
|
|
|
660
|
|
|
661 #tests whether there are bifurcations in the hairpin
|
|
|
662
|
|
|
663 #assembles the potential precursor sequence and structure from the expected Dicer products
|
|
|
664 #this is the expected biological precursor, in contrast with 'pri_seq' that includes
|
|
|
665 #some genomic flanks on both sides
|
|
|
666
|
|
|
667 my $pre_struct;
|
|
|
668 my $pre_seq;
|
|
|
669 if($hash_comp{"mature_arm"} eq "first"){
|
|
|
670 $pre_struct.=$hash_comp{"mature_struct"}.$hash_comp{"loop_struct"}.$hash_comp{"star_struct"};
|
|
|
671 $pre_seq.=$hash_comp{"mature_seq"}.$hash_comp{"loop_seq"}.$hash_comp{"star_seq"};
|
|
|
672 }else{
|
|
|
673 $pre_struct.=$hash_comp{"star_struct"}.$hash_comp{"loop_struct"}.$hash_comp{"mature_struct"};
|
|
|
674 $pre_seq.=$hash_comp{"star_seq"}.$hash_comp{"loop_seq"}.$hash_comp{"mature_seq"};
|
|
|
675 }
|
|
|
676
|
|
|
677 #read into hash
|
|
|
678 $hash_comp{"pre_struct"}=$pre_struct;
|
|
|
679 $hash_comp{"pre_seq"}=$pre_seq;
|
|
|
680
|
|
|
681 #simple pattern matching checks for bifurcations
|
|
|
682 unless($pre_struct=~/^((\.|\()+..(\.|\))+)$/){
|
|
|
683 filter_s("bifurcation in precursor");
|
|
|
684 # print STDERR "bifurcation in precursor\n";
|
|
|
685 return 0;
|
|
|
686 }
|
|
|
687
|
|
|
688 return 1;
|
|
|
689 }
|
|
|
690
|
|
|
691 sub bp_precursor{
|
|
|
692
|
|
|
693 #total number of bps in the precursor
|
|
|
694
|
|
|
695 my $pre_struct=$hash_comp{"pre_struct"};
|
|
|
696
|
|
|
697 #simple pattern matching
|
|
|
698 my $pre_bps=0;
|
|
|
699 while($pre_struct=~/\(/g){
|
|
|
700 $pre_bps++;
|
|
|
701 }
|
|
|
702 return $pre_bps;
|
|
|
703 }
|
|
|
704
|
|
|
705
|
|
|
706 sub bp_duplex{
|
|
|
707
|
|
|
708 #total number of bps in the duplex
|
|
|
709
|
|
|
710 my $duplex_bps=0;
|
|
|
711 my $mature_struct=$hash_comp{"mature_struct"};
|
|
|
712
|
|
|
713 #simple pattern matching
|
|
|
714 while($mature_struct=~/(\(|\))/g){
|
|
|
715 $duplex_bps++;
|
|
|
716 }
|
|
|
717 return $duplex_bps;
|
|
|
718 }
|
|
|
719
|
|
|
720 sub diff_lng{
|
|
|
721
|
|
|
722 #find difference between mature and star lengths
|
|
|
723
|
|
|
724 my $mature_lng=length $hash_comp{"mature_struct"};
|
|
|
725 my $star_lng=length $hash_comp{"star_struct"};
|
|
|
726 my $diff_lng=$mature_lng-$star_lng;
|
|
|
727 return $diff_lng;
|
|
|
728 }
|
|
|
729
|
|
|
730
|
|
|
731
|
|
|
732 sub do_test_assemble{
|
|
|
733
|
|
|
734 # not currently used, tests if the 'pri_struct' as assembled from the parts (Dicer products, lower flanks)
|
|
|
735 # is identical to 'pri_struct' before disassembly into parts
|
|
|
736
|
|
|
737 my $assemble_struct;
|
|
|
738
|
|
|
739 if($hash_comp{"flank_first_struct"} and $hash_comp{"mature_struct"} and $hash_comp{"loop_struct"} and $hash_comp{"star_struct"} and $hash_comp{"flank_second_struct"}){
|
|
|
740 if($hash_comp{"mature_arm"} eq "first"){
|
|
|
741 $assemble_struct.=$hash_comp{"flank_first_struct"}.$hash_comp{"mature_struct"}.$hash_comp{"loop_struct"}.$hash_comp{"star_struct"}.$hash_comp{"flank_second_struct"};
|
|
|
742 }else{
|
|
|
743 $assemble_struct.=$hash_comp{"flank_first_struct"}.$hash_comp{"star_struct"}.$hash_comp{"loop_struct"}.$hash_comp{"mature_struct"}.$hash_comp{"flank_second_struct"};
|
|
|
744 }
|
|
|
745 unless($assemble_struct eq $hash_comp{"pri_struct"}){
|
|
|
746 $hash_comp{"test_assemble"}=$assemble_struct;
|
|
|
747 print_hash_comp();
|
|
|
748 }
|
|
|
749 }
|
|
|
750 return;
|
|
|
751 }
|
|
|
752
|
|
|
753
|
|
|
754
|
|
|
755 sub fill_structure{
|
|
|
756
|
|
|
757 #reads the dot bracket structure into the 'bp' hash where each key and value are basepaired
|
|
|
758
|
|
|
759 my $struct=$hash_struct{$subject_old};
|
|
|
760 my $lng=length $struct;
|
|
|
761
|
|
|
762 #local stack for keeping track of basepairings
|
|
|
763 my @bps;
|
|
|
764
|
|
|
765 for(my $pos=1;$pos<=$lng;$pos++){
|
|
|
766 my $struct_pos=excise_struct($struct,$pos,$pos,"+");
|
|
|
767
|
|
|
768 if($struct_pos eq "("){
|
|
|
769 push(@bps,$pos);
|
|
|
770 }
|
|
|
771
|
|
|
772 if($struct_pos eq ")"){
|
|
|
773 my $pos_prev=pop(@bps);
|
|
|
774 $hash_bp{$pos_prev}=$pos;
|
|
|
775 $hash_bp{$pos}=$pos_prev;
|
|
|
776 }
|
|
|
777 }
|
|
|
778 return;
|
|
|
779 }
|
|
|
780
|
|
|
781
|
|
|
782
|
|
|
783 sub fill_star{
|
|
|
784
|
|
|
785 #fills specifics on the expected star strand into 'comp' hash ('component' hash)
|
|
|
786
|
|
|
787 #if the mature sequence is not plausible, don't look for the star arm
|
|
|
788 my $mature_arm=$hash_comp{"mature_arm"};
|
|
|
789 unless($mature_arm){$hash_comp{"star_arm"}=0; return;}
|
|
|
790
|
|
|
791 #if the star sequence is not plausible, don't fill into the hash
|
|
|
792 my($star_beg,$star_end)=find_star();
|
|
|
793 my $star_arm=arm_star($star_beg,$star_end);
|
|
|
794 unless($star_arm){return;}
|
|
|
795
|
|
|
796 #excise expected star sequence and structure
|
|
|
797 my $star_seq=excise_seq($hash_comp{"pri_seq"},$star_beg,$star_end,"+");
|
|
|
798 my $star_struct=excise_seq($hash_comp{"pri_struct"},$star_beg,$star_end,"+");
|
|
|
799
|
|
|
800 #fill into hash
|
|
|
801 $hash_comp{"star_beg"}=$star_beg;
|
|
|
802 $hash_comp{"star_end"}=$star_end;
|
|
|
803 $hash_comp{"star_seq"}=$star_seq;
|
|
|
804 $hash_comp{"star_struct"}=$star_struct;
|
|
|
805 $hash_comp{"star_arm"}=$star_arm;
|
|
|
806
|
|
|
807 return;
|
|
|
808 }
|
|
|
809
|
|
|
810
|
|
|
811 sub find_star{
|
|
|
812
|
|
|
813 #uses the 'bp' hash to find the expected star begin and end positions from the mature positions
|
|
|
814
|
|
|
815 #the -2 is for the overhang
|
|
|
816 my $mature_beg=$hash_comp{"mature_beg"};
|
|
|
817 my $mature_end=$hash_comp{"mature_end"}-2;
|
|
|
818 my $mature_lng=$mature_end-$mature_beg+1;
|
|
|
819
|
|
|
820 #in some cases, the last nucleotide of the mature sequence does not form a base pair,
|
|
|
821 #and therefore does not basepair with the first nucleotide of the star sequence.
|
|
|
822 #In this case, the algorithm searches for the last nucleotide of the mature sequence
|
|
|
823 #to form a base pair. The offset is the number of nucleotides searched through.
|
|
|
824 my $offset_star_beg=0;
|
|
|
825 my $offset_beg=0;
|
|
|
826
|
|
|
827 #the offset should not be longer than the length of the mature sequence, then it
|
|
|
828 #means that the mature sequence does not form any base pairs
|
|
|
829 while(!$offset_star_beg and $offset_beg<$mature_lng){
|
|
|
830 if($hash_bp{$mature_end-$offset_beg}){
|
|
|
831 $offset_star_beg=$hash_bp{$mature_end-$offset_beg};
|
|
|
832 }else{
|
|
|
833 $offset_beg++;
|
|
|
834 }
|
|
|
835 }
|
|
|
836 #when defining the beginning of the star sequence, compensate for the offset
|
|
|
837 my $star_beg=$offset_star_beg-$offset_beg;
|
|
|
838
|
|
|
839 #same as above
|
|
|
840 my $offset_star_end=0;
|
|
|
841 my $offset_end=0;
|
|
|
842 while(!$offset_star_end and $offset_end<$mature_lng){
|
|
|
843 if($hash_bp{$mature_beg+$offset_end}){
|
|
|
844 $offset_star_end=$hash_bp{$mature_beg+$offset_end};
|
|
|
845 }else{
|
|
|
846 $offset_end++;
|
|
|
847 }
|
|
|
848 }
|
|
|
849 #the +2 is for the overhang
|
|
|
850 my $star_end=$offset_star_end+$offset_end+2;
|
|
|
851
|
|
|
852 return($star_beg,$star_end);
|
|
|
853 }
|
|
|
854
|
|
|
855
|
|
|
856 sub fill_pri{
|
|
|
857
|
|
|
858 #fills basic specifics on the precursor into the 'comp' hash
|
|
|
859
|
|
|
860 my $seq=$hash_seq{$subject_old};
|
|
|
861 my $struct=$hash_struct{$subject_old};
|
|
|
862 my $mfe=$hash_mfe{$subject_old};
|
|
|
863 my $length=length $seq;
|
|
|
864
|
|
|
865 $hash_comp{"pri_id"}=$subject_old;
|
|
|
866 $hash_comp{"pri_seq"}=$seq;
|
|
|
867 $hash_comp{"pri_struct"}=$struct;
|
|
|
868 $hash_comp{"pri_mfe"}=$mfe;
|
|
|
869 $hash_comp{"pri_beg"}=1;
|
|
|
870 $hash_comp{"pri_end"}=$length;
|
|
|
871
|
|
|
872 return;
|
|
|
873 }
|
|
|
874
|
|
|
875
|
|
|
876 sub fill_mature{
|
|
|
877
|
|
|
878 #fills specifics on the mature sequence into the 'comp' hash
|
|
|
879
|
|
|
880 my $mature_query=find_mature_query();
|
|
|
881 my($mature_beg,$mature_end)=find_positions_query($mature_query);
|
|
|
882 my $mature_strand=find_strand_query($mature_query);
|
|
|
883 my $mature_seq=excise_seq($hash_comp{"pri_seq"},$mature_beg,$mature_end,$mature_strand);
|
|
|
884 my $mature_struct=excise_struct($hash_comp{"pri_struct"},$mature_beg,$mature_end,$mature_strand);
|
|
|
885 my $mature_arm=arm_mature($mature_beg,$mature_end,$mature_strand);
|
|
|
886
|
|
|
887 $hash_comp{"mature_query"}=$mature_query;
|
|
|
888 $hash_comp{"mature_beg"}=$mature_beg;
|
|
|
889 $hash_comp{"mature_end"}=$mature_end;
|
|
|
890 $hash_comp{"mature_strand"}=$mature_strand;
|
|
|
891 $hash_comp{"mature_struct"}=$mature_struct;
|
|
|
892 $hash_comp{"mature_seq"}=$mature_seq;
|
|
|
893 $hash_comp{"mature_arm"}=$mature_arm;
|
|
|
894
|
|
|
895 return;
|
|
|
896 }
|
|
|
897
|
|
|
898
|
|
|
899
|
|
|
900 sub fill_loop{
|
|
|
901
|
|
|
902 #fills specifics on the loop sequence into the 'comp' hash
|
|
|
903
|
|
|
904 #unless both mature and star sequences are plausible, do not look for the loop
|
|
|
905 unless($hash_comp{"mature_arm"} and $hash_comp{"star_arm"}){return;}
|
|
|
906
|
|
|
907 my $loop_beg;
|
|
|
908 my $loop_end;
|
|
|
909
|
|
|
910 #defining the begin and end positions of the loop from the mature and star positions
|
|
|
911 #excision depends on whether the mature or star sequence is 5' of the loop ('first')
|
|
|
912 if($hash_comp{"mature_arm"} eq "first"){
|
|
|
913 $loop_beg=$hash_comp{"mature_end"}+1;
|
|
|
914 }else{
|
|
|
915 $loop_end=$hash_comp{"mature_beg"}-1;
|
|
|
916 }
|
|
|
917
|
|
|
918 if($hash_comp{"star_arm"} eq "first"){
|
|
|
919 $loop_beg=$hash_comp{"star_end"}+1;
|
|
|
920 }else{
|
|
|
921 $loop_end=$hash_comp{"star_beg"}-1;
|
|
|
922 }
|
|
|
923
|
|
|
924 #unless the positions are plausible, do not fill into hash
|
|
|
925 unless(test_loop($loop_beg,$loop_end)){return;}
|
|
|
926
|
|
|
927 my $loop_seq=excise_seq($hash_comp{"pri_seq"},$loop_beg,$loop_end,"+");
|
|
|
928 my $loop_struct=excise_struct($hash_comp{"pri_struct"},$loop_beg,$loop_end,"+");
|
|
|
929
|
|
|
930 $hash_comp{"loop_beg"}=$loop_beg;
|
|
|
931 $hash_comp{"loop_end"}=$loop_end;
|
|
|
932 $hash_comp{"loop_seq"}=$loop_seq;
|
|
|
933 $hash_comp{"loop_struct"}=$loop_struct;
|
|
|
934
|
|
|
935 return;
|
|
|
936 }
|
|
|
937
|
|
|
938
|
|
|
939 sub fill_lower_flanks{
|
|
|
940
|
|
|
941 #fills specifics on the lower flanks and unpaired strands into the 'comp' hash
|
|
|
942
|
|
|
943 #unless both mature and star sequences are plausible, do not look for the flanks
|
|
|
944 unless($hash_comp{"mature_arm"} and $hash_comp{"star_arm"}){return;}
|
|
|
945
|
|
|
946 my $flank_first_end;
|
|
|
947 my $flank_second_beg;
|
|
|
948
|
|
|
949 #defining the begin and end positions of the flanks from the mature and star positions
|
|
|
950 #excision depends on whether the mature or star sequence is 5' in the potenitial precursor ('first')
|
|
|
951 if($hash_comp{"mature_arm"} eq "first"){
|
|
|
952 $flank_first_end=$hash_comp{"mature_beg"}-1;
|
|
|
953 }else{
|
|
|
954 $flank_second_beg=$hash_comp{"mature_end"}+1;
|
|
|
955 }
|
|
|
956
|
|
|
957 if($hash_comp{"star_arm"} eq "first"){
|
|
|
958 $flank_first_end=$hash_comp{"star_beg"}-1;
|
|
|
959 }else{
|
|
|
960 $flank_second_beg=$hash_comp{"star_end"}+1;
|
|
|
961 }
|
|
|
962
|
|
|
963 #unless the positions are plausible, do not fill into hash
|
|
|
964 unless(test_flanks($flank_first_end,$flank_second_beg)){return;}
|
|
|
965
|
|
|
966 $hash_comp{"flank_first_end"}=$flank_first_end;
|
|
|
967 $hash_comp{"flank_second_beg"}=$flank_second_beg;
|
|
|
968 $hash_comp{"flank_first_seq"}=excise_seq($hash_comp{"pri_seq"},$hash_comp{"pri_beg"},$hash_comp{"flank_first_end"},"+");
|
|
|
969 $hash_comp{"flank_second_seq"}=excise_seq($hash_comp{"pri_seq"},$hash_comp{"flank_second_beg"},$hash_comp{"pri_end"},"+");
|
|
|
970 $hash_comp{"flank_first_struct"}=excise_struct($hash_comp{"pri_struct"},$hash_comp{"pri_beg"},$hash_comp{"flank_first_end"},"+");
|
|
|
971 $hash_comp{"flank_second_struct"}=excise_struct($hash_comp{"pri_struct"},$hash_comp{"flank_second_beg"},$hash_comp{"pri_end"},"+");
|
|
|
972
|
|
|
973 if($options{z}){
|
|
|
974 fill_stems_drosha();
|
|
|
975 }
|
|
|
976
|
|
|
977 return;
|
|
|
978 }
|
|
|
979
|
|
|
980
|
|
|
981 sub fill_stems_drosha{
|
|
|
982
|
|
|
983 #scores the number of base pairings formed by the first ten nt of the lower stems
|
|
|
984 #in general, the more stems, the higher the score contribution
|
|
|
985 #warning: this options has not been thoroughly tested
|
|
|
986
|
|
|
987 my $flank_first_struct=$hash_comp{"flank_first_struct"};
|
|
|
988 my $flank_second_struct=$hash_comp{"flank_second_struct"};
|
|
|
989
|
|
|
990 my $stem_first=substr($flank_first_struct,-10);
|
|
|
991 my $stem_second=substr($flank_second_struct,0,10);
|
|
|
992
|
|
|
993 my $stem_bp_first=0;
|
|
|
994 my $stem_bp_second=0;
|
|
|
995
|
|
|
996 #find base pairings by simple pattern matching
|
|
|
997 while($stem_first=~/\(/g){
|
|
|
998 $stem_bp_first++;
|
|
|
999 }
|
|
|
1000
|
|
|
1001 while($stem_second=~/\)/g){
|
|
|
1002 $stem_bp_second++;
|
|
|
1003 }
|
|
|
1004
|
|
|
1005 my $stem_bp=min2($stem_bp_first,$stem_bp_second);
|
|
|
1006
|
|
|
1007 $hash_comp{"stem_first"}=$stem_first;
|
|
|
1008 $hash_comp{"stem_second"}=$stem_second;
|
|
|
1009 $hash_comp{"stem_bp_first"}=$stem_bp_first;
|
|
|
1010 $hash_comp{"stem_bp_second"}=$stem_bp_second;
|
|
|
1011 $hash_comp{"stem_bp"}=$stem_bp;
|
|
|
1012
|
|
|
1013 return;
|
|
|
1014 }
|
|
|
1015
|
|
|
1016
|
|
|
1017
|
|
|
1018
|
|
|
1019 sub arm_mature{
|
|
|
1020
|
|
|
1021 #tests whether the mature sequence is in the 5' ('first') or 3' ('second') arm of the potential precursor
|
|
|
1022
|
|
|
1023 my ($beg,$end,$strand)=@_;
|
|
|
1024
|
|
|
1025 #mature and star sequences should alway be on plus strand
|
|
|
1026 if($strand eq "-"){return 0;}
|
|
|
1027
|
|
|
1028 #there should be no bifurcations and minimum one base pairing
|
|
|
1029 my $struct=excise_seq($hash_comp{"pri_struct"},$beg,$end,$strand);
|
|
|
1030 if(defined($struct) and $struct=~/^(\(|\.)+$/ and $struct=~/\(/){
|
|
|
1031 return "first";
|
|
|
1032 }elsif(defined($struct) and $struct=~/^(\)|\.)+$/ and $struct=~/\)/){
|
|
|
1033 return "second";
|
|
|
1034 }
|
|
|
1035 return 0;
|
|
|
1036 }
|
|
|
1037
|
|
|
1038
|
|
|
1039 sub arm_star{
|
|
|
1040
|
|
|
1041 #tests whether the star sequence is in the 5' ('first') or 3' ('second') arm of the potential precursor
|
|
|
1042
|
|
|
1043 my ($beg,$end)=@_;
|
|
|
1044
|
|
|
1045 #unless the begin and end positions are plausible, test negative
|
|
|
1046 unless($beg>0 and $beg<=$hash_comp{"pri_end"} and $end>0 and $end<=$hash_comp{"pri_end"} and $beg<=$end){return 0;}
|
|
|
1047
|
|
|
1048 #no overlap between the mature and the star sequence
|
|
|
1049 if($hash_comp{"mature_arm"} eq "first"){
|
|
|
1050 ($hash_comp{"mature_end"}<$beg) or return 0;
|
|
|
1051 }elsif($hash_comp{"mature_arm"} eq "second"){
|
|
|
1052 ($end<$hash_comp{"mature_beg"}) or return 0;
|
|
|
1053 }
|
|
|
1054
|
|
|
1055 #there should be no bifurcations and minimum one base pairing
|
|
|
1056 my $struct=excise_seq($hash_comp{"pri_struct"},$beg,$end,"+");
|
|
|
1057 if($struct=~/^(\(|\.)+$/ and $struct=~/\(/){
|
|
|
1058 return "first";
|
|
|
1059 }elsif($struct=~/^(\)|\.)+$/ and $struct=~/\)/){
|
|
|
1060 return "second";
|
|
|
1061 }
|
|
|
1062 return 0;
|
|
|
1063 }
|
|
|
1064
|
|
|
1065
|
|
|
1066 sub test_loop{
|
|
|
1067
|
|
|
1068 #tests the loop positions
|
|
|
1069
|
|
|
1070 my ($beg,$end)=@_;
|
|
|
1071
|
|
|
1072 #unless the begin and end positions are plausible, test negative
|
|
|
1073 unless($beg>0 and $beg<=$hash_comp{"pri_end"} and $end>0 and $end<=$hash_comp{"pri_end"} and $beg<=$end){return 0;}
|
|
|
1074
|
|
|
1075 return 1;
|
|
|
1076 }
|
|
|
1077
|
|
|
1078
|
|
|
1079 sub test_flanks{
|
|
|
1080
|
|
|
1081 #tests the positions of the lower flanks
|
|
|
1082
|
|
|
1083 my ($beg,$end)=@_;
|
|
|
1084
|
|
|
1085 #unless the begin and end positions are plausible, test negative
|
|
|
1086 unless($beg>0 and $beg<=$hash_comp{"pri_end"} and $end>0 and $end<=$hash_comp{"pri_end"} and $beg<=$end){return 0;}
|
|
|
1087
|
|
|
1088 return 1;
|
|
|
1089 }
|
|
|
1090
|
|
|
1091
|
|
|
1092 sub comp{
|
|
|
1093
|
|
|
1094 #subroutine to retrive from the 'comp' hash
|
|
|
1095
|
|
|
1096 my $type=shift;
|
|
|
1097 my $component=$hash_comp{$type};
|
|
|
1098 return $component;
|
|
|
1099 }
|
|
|
1100
|
|
|
1101
|
|
|
1102 sub find_strand_query{
|
|
|
1103
|
|
|
1104 #subroutine to find the strand for a given query
|
|
|
1105
|
|
|
1106 my $query=shift;
|
|
|
1107 my $strand=$hash_query{$query}{"strand"};
|
|
|
1108 return $strand;
|
|
|
1109 }
|
|
|
1110
|
|
|
1111
|
|
|
1112 sub find_positions_query{
|
|
|
1113
|
|
|
1114 #subroutine to find the begin and end positions for a given query
|
|
|
1115
|
|
|
1116 my $query=shift;
|
|
|
1117 my $beg=$hash_query{$query}{"subject_beg"};
|
|
|
1118 my $end=$hash_query{$query}{"subject_end"};
|
|
|
1119 return ($beg,$end);
|
|
|
1120 }
|
|
|
1121
|
|
|
1122
|
|
|
1123
|
|
|
1124 sub find_mature_query{
|
|
|
1125
|
|
|
1126 #finds the query with the highest frequency of reads and returns it
|
|
|
1127 #is used to determine the positions of the potential mature sequence
|
|
|
1128
|
|
|
1129 my @queries=sort {$hash_query{$b}{"freq"} <=> $hash_query{$a}{"freq"}} keys %hash_query;
|
|
|
1130 my $mature_query=$queries[0];
|
|
|
1131 return $mature_query;
|
|
|
1132 }
|
|
|
1133
|
|
|
1134
|
|
|
1135
|
|
|
1136
|
|
|
1137 sub reset_variables{
|
|
|
1138
|
|
|
1139 #resets the hashes for the next potential precursor
|
|
|
1140
|
|
|
1141 # %hash_query=();
|
|
|
1142 # %hash_comp=();
|
|
|
1143 # %hash_bp=();
|
|
|
1144 foreach my $key (keys %hash_query) {delete($hash_query{$key});}
|
|
|
1145 foreach my $key (keys %hash_comp) {delete($hash_comp{$key});}
|
|
|
1146 foreach my $key (keys %hash_bp) {delete($hash_bp{$key});}
|
|
|
1147
|
|
|
1148 # $message_filter=();
|
|
|
1149 # $message_score=();
|
|
|
1150 # $lines=();
|
|
|
1151 undef($message_filter);
|
|
|
1152 undef($message_score);
|
|
|
1153 undef($lines);
|
|
|
1154 return;
|
|
|
1155 }
|
|
|
1156
|
|
|
1157
|
|
|
1158
|
|
|
1159 sub excise_seq{
|
|
|
1160
|
|
|
1161 #excise sub sequence from the potential precursor
|
|
|
1162
|
|
|
1163 my($seq,$beg,$end,$strand)=@_;
|
|
|
1164
|
|
|
1165 #begin can be equal to end if only one nucleotide is excised
|
|
|
1166 unless($beg<=$end){print STDERR "begin can not be smaller than end for $subject_old\n";exit;}
|
|
|
1167
|
|
|
1168 #rarely, permuted combinations of signature and structure cause out of bound excision errors.
|
|
|
1169 #this happens once appr. every two thousand combinations
|
|
|
1170 unless($beg<=length($seq)){$out_of_bound++;return 0;}
|
|
|
1171
|
|
|
1172 #if on the minus strand, the reverse complement should be excised
|
|
|
1173 if($strand eq "-"){$seq=revcom($seq);}
|
|
|
1174
|
|
|
1175 #the blast parsed format is 1-indexed, substr is 0-indexed
|
|
|
1176 my $sub_seq=substr($seq,$beg-1,$end-$beg+1);
|
|
|
1177
|
|
|
1178 return $sub_seq;
|
|
|
1179
|
|
|
1180 }
|
|
|
1181
|
|
|
1182 sub excise_struct{
|
|
|
1183
|
|
|
1184 #excise sub structure
|
|
|
1185
|
|
|
1186 my($struct,$beg,$end,$strand)=@_;
|
|
|
1187 my $lng=length $struct;
|
|
|
1188
|
|
|
1189 #begin can be equal to end if only one nucleotide is excised
|
|
|
1190 unless($beg<=$end){print STDERR "begin can not be smaller than end for $subject_old\n";exit;}
|
|
|
1191
|
|
|
1192 #rarely, permuted combinations of signature and structure cause out of bound excision errors.
|
|
|
1193 #this happens once appr. every two thousand combinations
|
|
|
1194 unless($beg<=length($struct)){return 0;}
|
|
|
1195
|
|
|
1196 #if excising relative to minus strand, positions are reversed
|
|
|
1197 if($strand eq "-"){($beg,$end)=rev_pos($beg,$end,$lng);}
|
|
|
1198
|
|
|
1199 #the blast parsed format is 1-indexed, substr is 0-indexed
|
|
|
1200 my $sub_struct=substr($struct,$beg-1,$end-$beg+1);
|
|
|
1201
|
|
|
1202 return $sub_struct;
|
|
|
1203 }
|
|
|
1204
|
|
|
1205
|
|
|
1206 sub create_hash_nuclei{
|
|
|
1207 #parses a fasta file with sequences of known miRNAs considered for conservation purposes
|
|
|
1208 #reads the nuclei into a hash
|
|
|
1209
|
|
|
1210 my ($file) = @_;
|
|
|
1211 my ($id, $desc, $sequence, $nucleus) = ();
|
|
|
1212
|
|
|
1213 open (FASTA, "<$file") or die "can not open $file\n";
|
|
|
1214 while (<FASTA>)
|
|
|
1215 {
|
|
|
1216 chomp;
|
|
|
1217 if (/^>(\S+)(.*)/)
|
|
|
1218 {
|
|
|
1219 $id = $1;
|
|
|
1220 $desc = $2;
|
|
|
1221 $sequence = "";
|
|
|
1222 $nucleus = "";
|
|
|
1223 while (<FASTA>){
|
|
|
1224 chomp;
|
|
|
1225 if (/^>(\S+)(.*)/){
|
|
|
1226 $nucleus = substr($sequence,1,$nucleus_lng);
|
|
|
1227 $nucleus =~ tr/[T]/[U]/;
|
|
|
1228 $hash_mirs{$nucleus} .="$id\t";
|
|
|
1229 $hash_nuclei{$nucleus} += 1;
|
|
|
1230
|
|
|
1231 $id = $1;
|
|
|
1232 $desc = $2;
|
|
|
1233 $sequence = "";
|
|
|
1234 $nucleus = "";
|
|
|
1235 next;
|
|
|
1236 }
|
|
|
1237 $sequence .= $_;
|
|
|
1238 }
|
|
|
1239 }
|
|
|
1240 }
|
|
|
1241 $nucleus = substr($sequence,1,$nucleus_lng);
|
|
|
1242 $nucleus =~ tr/[T]/[U]/;
|
|
|
1243 $hash_mirs{$nucleus} .="$id\t";
|
|
|
1244 $hash_nuclei{$nucleus} += 1;
|
|
|
1245 close FASTA;
|
|
|
1246 }
|
|
|
1247
|
|
|
1248
|
|
|
1249 sub parse_file_struct{
|
|
|
1250 #parses the output from RNAfoldand reads it into hashes
|
|
|
1251 my($file) = @_;
|
|
|
1252 my($id,$desc,$seq,$struct,$mfe) = ();
|
|
|
1253 open (FILE_STRUCT, "<$file") or die "can not open $file\n";
|
|
|
1254 while (<FILE_STRUCT>){
|
|
|
1255 chomp;
|
|
|
1256 if (/^>(\S+)\s*(.*)/){
|
|
|
1257 $id= $1;
|
|
|
1258 $desc= $2;
|
|
|
1259 $seq= "";
|
|
|
1260 $struct= "";
|
|
|
1261 $mfe= "";
|
|
|
1262 while (<FILE_STRUCT>){
|
|
|
1263 chomp;
|
|
|
1264 if (/^>(\S+)\s*(.*)/){
|
|
|
1265 $hash_desc{$id} = $desc;
|
|
|
1266 $hash_seq{$id} = $seq;
|
|
|
1267 $hash_struct{$id} = $struct;
|
|
|
1268 $hash_mfe{$id} = $mfe;
|
|
|
1269 $id = $1;
|
|
|
1270 $desc = $2;
|
|
|
1271 $seq = "";
|
|
|
1272 $struct = "";
|
|
|
1273 $mfe = "";
|
|
|
1274 next;
|
|
|
1275 }
|
|
|
1276 if(/^\w/){
|
|
|
1277 tr/uU/tT/;
|
|
|
1278 $seq .= $_;
|
|
|
1279 next;
|
|
|
1280 }
|
|
|
1281 if(/((\.|\(|\))+)/){$struct .=$1;}
|
|
|
1282 if(/\((\s*-\d+\.\d+)\)/){$mfe = $1;}
|
|
|
1283 }
|
|
|
1284 }
|
|
|
1285 }
|
|
|
1286 $hash_desc{$id} = $desc;
|
|
|
1287 $hash_seq{$id} = $seq;
|
|
|
1288 $hash_struct{$id} = $struct;
|
|
|
1289 $hash_mfe{$id} = $mfe;
|
|
|
1290 close FILE_STRUCT;
|
|
|
1291 return;
|
|
|
1292 }
|
|
|
1293
|
|
|
1294
|
|
|
1295 sub score_s{
|
|
|
1296
|
|
|
1297 #this score message is appended to the end of the string of score messages outputted for the potential precursor
|
|
|
1298
|
|
|
1299 my $message=shift;
|
|
|
1300 $message_score.=$message."\n";;
|
|
|
1301 return;
|
|
|
1302 }
|
|
|
1303
|
|
|
1304
|
|
|
1305
|
|
|
1306 sub score_p{
|
|
|
1307
|
|
|
1308 #this score message is appended to the beginning of the string of score messages outputted for the potential precursor
|
|
|
1309
|
|
|
1310 my $message=shift;
|
|
|
1311 $message_score=$message."\n".$message_score;
|
|
|
1312 return;
|
|
|
1313 }
|
|
|
1314
|
|
|
1315
|
|
|
1316
|
|
|
1317 sub filter_s{
|
|
|
1318
|
|
|
1319 #this filtering message is appended to the end of the string of filtering messages outputted for the potential precursor
|
|
|
1320
|
|
|
1321 my $message=shift;
|
|
|
1322 $message_filter.=$message."\n";
|
|
|
1323 return;
|
|
|
1324 }
|
|
|
1325
|
|
|
1326
|
|
|
1327 sub filter_p{
|
|
|
1328
|
|
|
1329 #this filtering message is appended to the beginning of the string of filtering messages outputted for the potential precursor
|
|
|
1330
|
|
|
1331 my $message=shift;
|
|
|
1332 if(defined $message_filter){$message_filter=$message."\n".$message_filter;}
|
|
|
1333 else{$message_filter=$message."\n";}
|
|
|
1334 return;
|
|
|
1335 }
|
|
|
1336
|
|
|
1337
|
|
|
1338 sub find_freq{
|
|
|
1339
|
|
|
1340 #finds the frequency of a given read query from its id.
|
|
|
1341
|
|
|
1342 my($query)=@_;
|
|
|
1343
|
|
|
1344 if($query=~/x(\d+)/i){
|
|
|
1345 my $freq=$1;
|
|
|
1346 return $freq;
|
|
|
1347 }else{
|
|
44
|
1348 #print STDERR "Problem with read format\n";
|
|
12
|
1349 return 0;
|
|
|
1350 }
|
|
|
1351 }
|
|
|
1352
|
|
|
1353
|
|
|
1354 sub print_hash_comp{
|
|
|
1355
|
|
|
1356 #prints the 'comp' hash
|
|
|
1357
|
|
|
1358 my @keys=sort keys %hash_comp;
|
|
|
1359 foreach my $key(@keys){
|
|
|
1360 my $value=$hash_comp{$key};
|
|
|
1361 print "$key \t$value\n";
|
|
|
1362 }
|
|
|
1363 }
|
|
|
1364
|
|
|
1365
|
|
|
1366
|
|
|
1367 sub print_hash_bp{
|
|
|
1368
|
|
|
1369 #prints the 'bp' hash
|
|
|
1370
|
|
|
1371 my @keys=sort {$a<=>$b} keys %hash_bp;
|
|
|
1372 foreach my $key(@keys){
|
|
|
1373 my $value=$hash_bp{$key};
|
|
|
1374 print "$key\t$value\n";
|
|
|
1375 }
|
|
|
1376 print "\n";
|
|
|
1377 }
|
|
|
1378
|
|
|
1379
|
|
|
1380
|
|
|
1381 sub find_strand{
|
|
|
1382
|
|
|
1383 #A subroutine to find the strand, parsing different blast formats
|
|
|
1384
|
|
|
1385 my($other)=@_;
|
|
|
1386
|
|
|
1387 my $strand="+";
|
|
|
1388
|
|
|
1389 if($other=~/-/){
|
|
|
1390 $strand="-";
|
|
|
1391 }
|
|
|
1392
|
|
|
1393 if($other=~/minus/i){
|
|
|
1394 $strand="-";
|
|
|
1395 }
|
|
|
1396 return($strand);
|
|
|
1397 }
|
|
|
1398
|
|
|
1399
|
|
|
1400 sub contained{
|
|
|
1401
|
|
|
1402 #Is the stretch defined by the first positions contained in the stretch defined by the second?
|
|
|
1403
|
|
|
1404 my($beg1,$end1,$beg2,$end2)=@_;
|
|
|
1405
|
|
|
1406 testbeginend($beg1,$end1,$beg2,$end2);
|
|
|
1407
|
|
|
1408 if($beg2<=$beg1 and $end1<=$end2){
|
|
|
1409 return 1;
|
|
|
1410 }else{
|
|
|
1411 return 0;
|
|
|
1412 }
|
|
|
1413 }
|
|
|
1414
|
|
|
1415
|
|
|
1416 sub testbeginend{
|
|
|
1417
|
|
|
1418 #Are the beginposition numerically smaller than the endposition for each pair?
|
|
|
1419
|
|
|
1420 my($begin1,$end1,$begin2,$end2)=@_;
|
|
|
1421
|
|
|
1422 unless($begin1<=$end1 and $begin2<=$end2){
|
|
|
1423 print STDERR "beg can not be larger than end for $subject_old\n";
|
|
|
1424 exit;
|
|
|
1425 }
|
|
|
1426 }
|
|
|
1427
|
|
|
1428
|
|
|
1429 sub rev_pos{
|
|
|
1430
|
|
|
1431 # The blast_parsed format always uses positions that are relative to the 5' of the given strand
|
|
|
1432 # This means that for a sequence of length n, the first nucleotide on the minus strand base pairs with
|
|
|
1433 # the n't nucleotide on the plus strand
|
|
|
1434
|
|
|
1435 # This subroutine reverses the begin and end positions of positions of the minus strand so that they
|
|
|
1436 # are relative to the 5' end of the plus strand
|
|
|
1437
|
|
|
1438 my($beg,$end,$lng)=@_;
|
|
|
1439
|
|
|
1440 my $new_end=$lng-$beg+1;
|
|
|
1441 my $new_beg=$lng-$end+1;
|
|
|
1442
|
|
|
1443 return($new_beg,$new_end);
|
|
|
1444 }
|
|
|
1445
|
|
|
1446 sub round {
|
|
|
1447
|
|
|
1448 #rounds to nearest integer
|
|
|
1449
|
|
|
1450 my($number) = shift;
|
|
|
1451 return int($number + .5);
|
|
|
1452
|
|
|
1453 }
|
|
|
1454
|
|
|
1455
|
|
|
1456 sub rev{
|
|
|
1457
|
|
|
1458 #reverses the order of nucleotides in a sequence
|
|
|
1459
|
|
|
1460 my($sequence)=@_;
|
|
|
1461
|
|
|
1462 my $rev=reverse $sequence;
|
|
|
1463
|
|
|
1464 return $rev;
|
|
|
1465 }
|
|
|
1466
|
|
|
1467 sub com{
|
|
|
1468
|
|
|
1469 #the complementary of a sequence
|
|
|
1470
|
|
|
1471 my($sequence)=@_;
|
|
|
1472
|
|
|
1473 $sequence=~tr/acgtuACGTU/TGCAATGCAA/;
|
|
|
1474
|
|
|
1475 return $sequence;
|
|
|
1476 }
|
|
|
1477
|
|
|
1478 sub revcom{
|
|
|
1479
|
|
|
1480 #reverse complement
|
|
|
1481
|
|
|
1482 my($sequence)=@_;
|
|
|
1483
|
|
|
1484 my $revcom=rev(com($sequence));
|
|
|
1485
|
|
|
1486 return $revcom;
|
|
|
1487 }
|
|
|
1488
|
|
|
1489
|
|
|
1490 sub max2 {
|
|
|
1491
|
|
|
1492 #max of two numbers
|
|
|
1493
|
|
|
1494 my($a, $b) = @_;
|
|
|
1495 return ($a>$b ? $a : $b);
|
|
|
1496 }
|
|
|
1497
|
|
|
1498 sub min2 {
|
|
|
1499
|
|
|
1500 #min of two numbers
|
|
|
1501
|
|
|
1502 my($a, $b) = @_;
|
|
|
1503 return ($a<$b ? $a : $b);
|
|
|
1504 }
|
|
|
1505
|
|
|
1506
|
|
|
1507
|
|
|
1508 sub score_freq{
|
|
|
1509
|
|
|
1510 # scores the count of reads that map to the potential precursor
|
|
|
1511 # Assumes geometric distribution as described in methods section of manuscript
|
|
|
1512
|
|
|
1513 my $freq=shift;
|
|
|
1514
|
|
|
1515 #parameters of known precursors and background hairpins
|
|
|
1516 my $parameter_test=0.999;
|
|
|
1517 my $parameter_control=0.6;
|
|
|
1518
|
|
|
1519 #log_odds calculated directly to avoid underflow
|
|
|
1520 my $intercept=log((1-$parameter_test)/(1-$parameter_control));
|
|
|
1521 my $slope=log($parameter_test/$parameter_control);
|
|
|
1522 my $log_odds=$slope*$freq+$intercept;
|
|
|
1523
|
|
|
1524 #if no strong evidence for 3' overhangs, limit the score contribution to 0
|
|
|
1525 unless($options{x} or $hash_comp{"star_read"}){$log_odds=min2($log_odds,0);}
|
|
|
1526
|
|
|
1527 return $log_odds;
|
|
|
1528 }
|
|
|
1529
|
|
|
1530
|
|
|
1531
|
|
|
1532 ##sub score_mfe{
|
|
|
1533
|
|
|
1534 # scores the minimum free energy in kCal/mol of the potential precursor
|
|
|
1535 # Assumes Gumbel distribution as described in methods section of manuscript
|
|
|
1536
|
|
|
1537 ## my $mfe=shift;
|
|
|
1538
|
|
|
1539 #numerical value, minimum 1
|
|
|
1540 ## my $mfe_adj=max2(1,-$mfe);
|
|
|
1541
|
|
|
1542 #parameters of known precursors and background hairpins, scale and location
|
|
|
1543 ## my $prob_test=prob_gumbel_discretized($mfe_adj,5.5,32);
|
|
|
1544 ## my $prob_background=prob_gumbel_discretized($mfe_adj,4.8,23);
|
|
|
1545
|
|
|
1546 ## my $odds=$prob_test/$prob_background;
|
|
|
1547 ## my $log_odds=log($odds);
|
|
|
1548
|
|
|
1549 ## return $log_odds;
|
|
|
1550 ##}
|
|
|
1551
|
|
|
1552 sub score_mfe{
|
|
|
1553 # use bignum;
|
|
|
1554
|
|
|
1555 # scores the minimum free energy in kCal/mol of the potential precursor
|
|
|
1556 # Assumes Gumbel distribution as described in methods section of manuscript
|
|
|
1557
|
|
|
1558 my ($mfe,$mlng)=@_;
|
|
|
1559
|
|
|
1560 #numerical value, minimum 1
|
|
|
1561 my $mfe_adj=max2(1,-$mfe);
|
|
|
1562 my $mfe_adj1=$mfe/$mlng;
|
|
|
1563 #parameters of known precursors and background hairpins, scale and location
|
|
|
1564 my $a=1.339e-12;my $b=2.778e-13;my $c=45.834;
|
|
|
1565 my $ev=$e**($mfe_adj1*$c);
|
|
44
|
1566 #print STDERR "\n***",$ev,"**\t",$ev+$b,"\t";
|
|
12
|
1567 my $log_odds=($a/($b+$ev));
|
|
|
1568
|
|
|
1569
|
|
|
1570 my $prob_test=prob_gumbel_discretized($mfe_adj,5.5,32);
|
|
|
1571 my $prob_background=prob_gumbel_discretized($mfe_adj,4.8,23);
|
|
|
1572
|
|
|
1573 my $odds=$prob_test/$prob_background;
|
|
|
1574 my $log_odds_2=log($odds);
|
|
44
|
1575 #print STDERR "log_odds :",$log_odds,"\t",$log_odds_2,"\n";
|
|
12
|
1576 return $log_odds;
|
|
|
1577 }
|
|
|
1578
|
|
|
1579
|
|
|
1580
|
|
|
1581 sub prob_gumbel_discretized{
|
|
|
1582
|
|
|
1583 # discretized Gumbel distribution, probabilities within windows of 1 kCal/mol
|
|
|
1584 # uses the subroutine that calculates the cdf to find the probabilities
|
|
|
1585
|
|
|
1586 my ($var,$scale,$location)=@_;
|
|
|
1587
|
|
|
1588 my $bound_lower=$var-0.5;
|
|
|
1589 my $bound_upper=$var+0.5;
|
|
|
1590
|
|
|
1591 my $cdf_lower=cdf_gumbel($bound_lower,$scale,$location);
|
|
|
1592 my $cdf_upper=cdf_gumbel($bound_upper,$scale,$location);
|
|
|
1593
|
|
|
1594 my $prob=$cdf_upper-$cdf_lower;
|
|
|
1595
|
|
|
1596 return $prob;
|
|
|
1597 }
|
|
|
1598
|
|
|
1599
|
|
|
1600 sub cdf_gumbel{
|
|
|
1601
|
|
|
1602 # calculates the cumulative distribution function of the Gumbel distribution
|
|
|
1603
|
|
|
1604 my ($var,$scale,$location)=@_;
|
|
|
1605
|
|
|
1606 my $cdf=$e**(-($e**(-($var-$location)/$scale)));
|
|
|
1607
|
|
|
1608 return $cdf;
|
|
|
1609 }
|
|
|
1610
|
