annotate miRDeep_plant.pl @ 18:a79212816cbc draft

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