|
0
|
1 #!/usr/bin/perl
|
|
|
2 #use feature ':5.10';
|
|
|
3 use File::Basename;
|
|
|
4 use lib dirname($0); # search skript directory for modules
|
|
|
5 use strict 'vars';
|
|
|
6 use warnings;
|
|
|
7 use Getopt::Long;
|
|
|
8 use Pod::Usage;
|
|
|
9 #use List::Util qw/ min max /;
|
|
|
10
|
|
|
11 #use StructureLibrary::RNAtoolsConfig;
|
|
|
12 #use StructureLibrary::Sequence;
|
|
|
13
|
|
|
14 use Cwd qw(getcwd abs_path);
|
|
|
15 use File::Temp qw(tempdir);
|
|
|
16 use File::Copy;
|
|
|
17 use POSIX qw(ceil);
|
|
|
18
|
|
|
19 use FindBin;
|
|
|
20 use lib "$FindBin::Bin";
|
|
|
21 #use GraphClust;
|
|
|
22
|
|
|
23 #use vars qw(%CONFIG);
|
|
|
24 #*CONFIG = \%GraphClust::CONFIG;
|
|
|
25
|
|
|
26 =head1 NAME
|
|
|
27
|
|
|
28
|
|
|
29 fasta2shrep_gspan.pl -fasta mysequences.fasta -wins "50,100,150,200" -shift 5 -M 8
|
|
|
30
|
|
|
31 =head1 SYNOPSIS
|
|
|
32
|
|
|
33 Options:
|
|
|
34
|
|
|
35 HELP
|
|
|
36 -help brief help message
|
|
|
37 -man full documentation
|
|
|
38
|
|
|
39 COMPULSORY
|
|
|
40 -fasta <STRING> e.g. "sequence.fasta"
|
|
|
41 All sequences in fasta format.
|
|
|
42
|
|
|
43 OPTIONS
|
|
|
44 -wins [INTEGER] e.g. "50,100,200"
|
|
|
45 A list of window sizes to use.
|
|
|
46 If none are given (empty string ''), then the entire sequence is
|
|
|
47 taken with no windows. Each window > 1 required!
|
|
|
48 -shift <INTEGER> e.g. 20
|
|
|
49 The shift of the window, relative to the window size given in
|
|
|
50 percent. So you give which percent of the window size shall be
|
|
|
51 used for the shift. Of course the shift is rounded down to the
|
|
|
52 nearest whole number.
|
|
|
53 Example 20 % of a window 150 would result in a step size of 30 nt.
|
|
|
54 It is a relative parameter, as you can give different window sizes.
|
|
|
55 If you do not give this parameter there is a default shift of 1 nt.
|
|
|
56 -cue Crop unpaired ends.
|
|
|
57 If you give this flag, then the unpaired ends of each
|
|
|
58 single structure are ignored. E.g. the structure
|
|
|
59 ...(((...))).. becomes just (((...)))
|
|
|
60 -stack Adds stacking information to graphs. This adds an additional
|
|
|
61 vertex (type P) for each pair of stacked base-pairs and four edges
|
|
|
62 (type p) from each of the involved bases to the new vertex.
|
|
|
63 -e <FLOAT> e.g. 5.0
|
|
|
64 Energy range in kcal/mol (RNAshapes)
|
|
|
65 Use only one of -e and -c!
|
|
|
66 -c <INTEGER> e.g. 10
|
|
|
67 Relative energy range, i.e. percentage (%) of MFE energy (RNAshapes)
|
|
|
68 Use only one of -e and -c!
|
|
|
69 -t <INTEGER> [1-5] e.g. 3 OR "3=0,4=100,5=200"
|
|
|
70 The shape type (RNAshapes). Default is 3.
|
|
|
71 With the list format, the shape level can be changed for different window length
|
|
|
72 "4=100" means that shape level 4 is used from length 100nt (window length)
|
|
|
73 The first given length has to be 0! Not continuous given levels are allowed!
|
|
|
74 -M <INTEGER> e.g. 10
|
|
|
75 Max number of shreps that should be taken per window.
|
|
|
76 -u Ignore unstable structures (RNAshapes).
|
|
|
77 This option filters out closed structures with positive free energy.
|
|
|
78 -r Calculate structure probabilities for shreps (RNAshapes)
|
|
|
79 -i <INT> e.g. 10
|
|
|
80 Turn on structure sampling and gives number of sampling iterations.
|
|
|
81 Default no sampling (i=0)
|
|
|
82 -sample-len <INT> e.g. 100
|
|
|
83 Only in sampling mode: Sampling is only used for seqs/windows >= given length,
|
|
|
84 Default: sample all lengths (0), if -i > 0
|
|
|
85 -q Turn on shape probabilities for RNAshapes, no sampling mode allowed
|
|
|
86 -Tp <FLOAT> e.g 0.001
|
|
|
87 Filter cutoff for shape probabilities, applied before -M filter!
|
|
|
88 -seq-graph-win add for each window a graph which contains no structure
|
|
|
89 -seq-graph-t add for each 't #' a graph which contains no structure
|
|
|
90 -seq-graph-alph change the alphabet of unstructured graphs
|
|
|
91 -annotate <STRING> annotation.tab
|
|
|
92 A file with annotations to be added as abstract graphs
|
|
|
93 on the sequence leven (if given) and on the structure
|
|
|
94 (SHREP) level. The format is has the following TAB-delimited
|
|
|
95 columns: SEQID, START, END, NAMESPACE#LABEL.
|
|
|
96 Labels with the same name-space and SEQID form connected
|
|
|
97 components, which is a sequence of label vertices ordered
|
|
|
98 by the START position in the sequence.
|
|
|
99 -abstr Add abstract structure graphs to the single shrep graph
|
|
|
100 instances.
|
|
|
101 -nostr Calculate no structures, only add sequence information,
|
|
|
102 if this is given, then -seq-graph-win AND/OR -seq-graph-t
|
|
|
103 are required.
|
|
|
104 -match-shape <SHAPE>
|
|
|
105 all seqs/windows will be constraint folded into that shape via
|
|
|
106 RNAshapes (if structure is given in another way this struct will be kept),
|
|
|
107 if this shape is not possible within given energy range, produce a
|
|
|
108 specific t graph with only one vertex 'X'. By this the instance
|
|
|
109 becomes very unsimilar to all other graphs (for knn)
|
|
|
110 -vp enable graph computation with viewpoints:
|
|
|
111 svmsgdnspdk will center on those nucleotides that are given
|
|
|
112 via capital letters and ignore those given as lowercase letters
|
|
|
113 -tmp <STRING> e.g. "/scratch/1/sita/tmp"
|
|
|
114 A directory for writing temporary files
|
|
|
115 -o <STRING> e.g. "ProjectX/MySequences/GSPAN/"
|
|
|
116 Output directory for gspan files containing graphs.
|
|
|
117 -group <INTEGER> e.g. 5
|
|
|
118 Combine/group that number of input seqs into 1 gspan file
|
|
|
119 output name is then '<INT>.group.gspan.bz2'
|
|
|
120 -sge Use SGE cluster for each sequence separately
|
|
|
121 -sge-logdir stdout directory for SGE call
|
|
|
122 -sge-errdir sdterr directory for SGE call
|
|
|
123 -stdout send graphs to stdout instead of writing to files
|
|
|
124 -ignore-header don't write fasta id part after first space to gspan
|
|
|
125 -debug additional debug output
|
|
|
126
|
|
|
127
|
|
|
128 DEFAULT VALUES
|
|
|
129 -wins ""
|
|
|
130 -shift 1 nt
|
|
|
131 -c 10
|
|
|
132 -t 3
|
|
|
133 -M 0 # selects all shreps
|
|
|
134 -tmp "/var/tmp/fasta2shrep"
|
|
|
135 -o "CURRENT_DIR/GSPAN/"
|
|
|
136
|
|
|
137 SGE mode
|
|
|
138 -sge-logdir "CURRENT_DIR/GSPAN/SGE_log"
|
|
|
139 -sge-errdir "CURRENT_DIR/GSPAN/SGE_log"
|
|
|
140 -task-id <NUM>
|
|
|
141
|
|
|
142 =head1 DESCRIPTION
|
|
|
143
|
|
|
144 =cut
|
|
|
145
|
|
|
146 ###############################################################################
|
|
|
147 # end handler for temporary directory
|
|
|
148 # adds an error handler that deletes the directory in case of error
|
|
|
149 # SIGUSR{1/2} are sent by the sge prior to the uncatchable SIGKILL if the
|
|
|
150 # option -notify was set
|
|
|
151 ###############################################################################
|
|
|
152 $SIG{'INT'} = 'end_handler';
|
|
|
153 $SIG{'TERM'} = 'end_handler';
|
|
|
154 $SIG{'ABRT'} = 'end_handler';
|
|
|
155 $SIG{'USR1'} = 'end_handler';
|
|
|
156 $SIG{'USR2'} = 'end_handler';
|
|
|
157
|
|
|
158 sub end_handler {
|
|
|
159 print STDERR "signal ", $_[0], " caught, cleaning up temporary files\n";
|
|
|
160
|
|
|
161 # change into home directory. deletion of the temporary directory will
|
|
|
162 # fail if it is the current working directory
|
|
|
163 chdir();
|
|
|
164 File::Temp::cleanup();
|
|
|
165 die();
|
|
|
166 }
|
|
|
167
|
|
|
168 ###############################################################################
|
|
|
169 # PARSE COMMAND LINE OPTIONS
|
|
|
170 ###############################################################################
|
|
|
171
|
|
|
172 # command line options
|
|
|
173 my ( $i_help, $i_man, $i_debug, $i_fas, $i_wins, $i_shift, $i_crop_unpaired_ends,
|
|
|
174 $i_r, $i_e, $i_c, $i_t, $i_u, $i_M, $i_o, $i_sge, $i_jobid, $i_tmp,
|
|
|
175 $i_ignore_seq_header, $i_stacks, $i_stdout, $i_q, $i_T, $i_i,
|
|
|
176 $i_sample_min_length, $i_sge_logDir, $i_sge_errDir, $i_groupsize, $i_annotate,
|
|
|
177 $i_abstr, $i_no_structure, $i_vp, $i_matchShape, $i_rnashapes_binpath );
|
|
|
178
|
|
|
179 my ( $i_add_seq_graph_win, $i_add_seq_graph_t, $i_change_seq_graph_alph );
|
|
|
180
|
|
|
181 my $options = GetOptions(
|
|
|
182 "help" => \$i_help,
|
|
|
183 "man" => \$i_man,
|
|
|
184 "debug" => \$i_debug,
|
|
|
185 "fasta=s" => \$i_fas,
|
|
|
186 "wins=s" => \$i_wins,
|
|
|
187 "shift=f" => \$i_shift,
|
|
|
188 "cue" => \$i_crop_unpaired_ends,
|
|
|
189 "stack" => \$i_stacks,
|
|
|
190 "r" => \$i_r,
|
|
|
191 "e=f" => \$i_e,
|
|
|
192 "c=i" => \$i_c,
|
|
|
193 "t=s" => \$i_t,
|
|
|
194 "u" => \$i_u,
|
|
|
195 "M=i" => \$i_M,
|
|
|
196 "tmp=s" => \$i_tmp,
|
|
|
197 "o=s" => \$i_o,
|
|
|
198 "i=i" => \$i_i,
|
|
|
199 "sample-len=i" => \$i_sample_min_length,
|
|
|
200 "q" => \$i_q,
|
|
|
201 "Tp=f" => \$i_T,
|
|
|
202 "seq-graph-win" => \$i_add_seq_graph_win,
|
|
|
203 "seq-graph-t" => \$i_add_seq_graph_t,
|
|
|
204 "-seq-graph-alph" => \$i_change_seq_graph_alph,
|
|
|
205 "sge" => \$i_sge,
|
|
|
206 "task-id=i" => \$i_jobid,
|
|
|
207 "ignore-header" => \$i_ignore_seq_header,
|
|
|
208 "stdout" => \$i_stdout,
|
|
|
209 "sge-logdir=s" => \$i_sge_logDir,
|
|
|
210 "sge-errdir=s" => \$i_sge_errDir,
|
|
|
211 "group=i" => \$i_groupsize,
|
|
|
212 "annotate=s" => \$i_annotate,
|
|
|
213 "abstr" => \$i_abstr,
|
|
|
214 "nostr" => \$i_no_structure,
|
|
|
215 "vp" => \$i_vp,
|
|
|
216 "match-shape=s" => \$i_matchShape,
|
|
|
217 "rnashapes-bin=s" => \$i_rnashapes_binpath,
|
|
|
218 );
|
|
|
219
|
|
|
220 pod2usage( -exitstatus => 1, -verbose => 1 ) if $i_help;
|
|
|
221 pod2usage( -exitstatus => 0, -verbose => 2 ) if $i_man;
|
|
|
222 ($options) or pod2usage(2);
|
|
|
223
|
|
|
224 # check compulsory options
|
|
|
225 ($i_fas) or pod2usage("Error: the option -fasta is compulsory!\n");
|
|
|
226 ( -e $i_fas ) or pod2usage("Error: no such file - $i_fas!\n");
|
|
|
227 $i_fas = abs_path($i_fas);
|
|
|
228
|
|
|
229 # check other options and set default values
|
|
|
230 pod2usage("Error: either -e OR -c can be given, but NOT both!\n") if ( $i_e && $i_c );
|
|
|
231 ( $i_e or $i_c ) or $i_c = 10; # set -c 10, if neither -e or -c are given
|
|
|
232 ($i_M) or $i_M = 0; # max number of shreps is 0 (=means take all computed)
|
|
|
233 ($i_i) or $i_i = 0; # default no sampling else sampling iterations
|
|
|
234 ($i_sample_min_length) or $i_sample_min_length = 0;
|
|
|
235 pod2usage("\nError: use --sample-len <INT> only with -i <1..INT> !\n") if ( not $i_i and $i_sample_min_length );
|
|
|
236 ($i_T) or $i_T = 0;
|
|
|
237 ($i_q) or $i_q = 0;
|
|
|
238 pod2usage("\nError: Sampling (-i) not possible with shape probabilities (-q)!\n") if ( $i_i and $i_q );
|
|
|
239 ($i_add_seq_graph_win) or $i_add_seq_graph_win = 0;
|
|
|
240 ($i_add_seq_graph_t) or $i_add_seq_graph_t = 0;
|
|
|
241 ($i_change_seq_graph_alph) or $i_change_seq_graph_alph = 0;
|
|
|
242 ($i_no_structure) or $i_no_structure = 0;
|
|
|
243
|
|
|
244 if ($i_change_seq_graph_alph) {
|
|
|
245 ($i_add_seq_graph_t) or ($i_add_seq_graph_win) or pod2usage( "Error: " .
|
|
|
246 "When giving the parameter -seq-graph-alph, then either -seq-graph-t or -seq-graph-win" .
|
|
|
247 " must also be given!\n" );
|
|
|
248 }
|
|
|
249 ( -e $i_annotate ) or pod2usage("Error: no such file - $i_annotate!\n")
|
|
|
250 if ($i_annotate);
|
|
|
251
|
|
|
252 $i_add_seq_graph_t = 1 if ($i_no_structure);
|
|
|
253
|
|
|
254 ($i_t) or $i_t = 3; # default abstraction type is 3
|
|
|
255 my $change_shape_level = 0;
|
|
|
256 my @level_lens = ( -1, -1, -1, -1, -1 ); ## array_idx-1=shape level, value=start length of this level
|
|
|
257
|
|
|
258 if ( $i_t !~ /^\d+$/ ) {
|
|
|
259
|
|
|
260 my @t_minlens = split( ",", $i_t ); ## -t "3=100,4=200"
|
|
|
261
|
|
|
262 foreach my $idx ( 1 .. @t_minlens ) {
|
|
|
263
|
|
|
264 my $level = $t_minlens[ $idx - 1 ]; ## $level = "3=100"
|
|
|
265
|
|
|
266 die "$level Wrong -t format! Example: -t 3=0,4=100,5=200\n" if ( $level !~ /^\d+\=\d+$/ );
|
|
|
267 my @lev_len = split( "=", $level ); ## $level = "3=100"
|
|
|
268
|
|
|
269 die "Wrong -t format! First level given needs to be length 0! Example: -t 3=0,4=100,5=200\n" if ( $idx == 1 && $lev_len[1] != 0 );
|
|
|
270 die "Wrong -t format! Only level 1-5 allowed! Example: -t 3=0,4=100,5=200\n" if ( $lev_len[0] < 1 or $lev_len[0] > 5 );
|
|
|
271 die "Wrong -t format! Length >= 0 expected! Example: -t 3=0,4=100,5=200\n" if ( $lev_len[1] < 0 );
|
|
|
272 $change_shape_level = 1;
|
|
|
273 $level_lens[ $lev_len[0] - 1 ] = $lev_len[1];
|
|
|
274 }
|
|
|
275
|
|
|
276 ($i_debug) and print STDERR "i_t = $i_t - change $change_shape_level - shape level lengths#" . join( ":", @level_lens ) . "#" . join( ":", @t_minlens ) . "#\n";
|
|
|
277 }
|
|
|
278
|
|
|
279 ## checks for match-shape
|
|
|
280 if ($i_matchShape){
|
|
|
281 die "Please provide correct match shape string like '[]'! Exit...\n\n"
|
|
|
282 if ($i_matchShape !~ /^[\[\]_]+$/);
|
|
|
283 ## RNAshapes prodices anyway only 1 structure, no suboptimal structs in match-shape folding
|
|
|
284 $i_M = 1;
|
|
|
285 }
|
|
|
286
|
|
|
287 my $CURRDIR = getcwd;
|
|
|
288
|
|
|
289 # set up tmp directory
|
|
|
290 # default tmp is /var/tmp, usually not NFS mounted!
|
|
|
291 ( defined $i_tmp ) or $i_tmp = '/var/tmp/';
|
|
|
292 my $tmp_template = 'fasta2shrep-XXXXXX';
|
|
|
293
|
|
|
294 # CLEANUP => 1 : automatically delete at exit
|
|
|
295 $i_tmp = tempdir( $tmp_template, DIR => $i_tmp, CLEANUP => 1 );
|
|
|
296
|
|
|
297 # create GSPAN directory when not printing to stdout
|
|
|
298 if ( not $i_stdout ) {
|
|
|
299 if ($i_o) {
|
|
|
300 ( -e $i_o ) or system("mkdir -p $i_o");
|
|
|
301 } else {
|
|
|
302 system("mkdir -p GSPAN");
|
|
|
303 $i_o = $CURRDIR . "/GSPAN/";
|
|
|
304 }
|
|
|
305 }
|
|
|
306
|
|
|
307 ###############################################################################
|
|
|
308 # GLOBAL VARIABLES
|
|
|
309 ###############################################################################
|
|
|
310 $i_rnashapes_binpath or $i_rnashapes_binpath = "RNAshapes";
|
|
|
311 die "Please provide full RNAshapes path (e.g. /usr/bin/RNAshapes) !\n" if (!-f "$i_rnashapes_binpath");
|
|
|
312 my $rnashapes_loc = $i_rnashapes_binpath;
|
|
|
313 if ( !$rnashapes_loc || !-e $rnashapes_loc ) {
|
|
|
314 my $loc = `which RNAshapes`;
|
|
|
315 chomp($loc);
|
|
|
316 die "\nCannot find RNAshapes binary! Exit...\n\n" if ( !$loc );
|
|
|
317 die "\nCannot find RNAshapes binary! Exit...\n\n" if ( !-e $loc );
|
|
|
318 $rnashapes_loc = $loc;
|
|
|
319 }
|
|
|
320
|
|
|
321 my @WINDOWS = ();
|
|
|
322 @WINDOWS = split( ",", $i_wins ) if ($i_wins);
|
|
|
323 my $globalFolding;
|
|
|
324 $globalFolding = 1 unless @WINDOWS;
|
|
|
325 my $CURRUSER = getlogin;
|
|
|
326 my $SEQNO = 0; # used to id the sequences
|
|
|
327 my $GSPANNO = 0; # used for gspan filenames
|
|
|
328 # minimum length of sequence that will work with RNAshapes
|
|
|
329 # sequences with one or two nucleotides will be restricted to sequence-only graphs
|
|
|
330 # this ensures that no sequences are skipped and external info kept synchronized
|
|
|
331 my $GSPAN_SEQ_MINLEN = 3;
|
|
|
332
|
|
|
333 # name spaces
|
|
|
334 my $ABSTRUCT = "AS";
|
|
|
335 ###############################################################################
|
|
|
336 # EXECUTION CODE
|
|
|
337 ###############################################################################
|
|
|
338
|
|
|
339 # read fasta file into hash
|
|
|
340 my ( $headers_aref, $sequences_aref, $metainfo_aref ) = read_fasta_with_nonunique_headers_meta($i_fas);
|
|
|
341
|
|
|
342 # call script again on the sge cluster in a batch job
|
|
|
343 if ($i_sge) {
|
|
|
344 my $SCRIPTNAME = $0;
|
|
|
345 my $CLUSTERSUBMIT = $FindBin::Bin."/fasta2shrep_gspan.sge";
|
|
|
346
|
|
|
347 die "Cannot find SGE submit script $CLUSTERSUBMIT! Exit...\n\n" if ( !-e $CLUSTERSUBMIT );
|
|
|
348
|
|
|
349 my $SGE_jobs = @{$sequences_aref};
|
|
|
350 $SGE_jobs = ceil( @{$sequences_aref} / $i_groupsize ) if ($i_groupsize);
|
|
|
351
|
|
|
352 my $params = "-fasta $i_fas ";
|
|
|
353 $params .= "-wins $i_wins " if ($i_wins);
|
|
|
354 $params .= "-shift $i_shift " if ($i_shift);
|
|
|
355 $params .= "-e $i_e " if ($i_e);
|
|
|
356 $params .= "-c $i_c " if ($i_c);
|
|
|
357 $params .= "-t $i_t " if ($i_t);
|
|
|
358 $params .= "-u " if ($i_u);
|
|
|
359 $params .= "-r " if ($i_r);
|
|
|
360 $params .= "-M $i_M " if ($i_M);
|
|
|
361 $params .= "-o $i_o " if ($i_o);
|
|
|
362 $params .= "-i $i_i " if ($i_i);
|
|
|
363 $params .= "--sample-len $i_sample_min_length" if ($i_sample_min_length);
|
|
|
364 $params .= "-q " if ($i_q);
|
|
|
365 $params .= "-Tp $i_T " if ($i_T);
|
|
|
366 $params .= "--seq-grap-win " if ($i_add_seq_graph_win);
|
|
|
367 $params .= "--seq-graph-t " if ($i_add_seq_graph_t);
|
|
|
368 $params .= "--seq-graph-alph " if ($i_change_seq_graph_alph);
|
|
|
369 $params .= "-ignore-header " if ($i_ignore_seq_header);
|
|
|
370 $params .= "-cue " if ($i_crop_unpaired_ends);
|
|
|
371 $params .= "-stack " if ($i_stacks);
|
|
|
372 $params .= "--group $i_groupsize " if ($i_groupsize);
|
|
|
373 $params .= "-annotate $i_annotate" if ($i_annotate);
|
|
|
374 $params .= "-abstr" if ($i_abstr);
|
|
|
375 $params .= "-nostr " if ($i_no_structure);
|
|
|
376 $params .= "--debug " if ($i_debug);
|
|
|
377 $params .= "-vp " if ($i_vp);
|
|
|
378 $params .= "--match-shape $i_matchShape " if ($i_matchShape);
|
|
|
379
|
|
|
380 print "used script:" . $SCRIPTNAME . "\n";
|
|
|
381 print "used submit script:" . $CLUSTERSUBMIT . "\n";
|
|
|
382
|
|
|
383 $i_sge_logDir = "$i_o/SGE_log" if ( !$i_sge_logDir );
|
|
|
384 mkdir($i_sge_logDir);
|
|
|
385
|
|
|
386 $i_sge_errDir = $i_sge_logDir if ( !$i_sge_errDir );
|
|
|
387 mkdir($i_sge_errDir);
|
|
|
388
|
|
|
389 my $ssh = 1; ## can be used for debug shell script call
|
|
|
390 if ($ssh) {
|
|
|
391 system( "ssh $CURRUSER\@biui.informatik.uni-freiburg.de "
|
|
|
392 . "'export SGE_ROOT=/opt/sge-6.0/; cd $CURRDIR; "
|
|
|
393 . "/opt/sge-6.0/bin/lx24-amd64/qsub -t 1-$SGE_jobs -o $i_sge_errDir/ -e $i_sge_errDir/ "
|
|
|
394 . "$CLUSTERSUBMIT $CURRDIR $SCRIPTNAME \"$params\" ' " );
|
|
|
395 } else {
|
|
|
396 system("$CLUSTERSUBMIT $CURRDIR $SCRIPTNAME '$params' ");
|
|
|
397 }
|
|
|
398 exit;
|
|
|
399 }
|
|
|
400
|
|
|
401 ## compute shreps and gspan, either local or after SGE submission
|
|
|
402
|
|
|
403 # TODO read and process annotations
|
|
|
404
|
|
|
405 my @used_seq_headers;
|
|
|
406 my @used_seqs;
|
|
|
407 my @used_meta;
|
|
|
408 my $group_idx;
|
|
|
409
|
|
|
410 if ($i_jobid) {
|
|
|
411 ## just process the one sequence as given by the jobid number
|
|
|
412 my $used_grouping = 1; ## if no group is given, make 1 seq per job = group=1
|
|
|
413 $used_grouping = $i_groupsize if ($i_groupsize);
|
|
|
414 my $st = ( $i_jobid - 1 ) * $used_grouping + 1;
|
|
|
415 my $end = $st + $used_grouping - 1;
|
|
|
416 $end = @{$sequences_aref} if ( $end > @{$sequences_aref} );
|
|
|
417
|
|
|
418 foreach my $idx ( $st .. $end ) {
|
|
|
419 push( @used_seq_headers, $headers_aref->[ $idx - 1 ] );
|
|
|
420 push( @used_seqs, $sequences_aref->[ $idx - 1 ] );
|
|
|
421 push( @used_meta, $metainfo_aref->[ $idx - 1 ] );
|
|
|
422 }
|
|
|
423
|
|
|
424 $group_idx = $i_jobid;
|
|
|
425
|
|
|
426 #($i_debug) and print STDERR "st $st end $end gr $group_idx job $i_jobid gs $i_groupsize\n";
|
|
|
427
|
|
|
428 $GSPANNO = $i_jobid-1 if (!$i_groupsize);
|
|
|
429
|
|
|
430 } else {
|
|
|
431 ## process all sequences at once
|
|
|
432 @used_seq_headers = @{$headers_aref};
|
|
|
433 @used_seqs = @{$sequences_aref};
|
|
|
434 @used_meta = @{$metainfo_aref};
|
|
|
435 }
|
|
|
436
|
|
|
437 my $out;
|
|
|
438 my $gspanfile;
|
|
|
439 my $out_no_match_shape;
|
|
|
440
|
|
|
441 if ($i_matchShape && !$i_stdout && !$i_groupsize){
|
|
|
442 open($out_no_match_shape,">$i_o/fasta2shrep.no_match");
|
|
|
443 }
|
|
|
444
|
|
|
445 # for each sequence in the fasta file
|
|
|
446 while ( my $seq = shift @used_seqs ) {
|
|
|
447 my $tmp_header = shift @used_seq_headers;
|
|
|
448 my $tmp_meta = shift @used_meta;
|
|
|
449
|
|
|
450 my ( $seq_id, $seq_header ) = ( $tmp_header =~ /(\S+)\s*([\S*\s*]*)/ );
|
|
|
451 $i_ignore_seq_header and $seq_header = '';
|
|
|
452
|
|
|
453 my $seq_fasta = generate_single_fasta_from_sequence_X( $seq_id, $seq );
|
|
|
454 my $seq_len = length($seq);
|
|
|
455
|
|
|
456 # only print sequence graphs for sequences below this threshold
|
|
|
457 my $no_structure_override = ($seq_len < $GSPAN_SEQ_MINLEN) ? 1 : 0;
|
|
|
458
|
|
|
459 $GSPANNO++;
|
|
|
460
|
|
|
461 # set outstream for gspan output to correct file/STDOUT
|
|
|
462 if ($i_stdout) {
|
|
|
463 $out = \*STDOUT;
|
|
|
464 } elsif ( !$i_groupsize ) {
|
|
|
465 $gspanfile = $i_tmp . '/' . $GSPANNO . '.gspan';
|
|
|
466 open( $out, "| bzip2 -f > $gspanfile.bz2" );
|
|
|
467 } elsif ( ( $GSPANNO - 1 ) % $i_groupsize == 0 ) {
|
|
|
468
|
|
|
469 if ( $GSPANNO > 1 ) {
|
|
|
470 close($out);
|
|
|
471 system("mv $gspanfile.bz2 $i_o/$group_idx.group.gspan.bz2");
|
|
|
472 if ($out_no_match_shape){
|
|
|
473 close($out_no_match_shape);
|
|
|
474 system("mv $gspanfile.no_match $i_o/$group_idx.group.gspan.no_match");
|
|
|
475 }
|
|
|
476 }
|
|
|
477
|
|
|
478 if ( !$i_jobid ) {
|
|
|
479 $group_idx = int( ( $GSPANNO - 1 ) / $i_groupsize ) + 1;
|
|
|
480 }
|
|
|
481
|
|
|
482 $gspanfile = "$i_tmp/$group_idx.group.gspan";
|
|
|
483 open( $out, "| bzip2 -f > $gspanfile.bz2" );
|
|
|
484 open( $out_no_match_shape, ">$gspanfile.no_match" ) if ($out_no_match_shape);
|
|
|
485 }
|
|
|
486
|
|
|
487 ## do not use folding windows in special cases
|
|
|
488 if ( $globalFolding || $i_no_structure || $no_structure_override) {
|
|
|
489 @WINDOWS = ();
|
|
|
490 push( @WINDOWS, $seq_len );
|
|
|
491 }
|
|
|
492
|
|
|
493 ##check win sizes for global folding
|
|
|
494 my @WINDOWS_used = sort { $a <=> $b } @WINDOWS;
|
|
|
495 foreach my $w_idx ( 0 .. ( @WINDOWS_used - 1 ) ) {
|
|
|
496 if ( $WINDOWS_used[$w_idx] >= $seq_len && $WINDOWS_used[$w_idx] > 1 ) {
|
|
|
497 if ( $w_idx < $#WINDOWS_used ) {
|
|
|
498 @WINDOWS_used = @WINDOWS_used[ 0 .. $w_idx ];
|
|
|
499 last;
|
|
|
500 }
|
|
|
501 }
|
|
|
502 }
|
|
|
503
|
|
|
504 ## use seq graph only if no shape folding wanted
|
|
|
505 $i_add_seq_graph_t = 1 if ($i_no_structure || $no_structure_override);
|
|
|
506
|
|
|
507 ## no shape info in graphheader if we have fixed structure (according LocaRNA handling)
|
|
|
508 ## use tags #FS or #S for provided structure
|
|
|
509 my $graph_header;
|
|
|
510 my @struct_meta = grep { $_ =~ /#FS/ || $_ =~ /#S/ } keys %{$tmp_meta};
|
|
|
511
|
|
|
512 if (@struct_meta) {
|
|
|
513 $graph_header = getGraphHeader( $seq_id, $seq_header, \@WINDOWS, $i_shift, $i_e, $i_c, $i_t, $i_u, $i_r, $i_M, $i_crop_unpaired_ends, $i_i, $i_sample_min_length, $i_q, $i_T, $seq_len, 1 );
|
|
|
514 } else {
|
|
|
515 $graph_header = getGraphHeader( $seq_id, $seq_header, \@WINDOWS, $i_shift, $i_e, $i_c, $i_t, $i_u, $i_r, $i_M, $i_crop_unpaired_ends, $i_i, $i_sample_min_length, $i_q, $i_T, $seq_len, $i_no_structure || $no_structure_override );
|
|
|
516 }
|
|
|
517
|
|
|
518 print $out $graph_header;
|
|
|
519
|
|
|
520 my $gi = 0; # current graph index
|
|
|
521
|
|
|
522 ## add graph with no structure at all depending on $i_add_seq_graph_t
|
|
|
523 if ($i_add_seq_graph_t) {
|
|
|
524
|
|
|
525 ($gi) = convertSeqWindow( $seq, $seq_len, 1, $gi, $graph_header, $out, $i_annotate, $seq );
|
|
|
526
|
|
|
527 }
|
|
|
528
|
|
|
529 ## encode fixed structure only if provided and structures wanted in general
|
|
|
530 if ( @struct_meta && !$i_no_structure && !$no_structure_override) {
|
|
|
531
|
|
|
532 my $struct_meta = "#FS";
|
|
|
533 $struct_meta = "#S" if ( !exists $tmp_meta->{$struct_meta} );
|
|
|
534
|
|
|
535 my $seq_shrep = [ $tmp_meta->{$struct_meta}, "ENERGY", "0.00", "SHAPE", $struct_meta ];
|
|
|
536 $gi = convertShapeWindow( [$seq_shrep], $seq, $seq_len, 1, $gi, $out,
|
|
|
537 $graph_header, $i_annotate, $i_abstr, $i_crop_unpaired_ends, $i_stacks, $seq );
|
|
|
538 @WINDOWS_used = (); ## no shape folding if we have a fixed structure
|
|
|
539
|
|
|
540 }
|
|
|
541
|
|
|
542 ## ignore RNAshapes folding if wanted, but do correct file move afterwards
|
|
|
543 ## (just "next" does not work due to output)
|
|
|
544 @WINDOWS_used = () if ($i_no_structure or $no_structure_override);
|
|
|
545
|
|
|
546 #for each window size in list
|
|
|
547 foreach my $win_size (@WINDOWS_used) {
|
|
|
548
|
|
|
549 # calculate shift size from percentage
|
|
|
550 my $curr_shift = 1;
|
|
|
551 if ($i_shift) {
|
|
|
552 $curr_shift = ( $i_shift / 100 ) * $win_size;
|
|
|
553 $curr_shift = int($curr_shift); #round down
|
|
|
554 $curr_shift = 1 unless ($curr_shift); # just in case it is 0
|
|
|
555 }
|
|
|
556 ($i_debug) and print STDERR "winsize: $win_size curr_shift: $curr_shift\n";
|
|
|
557 ($i_debug) and print STDERR "\nNext: $seq_id\t winsize:$win_size \n";
|
|
|
558
|
|
|
559 # choose current shape level, depending on $i_t
|
|
|
560 my $curr_t = 0;
|
|
|
561 if ($change_shape_level) {
|
|
|
562 for ( my $i = 0 ; $i < @level_lens ; $i++ ) {
|
|
|
563 $curr_t = $i + 1 if ( $level_lens[$i] != -1 && ( $level_lens[$i] <= $win_size ) );
|
|
|
564 }
|
|
|
565 ($i_debug) and print STDERR "$win_size curr type $curr_t\n";
|
|
|
566 } else {
|
|
|
567 $curr_t = $i_t;
|
|
|
568 }
|
|
|
569
|
|
|
570 my $rnashapesoutput_fh;
|
|
|
571
|
|
|
572 # call RNAshapes and write to $rnashapesoutput_fh
|
|
|
573 $rnashapesoutput_fh = call_RNAshapes( $seq_fasta, $rnashapes_loc, $win_size,
|
|
|
574 $curr_shift, $i_e, $i_c, $curr_t, $i_u, $i_r, $i_q, $i_T, $i_i, $i_sample_min_length, $seq_len, $i_matchShape );
|
|
|
575
|
|
|
576 # read RNAshapes output from $rnashapesoutput_fh and write subgraph
|
|
|
577 # to gspan file
|
|
|
578 my $gi_old = $gi;
|
|
|
579 $gi = convert_RNAshapes_output( $rnashapesoutput_fh, $gi, $i_M, $out, $graph_header,
|
|
|
580 $win_size, $seq_len, $curr_t, $i_annotate, $i_abstr, $i_crop_unpaired_ends, $i_stacks, $seq );
|
|
|
581
|
|
|
582 ## no (match) shape found at all for this seq
|
|
|
583 if ($gi == $gi_old+1){
|
|
|
584 $gi = convertSeqWindow( "X", 1, 1, $gi, $graph_header, $out, $i_annotate, $seq );
|
|
|
585 print $out_no_match_shape $seq_id."\n" if (!$i_stdout && $out_no_match_shape);
|
|
|
586 }
|
|
|
587 } ## foreach WINDOW_used
|
|
|
588
|
|
|
589 if ( !$i_stdout && !$i_groupsize ) {
|
|
|
590 close($out);
|
|
|
591 move "$gspanfile.bz2", "$i_o/$GSPANNO.gspan.bz2";
|
|
|
592 }
|
|
|
593
|
|
|
594 system("rm $seq_fasta");
|
|
|
595
|
|
|
596 } ## while @used_seqs
|
|
|
597
|
|
|
598 if ($i_groupsize) {
|
|
|
599 close($out);
|
|
|
600 system("mv $gspanfile.bz2 $i_o/$group_idx.group.gspan.bz2");
|
|
|
601 close($out_no_match_shape) if ($out_no_match_shape);
|
|
|
602 move "$gspanfile.no_match", "$i_o/$group_idx.group.gspan.no_match" if ($out_no_match_shape);
|
|
|
603 } elsif ( $out_no_match_shape && !$i_stdout && !$i_groupsize ) {
|
|
|
604 close($out_no_match_shape);
|
|
|
605 }
|
|
|
606
|
|
|
607 ###############################################################################
|
|
|
608 # METHODS
|
|
|
609 ###############################################################################
|
|
|
610
|
|
|
611 ############################################################################
|
|
|
612 # Generates fasta file for a single sequence (one-lined-fasta). This
|
|
|
613 # fasta is stored in the temp directory and should be deleted at the end.
|
|
|
614 # Input:
|
|
|
615 # seq_id : the sequence ID
|
|
|
616 # seq : the sequence
|
|
|
617 #
|
|
|
618 # Output:
|
|
|
619 # The fasta file name
|
|
|
620 ############################################################################
|
|
|
621 sub generate_single_fasta_from_sequence_X {
|
|
|
622 my ( $seq_id, $seq ) = @_;
|
|
|
623
|
|
|
624 $seq = uc($seq);
|
|
|
625 $seq =~ tr/T/U/;
|
|
|
626 $seq =~ s/[^AUCGN]/N/g;
|
|
|
627
|
|
|
628 my $outfas = $i_tmp . "/seq_" . $SEQNO++ . ".fasta";
|
|
|
629 my $host = readpipe("hostname");
|
|
|
630 open( FAS, ">$outfas" ) or die "$host Cannot open file $outfas! Exit...\n\n";
|
|
|
631 print FAS ">$seq_id\n$seq";
|
|
|
632 close(FAS);
|
|
|
633
|
|
|
634 return $outfas;
|
|
|
635 }
|
|
|
636
|
|
|
637 ############################################################################
|
|
|
638 # RNAshapes is called with the given input or default parameters.
|
|
|
639 # Input:
|
|
|
640 # seq_fasta : the sequence fasta file
|
|
|
641 # rnashapes_location : the location of the installation files for RNAshapes
|
|
|
642 # win_size : the current window size
|
|
|
643 # shift : the input parameter -shift
|
|
|
644 # e : the input parameter -e
|
|
|
645 # c : the input parameter -c
|
|
|
646 # t : the input parameter -t
|
|
|
647 # u : the input parameter -u
|
|
|
648 # r : the input parameter -r
|
|
|
649 #
|
|
|
650 # Output: none
|
|
|
651 ############################################################################
|
|
|
652 sub call_RNAshapes {
|
|
|
653 my ( $seq_fasta, $rnashapes_location, $win_size, $shift, $e, $c, $t, $u, $r, $q, $T, $i, $sample_length, $seqLen, $matchShape ) = @_;
|
|
|
654 my $FUNCTION = "call_RNAshapes in fasta2shrep_gspan.pl";
|
|
|
655
|
|
|
656 ($seq_fasta) or die("INPUT ERROR in $FUNCTION: the fasta file is compulsory!\n");
|
|
|
657 ($rnashapes_location) or die( "INPUT ERROR in $FUNCTION: the RNAshapes location" . " is compulsory!\n" );
|
|
|
658 die "$rnashapes_location does not exist! Exit...\n\n" if ( !-e $rnashapes_location );
|
|
|
659 my $call = $rnashapes_location . " -o 1 "; # the output format is of type 1
|
|
|
660 $call .= "-q " if ($q);
|
|
|
661 $call .= "-T $T " if ( $q and $T );
|
|
|
662 $call .= "-w $win_size ";
|
|
|
663 $call .= "-W $shift " if ($i_shift);
|
|
|
664 die("ERROR in $FUNCTION: Give only one of the options -c or -e (RNAshapes)!\n")
|
|
|
665
|
|
|
666 if ( $e && $c );
|
|
|
667 $call .= "-e $e " if ($e);
|
|
|
668 $call .= "-c $c " if ($c);
|
|
|
669 $call .= "-t $t " if ($t);
|
|
|
670 $call .= "-u " if ( $u and not $i ); ## not possible in sampling mode
|
|
|
671 $call .= "-r " if ($r);
|
|
|
672 $call .= "-m $matchShape " if ($matchShape);
|
|
|
673
|
|
|
674 ## check is a bit long but : we want to sample if the window is larger than $sample_length and full seq is larger than window or sample_len
|
|
|
675 ## necessary to do sampling a large window is given, sample_length is shorter than window, but seq is longer than sample_len
|
|
|
676 $call .= "-i $i -A " if ( not $q and $i and ( $win_size >= $sample_length and ( $win_size <= $seqLen or $seqLen >= $sample_length ) ) ); ## -A is to omit samples and print only combined shape probs
|
|
|
677 $call .= " < $seq_fasta";
|
|
|
678
|
|
|
679 ($i_debug) and print STDERR "$seqLen $sample_length $win_size $call\n";
|
|
|
680
|
|
|
681 open my $rnashapesoutput, "$call |" or die( "ERROR in $FUNCTION: The following call " . "could not be carried out!\n$call\n" );
|
|
|
682
|
|
|
683 return $rnashapesoutput;
|
|
|
684 }
|
|
|
685
|
|
|
686 ############################################################################
|
|
|
687 # The output of RNAshapes for one sequence and one window size is read
|
|
|
688 # and converted into graph format.
|
|
|
689 # Input:
|
|
|
690 # rnashapeoutput : filehandle for RNAshapes output in format -o 1
|
|
|
691 # curr_gi : current graph index (for vertices)
|
|
|
692 # maxShreps : max number of shreps to convert to graphs
|
|
|
693 # graph_file_hdl : the output handler for the graph file
|
|
|
694 # graphHead : the header line for the complete graph (for sequence)
|
|
|
695 # winSize : current window size in input for RNAshapes
|
|
|
696 # seqLen : full input seq length
|
|
|
697 # used_t :
|
|
|
698 # annotate:
|
|
|
699 # abstr :
|
|
|
700 # cue :
|
|
|
701 # stacks :
|
|
|
702 # orig_seq : the nucleotide sequence as read from fasta
|
|
|
703 #
|
|
|
704 # Output:
|
|
|
705 # The current graph index
|
|
|
706 ############################################################################
|
|
|
707 sub convert_RNAshapes_output {
|
|
|
708 my ( $rnashapesoutput, $curr_gi, $maxShreps, $graph_file_hdl, $graphHead, $winSize, $seqLen, $used_t, $annotate, $abstr, $cue, $stacks, $orig_seq ) = @_;
|
|
|
709
|
|
|
710 ## omit first line in output, contains fasta header line
|
|
|
711 my $line = <$rnashapesoutput>;
|
|
|
712 my $win_shreps = [];
|
|
|
713 my $win_start;
|
|
|
714 my $win_end;
|
|
|
715 my $win_seq;
|
|
|
716 my $win_shrep_count = 0;
|
|
|
717 my $win_sample = 0;
|
|
|
718 my $winHead;
|
|
|
719 my $win_globalFolding;
|
|
|
720 my $win_size_real;
|
|
|
721
|
|
|
722 # reading RNAshapes output
|
|
|
723 while ( $line = <$rnashapesoutput> ) {
|
|
|
724
|
|
|
725 if ( $line =~ /^(\d+)\s+(\d+)$/ ) {
|
|
|
726 ## line: "<start> <end>"
|
|
|
727
|
|
|
728 if ( @{$win_shreps} > 0 ) {
|
|
|
729
|
|
|
730 print $graph_file_hdl $winHead;
|
|
|
731
|
|
|
732 ## remove SHAPE="_" shrep if it exists in $win_shreps
|
|
|
733 ## do this only when have a sequence graph already
|
|
|
734 if ($i_add_seq_graph_t) {
|
|
|
735 my @new_win_shreps = ();
|
|
|
736 map { push( @new_win_shreps, $_ ) if ( $_->[ $#{$_} ] ne "_" ) } @{$win_shreps};
|
|
|
737 $win_shreps = \@new_win_shreps;
|
|
|
738 }
|
|
|
739
|
|
|
740 ## add graph with no structure depending on $i_add_seq_graph_win to win
|
|
|
741 if ($i_add_seq_graph_win ) {
|
|
|
742 ( $curr_gi ) = convertSeqWindow( $win_seq, $win_size_real,
|
|
|
743 $win_start, $curr_gi, $winHead, $graph_file_hdl, $annotate, $orig_seq );
|
|
|
744 }
|
|
|
745
|
|
|
746 $curr_gi = convertShapeWindow( $win_shreps, $win_seq, $win_size_real,
|
|
|
747 $win_start, $curr_gi, $graph_file_hdl, $winHead, $annotate, $abstr,
|
|
|
748 $cue, $stacks, $orig_seq );
|
|
|
749 }
|
|
|
750
|
|
|
751 ## set new window params
|
|
|
752 $win_shreps = [];
|
|
|
753 $win_start = $1;
|
|
|
754 $win_end = $2;
|
|
|
755 $win_shrep_count = 0;
|
|
|
756 $win_size_real = $win_end - $win_start + 1;
|
|
|
757
|
|
|
758 if ( ($win_size_real) >= $seqLen ) {
|
|
|
759 $win_globalFolding = 1;
|
|
|
760 } else {
|
|
|
761 $win_globalFolding = 0;
|
|
|
762 }
|
|
|
763 my $win_center = $win_start + ( ( $win_size_real + 1 ) / 2 );
|
|
|
764
|
|
|
765 $winHead = getWindowHeader( $graphHead, $winSize, $win_start, $win_end, $win_center, $win_globalFolding, $win_sample, $used_t );
|
|
|
766
|
|
|
767 } elsif ( $line =~ /^(\S+)$/ ) {
|
|
|
768 ## line: "CUUAUGAGUAAGGAAAAUAACGAUUCGGGGUGACGCCCGAAUCCUCACUG"
|
|
|
769 $win_seq = uc($1); ## to be 101% sure that we have by default uppercase chars
|
|
|
770
|
|
|
771 } elsif ( $line =~ /^Results for (\d+) iterations:$/ ) {
|
|
|
772 ## line: "Results for 10 iterations:"
|
|
|
773 $win_sample = $1;
|
|
|
774
|
|
|
775 } elsif ( $line =~ /^Shape\s+\S+\s+not found within energy range.*$/ ) {
|
|
|
776 ## line: "Shape [] not found within energy range (-24.75 to -27.50). Try -c or -e to increase range."
|
|
|
777 $win_shreps = [];
|
|
|
778
|
|
|
779 } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+(\S+)$/ ) {
|
|
|
780 ## line:"...((((..(((....)))))))...........(((((......))))) (-10.10) [[]][]"
|
|
|
781 ## take only $maxShreps shreps per window if set
|
|
|
782 next if ( $maxShreps && $win_shrep_count >= $maxShreps );
|
|
|
783 push( @{$win_shreps}, [ $1, "ENERGY", $2, "SHAPE", $3 ] );
|
|
|
784 $win_shrep_count++;
|
|
|
785
|
|
|
786 } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+\((\S+)\)\s+(\S+)$/ ) {
|
|
|
787 ## line:"((((..((((...((.((.((.....)).)).))...))))..))))... (-10.60) (0.7795360) [[[[[]]]]]"
|
|
|
788 ## take only $maxShreps shreps per window if set
|
|
|
789 next if ( $maxShreps && $win_shrep_count >= $maxShreps );
|
|
|
790 push( @{$win_shreps}, [ $1, "ENERGY", $2, "PROB", $3, "SHAPE", $4 ] );
|
|
|
791 $win_shrep_count++;
|
|
|
792
|
|
|
793 } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+(\S+)\s+(\S+)$/ ) {
|
|
|
794 ## line:"((((..((((...((.((.((.....)).)).))...))))..))))... (-10.60) 0.3000000 [[[[[]]]]]"
|
|
|
795 ## take only $maxShreps shreps per window if set
|
|
|
796 next if ( $maxShreps && $win_shrep_count >= $maxShreps );
|
|
|
797 push( @{$win_shreps}, [ $1, "ENERGY", $2, "SHAPEPROB", $3, "SHAPE", $4 ] );
|
|
|
798 $win_shrep_count++;
|
|
|
799
|
|
|
800 } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+\((\S+)\)\s+(\S+)\s+(\S+)$/ ) {
|
|
|
801 ## line:"((((..((((...((.((.((.....)).)).))...))))..))))... (-10.60) (0.7795360) 0.3000000 [[[[[]]]]]"
|
|
|
802 ## take only $maxShreps shreps per window if set
|
|
|
803 next if ( $maxShreps && $win_shrep_count >= $maxShreps );
|
|
|
804 push( @{$win_shreps}, [ $1, "ENERGY", $2, "PROB", $3, "SHAPEPROB", $4, "SHAPE", $5 ] );
|
|
|
805 $win_shrep_count++;
|
|
|
806
|
|
|
807 } else {
|
|
|
808 next if ( $line =~ /^$/ );
|
|
|
809 die "Unexpected shape output format!\nline=$line\n\nExit...\n\n";
|
|
|
810 }
|
|
|
811 }
|
|
|
812
|
|
|
813 ## convert last windows
|
|
|
814 if ( @{$win_shreps} > 0 ) {
|
|
|
815
|
|
|
816 print $graph_file_hdl $winHead;
|
|
|
817
|
|
|
818 ## remove SHAPE="_" shrep if it exists in $win_shreps
|
|
|
819 ## do this only when have a sequence graph already
|
|
|
820 if ($i_add_seq_graph_t) {
|
|
|
821 my @new_win_shreps = ();
|
|
|
822 map { push( @new_win_shreps, $_ ) if ( $_->[ $#{$_} ] ne "_" ) } @{$win_shreps};
|
|
|
823 $win_shreps = \@new_win_shreps;
|
|
|
824 }
|
|
|
825
|
|
|
826 ## add graph with no structure depending on $i_add_seq_graph_win to win
|
|
|
827 if ($i_add_seq_graph_win) {
|
|
|
828 ($curr_gi) = convertSeqWindow( $win_seq, $win_size_real, $win_start,
|
|
|
829 $curr_gi, $winHead, $graph_file_hdl, $annotate, $orig_seq );
|
|
|
830 }
|
|
|
831
|
|
|
832 $curr_gi = convertShapeWindow( $win_shreps, $win_seq, $win_size_real,
|
|
|
833 $win_start, $curr_gi, $graph_file_hdl, $winHead, $annotate, $abstr, $cue,
|
|
|
834 $stacks, $orig_seq );
|
|
|
835 }
|
|
|
836
|
|
|
837 close($rnashapesoutput);
|
|
|
838
|
|
|
839 return $curr_gi + 1; # return the gi (graph index) for the next subgraph
|
|
|
840 }
|
|
|
841
|
|
|
842 ## Sub to create graph for a complete unstructured sequence
|
|
|
843 # TODO: document function
|
|
|
844 sub convertSeqWindow {
|
|
|
845 my ( $win_seq, $win_size_real, $win_start, $curr_gi, $winHead, $graph_file_hdl, $annotate, $orig_seq ) = @_;
|
|
|
846
|
|
|
847 my $seq_shrep;
|
|
|
848
|
|
|
849 # use a different alphabet (lowercase) if option seq-graph-alph set
|
|
|
850 my $seq_graph_sequence = $i_change_seq_graph_alph ? lc($win_seq) : uc($win_seq);
|
|
|
851 $seq_shrep = [ "." x $win_size_real, "ENERGY", "0.00", "SHAPE", "_", "STRUCT", "." x $win_size_real, "SEQ", $seq_graph_sequence ];
|
|
|
852 my $backboneGraph_ref = getBackboneGraph( $seq_graph_sequence, $curr_gi, $win_start, 0, ( $win_size_real - 1 ), $orig_seq );
|
|
|
853
|
|
|
854 print $graph_file_hdl getSeqHeader( $winHead, $seq_shrep );
|
|
|
855 print $graph_file_hdl join( "\n", @{$backboneGraph_ref} ) . "\n";
|
|
|
856
|
|
|
857 $curr_gi += $win_size_real;
|
|
|
858
|
|
|
859 if ($annotate) {
|
|
|
860 ## TODO add annotation graphs
|
|
|
861 }
|
|
|
862
|
|
|
863 return ($curr_gi);
|
|
|
864 }
|
|
|
865
|
|
|
866 ############################################################################
|
|
|
867 # The results for one window are converted in this method to GSPAN graphs.
|
|
|
868 # Input:
|
|
|
869 # win_shreps_aref : the array ref of shreps (dot-bracket structures) for the
|
|
|
870 # current window
|
|
|
871 # win_seq : the nucleotide sequence for the current window
|
|
|
872 # win_size_real : the current (true) window size
|
|
|
873 # win_start : the starting position of win_seq in the original seq (1-n)
|
|
|
874 # curr_gi : the current graph index
|
|
|
875 # graph_file_hdl : the graph file handler
|
|
|
876 # winHead : the header line for the current sequence window
|
|
|
877 # TODO annotate
|
|
|
878 # abstr : input parameter -abstr
|
|
|
879 # cue : input parameter -cue
|
|
|
880 # stacks : input parameter -stacks
|
|
|
881 # orig_seq : the nucleotide sequence as read from fasta
|
|
|
882 #
|
|
|
883 # Output:
|
|
|
884 # The current graph index
|
|
|
885 ############################################################################
|
|
|
886 sub convertShapeWindow {
|
|
|
887 my ( $win_shreps_aref, $win_seq, $win_size_real, $win_start, $curr_gi,
|
|
|
888 $graph_file_hdl, $winHead, $annotate, $abstr, $cue, $stacks, $orig_seq ) = @_;
|
|
|
889
|
|
|
890 ## generate for each shrep a connected component in gspan file
|
|
|
891 foreach my $shrep ( @{$win_shreps_aref} ) {
|
|
|
892
|
|
|
893 # get the current gi as it was at the beginning of this shrep
|
|
|
894 my $shrep_gi = $curr_gi;
|
|
|
895
|
|
|
896 # cut off unpaired ends, if option is given
|
|
|
897 my $shrep_struct = $shrep->[0];
|
|
|
898 my $crop_index_left = 0;
|
|
|
899 my $crop_index_right = length($shrep_struct) - 1;
|
|
|
900
|
|
|
901 # find croping indices with numbering from 0 to n-1
|
|
|
902 if ($cue) {
|
|
|
903 $crop_index_left = index( $shrep_struct, "(" ); # find 1st occ of "("
|
|
|
904 $crop_index_right = rindex( $shrep_struct, ")" ); # find last occ of ")"
|
|
|
905
|
|
|
906 # if the complete window is unpaired, then don't crop
|
|
|
907 if ( $crop_index_left == -1 ) {
|
|
|
908 $crop_index_left = 0;
|
|
|
909 $crop_index_right = length($shrep_struct) - 1;
|
|
|
910 }
|
|
|
911 }
|
|
|
912
|
|
|
913 # create structure graph
|
|
|
914 my $backboneGraph_ref = getBackboneGraph( $win_seq, $curr_gi, $win_start, $crop_index_left, $crop_index_right, $orig_seq );
|
|
|
915 my $structGraph_ref;
|
|
|
916
|
|
|
917 # add both basic edges, stacks and abstract graph
|
|
|
918 if ($abstr) {
|
|
|
919 ( $structGraph_ref, $curr_gi ) = getStructPlusAbstractStructGraph( $shrep, $curr_gi, $win_size_real, $cue, $stacks );
|
|
|
920
|
|
|
921 # just add basic structure graph edges, and stacks
|
|
|
922 } else {
|
|
|
923 ( $structGraph_ref, $curr_gi ) = getStructGraph( $shrep, $curr_gi, $win_size_real, $stacks );
|
|
|
924 }
|
|
|
925
|
|
|
926 # additional information for shrep header: sequence and dot-bracket
|
|
|
927 my $crop_length = $crop_index_right - $crop_index_left + 1;
|
|
|
928 my $shrepheader_struct = substr( $shrep_struct, $crop_index_left, $crop_length );
|
|
|
929 my $shrepheader_seq = substr( $win_seq, $crop_index_left, $crop_length );
|
|
|
930 push( @{$shrep}, 'STRUCT', $shrepheader_struct, 'SEQ', $shrepheader_seq );
|
|
|
931
|
|
|
932 # print structure graph to file
|
|
|
933 print $graph_file_hdl getShapeHeader( $winHead, $shrep );
|
|
|
934 print $graph_file_hdl join( "\n", @{$backboneGraph_ref} ) . "\n";
|
|
|
935
|
|
|
936 # don't print empty array; sgdnspdk breaks with empty lines
|
|
|
937 if ( @{$structGraph_ref} > 0 ) {
|
|
|
938 print $graph_file_hdl join( "\n", @{$structGraph_ref} ) . "\n";
|
|
|
939
|
|
|
940 if ($annotate) {
|
|
|
941
|
|
|
942 # TODO add annotations to shrep structure
|
|
|
943 # $win_start, $win_end (1-n)
|
|
|
944 }
|
|
|
945
|
|
|
946 }
|
|
|
947 }
|
|
|
948
|
|
|
949 return $curr_gi;
|
|
|
950 }
|
|
|
951
|
|
|
952 ###########################################################################
|
|
|
953 # Here the vertices for the nucleotides are generated and also the backbone
|
|
|
954 # edges, which connect the nucleotides in their correct order (i.e order
|
|
|
955 # of the given sequence)
|
|
|
956 # Input:
|
|
|
957 # win_seq : the sequence for the current window
|
|
|
958 # win_start : the starting position for the current window
|
|
|
959 # curr_crop_i_left : the cropping index for the left end of sequence (-cue)
|
|
|
960 # curr_crop_i_right : the cropping index for the right end of sequence (-cue)
|
|
|
961 # orig_seq : the nucleotide sequence as read from fasta
|
|
|
962 #
|
|
|
963 # Output:
|
|
|
964 # The lines of the graph in an array reference to be printed to a file with
|
|
|
965 # one element per line.
|
|
|
966 ############################################################################
|
|
|
967 sub getBackboneGraph {
|
|
|
968 my ( $win_seq, $curr_gi, $win_start, $curr_crop_i_left, $curr_crop_i_right, $orig_seq ) = @_;
|
|
|
969
|
|
|
970 # RNAshapes substitutes T -> U, sequence only does not
|
|
|
971 # thus, just transform all t/T -> u/U to remain consistent
|
|
|
972 $orig_seq =~ tr /tT/uU/;
|
|
|
973 $win_seq =~ tr /tT/uU/;
|
|
|
974
|
|
|
975 my @seq;
|
|
|
976 @seq = split( "", $win_seq );
|
|
|
977
|
|
|
978 # if the vp option is set we need to obtain the original capitalization from
|
|
|
979 # $orig_seq and extract the windowed sequence manually
|
|
|
980 my @seq_vp;
|
|
|
981 if ( defined $i_vp ) {
|
|
|
982 my $win_len = length($win_seq);
|
|
|
983 my $capitalized_win_seq = substr( $orig_seq, $win_start - 1, $win_len );
|
|
|
984 ( uc($capitalized_win_seq) eq uc($win_seq) ) or
|
|
|
985 die( "error: windowed sequence generated due to vp option not equal " .
|
|
|
986 "to sequence reported by RNASHAPES.\n" .
|
|
|
987 "'${capitalized_win_seq}' != '${win_seq}'" );
|
|
|
988 @seq_vp = split( "", $capitalized_win_seq );
|
|
|
989 }
|
|
|
990
|
|
|
991 my @vert = ();
|
|
|
992 my @edg = ();
|
|
|
993
|
|
|
994 my $curr_abs_pos = $win_start + $curr_crop_i_left;
|
|
|
995 $curr_gi += $curr_crop_i_left;
|
|
|
996
|
|
|
997 # set vertice labeled with 'v' or 'V' according to sequence and vp option
|
|
|
998 # when vp set,
|
|
|
999 # uppercase nucleotides are annotated with 'v'
|
|
|
1000 # lowercase nucleotides are annotated with 'V'
|
|
|
1001 my $vertice_label = ( defined $i_vp and ( $seq_vp[$curr_crop_i_left] =~ /[a-z]/ ) ) ? 'V' : 'v';
|
|
|
1002
|
|
|
1003 # create backbone vertice of first nucleotide
|
|
|
1004 push( @vert, join( ' ', $vertice_label, $curr_gi, $seq[$curr_crop_i_left], $curr_abs_pos ) );
|
|
|
1005
|
|
|
1006 foreach my $idx ( ( $curr_crop_i_left + 1 ) .. $curr_crop_i_right ) {
|
|
|
1007 $curr_abs_pos++;
|
|
|
1008 $curr_gi++;
|
|
|
1009
|
|
|
1010 # set vertice label as described above
|
|
|
1011 $vertice_label = ( defined $i_vp and ( $seq_vp[$idx] =~ /[a-z]/ ) ) ? 'V' : 'v';
|
|
|
1012 push( @vert, join( ' ', $vertice_label, $curr_gi, $seq[$idx], $curr_abs_pos ) );
|
|
|
1013 push( @edg, join( ' ', "e", $curr_gi - 1, $curr_gi, '> 1' ) );
|
|
|
1014 }
|
|
|
1015
|
|
|
1016 my @ret = ( @vert, @edg );
|
|
|
1017
|
|
|
1018 return \@ret;
|
|
|
1019 }
|
|
|
1020
|
|
|
1021 ###########################################################################
|
|
|
1022 # This method does the same as getStructGraph, but is extended to identify
|
|
|
1023 # the abstract parts of the structure and also add this to the graph
|
|
|
1024 # via abstract relations, see information below.
|
|
|
1025 # We already have the backbone graph, now the base-pair edges need to be
|
|
|
1026 # added to this graph and the abstract graph is added after this. Finally
|
|
|
1027 # the abstract relations, pointing from the abstract level to the basic
|
|
|
1028 # structure level is added. All the while, we keep track of the current
|
|
|
1029 # index graph. Furthermore, in the basic structure graph, we have vertices
|
|
|
1030 # stacks that connects the four stacked base-pairs involved. (symbols P
|
|
|
1031 # for vertex and p for edges). The abstract vertices are labelled according
|
|
|
1032 # to the given name-space, then #, then the abstract structure type, e.g.
|
|
|
1033 # HL for hairpin-loop. The edges between the structure types are denoted
|
|
|
1034 # with n, and the relation vertex is ^R, with ^r going from the abstract
|
|
|
1035 # structure to the relation and from the relation to the basic structure
|
|
|
1036 # is @r.
|
|
|
1037 #
|
|
|
1038 # Input:
|
|
|
1039 # curr_shrep : the current shrep as dot-bracket structure
|
|
|
1040 # curr_gi : the current graph index
|
|
|
1041 # win_size : the window size
|
|
|
1042 # cue : the input parameter -cue (whether to crop unpaired ends)
|
|
|
1043 # stacks : the input parameter -stacks (whether to add stack information)
|
|
|
1044 #
|
|
|
1045 # Output:
|
|
|
1046 # The structure graph information line-by-line as an array reference and
|
|
|
1047 # the current graph index.
|
|
|
1048 #
|
|
|
1049 ############################################################################
|
|
|
1050 sub getStructPlusAbstractStructGraph {
|
|
|
1051 my ( $curr_shrep, $curr_gi, $win_size, $cue, $stacks ) = @_;
|
|
|
1052
|
|
|
1053 # $win_size = length($curr_shrep->[0]);
|
|
|
1054
|
|
|
1055 my @struct = split( "", $curr_shrep->[0] );
|
|
|
1056
|
|
|
1057 # print "=================testing: new shrep : gi=$curr_gi ===============================\n" if ($i_debug);
|
|
|
1058
|
|
|
1059 # OBJECTS and VARIABLES
|
|
|
1060 # all indices are saved according to the current graph index, $curr_gi,
|
|
|
1061 # which is at the beginning of the current window, so that the nucleotide
|
|
|
1062 # index can be inferred using $idx + $curr_gi
|
|
|
1063 # opening brackets
|
|
|
1064 my @open_blks = (); # open blocks of consecutive opening brackets (array of arrays)
|
|
|
1065 my @p_open = (); # currently open block of the last identified complete block
|
|
|
1066 my @c_open = (); # current "incomplete" block of consecutive opening brackets
|
|
|
1067 # base-pair objects
|
|
|
1068 my %bps = (); # hash of all base-pairs with stem for a value, if it is at
|
|
|
1069 # the BP is opening or closing a stem object
|
|
|
1070 my @c_bp = (); # current BP, just being closed at current index pos (i,j)
|
|
|
1071 my @p_bp = (); # previously closed BP at position before current index
|
|
|
1072 # stem objects
|
|
|
1073 my %stems = (); # all stems, key="i:j,k:l", value=stem array, (i,j) outer BP
|
|
|
1074 my @c_stem = (); # current incomplete stem object
|
|
|
1075 my @stmbrks = (); # all stem-breaks in the form of (i,k,l), where i is the remaining
|
|
|
1076 # open base and (k,l) is the subsequent closing base-pair
|
|
|
1077
|
|
|
1078 # my @p_stmbrk = (); # previous stem-break in the form of (i,k,l), where i is the remaining
|
|
|
1079 # # open base and (k,l) is the subsequent closing base-pair
|
|
|
1080 # loop objects
|
|
|
1081 my @up = (); # all unpaired regions (i,j) that are of unkown type
|
|
|
1082 my @c_up = (); # current incomplete unpaired object
|
|
|
1083 my @hls = (); # all hairpin loop objects
|
|
|
1084 my @bls = (); # all bulge-loop objects
|
|
|
1085 my @ils = (); # all internal loops
|
|
|
1086 my @mls = (); # all multi-loops
|
|
|
1087 my @c_mls = (); # current incomplete multi-loops! (there can be more than one)
|
|
|
1088 my @els = (); # all external loops
|
|
|
1089 my @c_el = (); # current incomplete external loop
|
|
|
1090
|
|
|
1091 # iterate through shrep structure an identify abstract parts
|
|
|
1092 foreach my $idx ( 0 .. @struct - 1 ) {
|
|
|
1093
|
|
|
1094 #===================================================
|
|
|
1095 # current char is the opening bracket of a base-pair
|
|
|
1096 if ( $struct[$idx] eq "(" ) {
|
|
|
1097
|
|
|
1098 #===================================================
|
|
|
1099
|
|
|
1100 # update current window index to current graph index
|
|
|
1101 $idx += $curr_gi;
|
|
|
1102
|
|
|
1103 # case: '((' currently there is an open block, extend it
|
|
|
1104 if (@c_open) {
|
|
|
1105 push( @c_open, $idx );
|
|
|
1106
|
|
|
1107 # case: '.(' or ')(' or BOS'(' there is no open block, create a new one ".(" or ")("
|
|
|
1108 } else {
|
|
|
1109
|
|
|
1110 # begin a new current open block
|
|
|
1111 push( @c_open, $idx );
|
|
|
1112
|
|
|
1113 # case: ')('
|
|
|
1114 if (@c_stem) {
|
|
|
1115
|
|
|
1116 # CLOSE STEM
|
|
|
1117 print STDERR "TEST - closed stem:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_stem ) ) if ($i_debug);
|
|
|
1118 close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps );
|
|
|
1119
|
|
|
1120 # case (x(x)(
|
|
|
1121 if (@open_blks) {
|
|
|
1122
|
|
|
1123 # EXTEND_ML, if one exists
|
|
|
1124 if (@c_mls) {
|
|
|
1125 extend_or_open_nested_ML( $curr_shrep, \@c_mls, \@p_bp, \@up, \@stmbrks, \@open_blks, $idx, $curr_gi );
|
|
|
1126
|
|
|
1127 # OPEN_ML
|
|
|
1128 } else {
|
|
|
1129
|
|
|
1130 # case 'x((x)(' the opening of the ML was a stem-break
|
|
|
1131 if ( @stmbrks
|
|
|
1132 && $stmbrks[-1]->[1] == $p_bp[0]
|
|
|
1133 && $stmbrks[-1]->[2] == $p_bp[1] ) {
|
|
|
1134 my @c_stmbrk = @{ pop(@stmbrks) };
|
|
|
1135 my @a = ( $c_stmbrk[0], -1 );
|
|
|
1136 my @ml = ();
|
|
|
1137 push( @ml, \@a ); # 1st base-pair is in-waiting
|
|
|
1138 my @a2 = ( $c_stmbrk[1], $c_stmbrk[2] );
|
|
|
1139 push( @ml, \@a2 ); #2nd BP
|
|
|
1140 my @a3 = ($idx);
|
|
|
1141 push( @ml, \@a3 ); #first of current BP
|
|
|
1142 push( @c_mls, \@ml );
|
|
|
1143 print STDERR "TEST - opened ML with SB:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
|
|
|
1144
|
|
|
1145 # case 'x(.(x)('
|
|
|
1146 } elsif (@up) {
|
|
|
1147 my @tmp_up = @{ pop(@up) };
|
|
|
1148
|
|
|
1149 # compare adjacent BP to UP to see if they fit
|
|
|
1150 if ( $tmp_up[1] + 1 == $p_bp[0]
|
|
|
1151 && $open_blks[-1]->[-1] == $tmp_up[0] - 1 ) {
|
|
|
1152 my @a = ( $open_blks[-1]->[-1], -1 );
|
|
|
1153 my @ml = ();
|
|
|
1154 push( @ml, \@a ); # 1st BP awaits closing
|
|
|
1155 my @newbp = @p_bp;
|
|
|
1156 push( @ml, \@newbp ); # 2nd BP
|
|
|
1157 my @a2 = ($idx);
|
|
|
1158 push( @ml, \@a2 ); # current open BP, awaits closing
|
|
|
1159 push( @c_mls, \@ml );
|
|
|
1160 print STDERR "TEST - opened ML with prev UP1:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) )
|
|
|
1161
|
|
|
1162 if ($i_debug);
|
|
|
1163 } else {
|
|
|
1164 die "ERROR: the base-pairs to not match the " . "previous unpaired region?!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1165 }
|
|
|
1166 } else {
|
|
|
1167 die "ERROR: in case ML, but there is no initial" . "part of the ML to be identified\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1168 }
|
|
|
1169 }
|
|
|
1170
|
|
|
1171 # case: 'x(x)('
|
|
|
1172 } else {
|
|
|
1173
|
|
|
1174 # case: '.(x)(' => CLOSE_EL
|
|
|
1175 if (@c_el) {
|
|
|
1176 if ( $c_el[-1] == $p_bp[0] ) {
|
|
|
1177 push( @c_el, $p_bp[1] );
|
|
|
1178 my @newel = @c_el;
|
|
|
1179 push( @els, \@newel );
|
|
|
1180 @c_el = ();
|
|
|
1181 print STDERR "TEST - closed EL1:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newel ) ) if ($i_debug);
|
|
|
1182 } else {
|
|
|
1183 die "ERROR: the previous BP must match the current EL\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1184 }
|
|
|
1185 }
|
|
|
1186
|
|
|
1187 # in both cases OPEN_EL
|
|
|
1188 push( @c_el, @p_bp );
|
|
|
1189 push( @c_el, $idx );
|
|
|
1190 print STDERR "TEST - opened EL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_el ) ) if ($i_debug);
|
|
|
1191 }
|
|
|
1192
|
|
|
1193 # case '.('
|
|
|
1194 } elsif (@c_up) {
|
|
|
1195
|
|
|
1196 # case (x(x).( OR 'x(.('
|
|
|
1197 if (@open_blks) {
|
|
|
1198
|
|
|
1199 # case 'x(x).(' ML
|
|
|
1200 if ( @p_bp && ( $p_bp[1] == $c_up[0] - 1 ) ) {
|
|
|
1201
|
|
|
1202 # case '((x)x(x).(' EXTEND_ML, if one exists
|
|
|
1203 if (@c_mls) {
|
|
|
1204
|
|
|
1205 extend_or_open_nested_ML( $curr_shrep, \@c_mls, \@p_bp, \@up, \@stmbrks, \@open_blks, $idx, $curr_gi );
|
|
|
1206
|
|
|
1207 # case 'x((x).(' or 'x(.(x).(' OPEN_ML
|
|
|
1208 } else {
|
|
|
1209
|
|
|
1210 # case 'x((x).(' the opening of the ML was a stem-break
|
|
|
1211 if ( @stmbrks
|
|
|
1212 && $stmbrks[-1]->[1] == $p_bp[0]
|
|
|
1213 && $stmbrks[-1]->[2] == $p_bp[1] ) {
|
|
|
1214 my @c_stmbrk = @{ pop(@stmbrks) };
|
|
|
1215 my @a = ( $c_stmbrk[0], -1 );
|
|
|
1216 my @ml = ();
|
|
|
1217 push( @ml, \@a ); # 1st base-pair is in-waiting
|
|
|
1218 my @a2 = ( $c_stmbrk[1], $c_stmbrk[2] );
|
|
|
1219 push( @ml, \@a2 ); #2nd BP
|
|
|
1220 my @a3 = ($idx);
|
|
|
1221 push( @ml, \@a3 ); #first of current BP
|
|
|
1222 push( @c_mls, \@ml );
|
|
|
1223 print STDERR "TEST - open ML with SB:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
|
|
|
1224
|
|
|
1225 # case 'x(.(x).(' opening the ML with initial UP region
|
|
|
1226 } elsif (@up) {
|
|
|
1227
|
|
|
1228 # get previous unpaired region
|
|
|
1229 my $tmp = pop(@up);
|
|
|
1230 my @p_up = @{$tmp};
|
|
|
1231
|
|
|
1232 # case: '(.(x).(' compare adjacent BP to UP to see if they create a ML
|
|
|
1233 if ( $p_up[1] + 1 == $p_bp[0]
|
|
|
1234 && $open_blks[-1]->[-1] == $p_up[0] - 1 ) {
|
|
|
1235 my @a = ( $open_blks[-1]->[-1], -1 );
|
|
|
1236 my @ml = ();
|
|
|
1237 push( @ml, \@a ); # 1st BP awaits closing
|
|
|
1238 my @newbp = @p_bp;
|
|
|
1239 push( @ml, \@newbp ); # 2nd BP
|
|
|
1240 my @a2 = ($idx);
|
|
|
1241 push( @ml, \@a2 ); # current open BP, awaits closing
|
|
|
1242 push( @c_mls, \@ml );
|
|
|
1243 print STDERR "TEST - opened ML with prev UP2:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) )
|
|
|
1244
|
|
|
1245 if ($i_debug);
|
|
|
1246 }
|
|
|
1247
|
|
|
1248 # has to be an ML, but can't find first part
|
|
|
1249 } else {
|
|
|
1250 die "ERROR: in case ML, but there is no initial " . "part of the ML to be identified\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1251 }
|
|
|
1252 }
|
|
|
1253
|
|
|
1254 # case '(.(
|
|
|
1255 } else {
|
|
|
1256 my @newup = @c_up;
|
|
|
1257 push( @up, \@newup );
|
|
|
1258 print STDERR "TEST - closed UP:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newup ) ) if ($i_debug);
|
|
|
1259 }
|
|
|
1260
|
|
|
1261 # case: 'x(x).(' OR 'EOS.('
|
|
|
1262 } else {
|
|
|
1263
|
|
|
1264 # case: '.(x)(' => CLOSE_EL
|
|
|
1265 if (@c_el) {
|
|
|
1266 if ( $c_el[-1] == $c_bp[0] ) {
|
|
|
1267 push( @c_el, $p_bp[1] );
|
|
|
1268 my @newel = @c_el;
|
|
|
1269 push( @els, \@newel );
|
|
|
1270 @c_el = ();
|
|
|
1271 print STDERR "TEST - closed EL2:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newel ) ) if ($i_debug);
|
|
|
1272 } else {
|
|
|
1273 die "ERROR: the previous BP must match the current EL\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1274 }
|
|
|
1275 }
|
|
|
1276
|
|
|
1277 # in both cases OPEN_EL
|
|
|
1278 if (@p_bp) {
|
|
|
1279 push( @c_el, @p_bp );
|
|
|
1280 push( @c_el, $idx );
|
|
|
1281 } else {
|
|
|
1282
|
|
|
1283 # ignore this if -cue is given
|
|
|
1284 # TODO check the gi index
|
|
|
1285 if ($cue) {
|
|
|
1286 @c_el = ();
|
|
|
1287 } else {
|
|
|
1288 push( @c_el, ( $curr_gi, $curr_gi ) );
|
|
|
1289 push( @c_el, $idx );
|
|
|
1290 print STDERR "TEST - opened EL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_el ) ) if ($i_debug);
|
|
|
1291 }
|
|
|
1292 }
|
|
|
1293 }
|
|
|
1294 @c_up = ();
|
|
|
1295
|
|
|
1296 } else {
|
|
|
1297
|
|
|
1298 # no more cases except opening bracket at beginning of sequence
|
|
|
1299 # in this case do nothing, as the index has already been added
|
|
|
1300 die "ERROR: '((', ')(', '.(' have all been covered\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) )
|
|
|
1301 unless ( $idx == $curr_gi );
|
|
|
1302
|
|
|
1303 }
|
|
|
1304 }
|
|
|
1305
|
|
|
1306 #===================================================
|
|
|
1307 # current char is the unpaired base
|
|
|
1308 } elsif ( $struct[$idx] eq "." ) {
|
|
|
1309
|
|
|
1310 #===================================================
|
|
|
1311
|
|
|
1312 # update current window index to current graph index
|
|
|
1313 $idx += $curr_gi;
|
|
|
1314
|
|
|
1315 # case ').' => CLOSE_STEM
|
|
|
1316 if (@c_stem) {
|
|
|
1317 print STDERR "TEST - closed stem:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_stem ) ) if ($i_debug);
|
|
|
1318 close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps );
|
|
|
1319 }
|
|
|
1320
|
|
|
1321 # case '..' extend UP, @c_up
|
|
|
1322 if (@c_up) {
|
|
|
1323 $c_up[1] = $idx;
|
|
|
1324
|
|
|
1325 # case 'x.' open new UP
|
|
|
1326 } else {
|
|
|
1327
|
|
|
1328 # OPEN_UP, @c_up
|
|
|
1329 @c_up = ( $idx, $idx );
|
|
|
1330
|
|
|
1331 # case '(.' just come to the end of an open block, push it onto stack
|
|
|
1332 if (@c_open) {
|
|
|
1333 my @newopen = @c_open;
|
|
|
1334 push( @open_blks, \@newopen );
|
|
|
1335 @c_open = ();
|
|
|
1336
|
|
|
1337 # case ').' or '.'
|
|
|
1338 # (have already closed stem and created stem-break and emptied p_open)
|
|
|
1339 } else {
|
|
|
1340
|
|
|
1341 # case 'x(x).' or '.' => CLOSE_EL
|
|
|
1342 if ( !@open_blks ) {
|
|
|
1343
|
|
|
1344 # case 'x().().'
|
|
|
1345 if (@c_el) {
|
|
|
1346
|
|
|
1347 # double check
|
|
|
1348 if ( $c_el[-1] == $p_bp[0] ) {
|
|
|
1349 push( @c_el, $p_bp[1] );
|
|
|
1350 my @newel = @c_el;
|
|
|
1351 push( @els, \@newel );
|
|
|
1352 @c_el = ();
|
|
|
1353 print STDERR "TEST - closed EL3:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newel ) ) if ($i_debug);
|
|
|
1354 } else {
|
|
|
1355 die "This case should not occur!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1356 }
|
|
|
1357 }
|
|
|
1358 }
|
|
|
1359 }
|
|
|
1360 }
|
|
|
1361
|
|
|
1362 #===================================================
|
|
|
1363 # current char is the closing bracket of a base-pair
|
|
|
1364 } elsif ( $struct[$idx] eq ")" ) {
|
|
|
1365
|
|
|
1366 #===================================================
|
|
|
1367
|
|
|
1368 # update current window index to current graph index
|
|
|
1369 $idx += $curr_gi;
|
|
|
1370
|
|
|
1371 # save previous base-pair
|
|
|
1372 if (@c_bp) {
|
|
|
1373 @p_bp = @c_bp;
|
|
|
1374 }
|
|
|
1375
|
|
|
1376 # case: '((x))', extend stem
|
|
|
1377 if (@p_open) {
|
|
|
1378
|
|
|
1379 # get current base-pair
|
|
|
1380 @c_bp = ( pop(@p_open), $idx );
|
|
|
1381
|
|
|
1382 # add it to the base-pair hash
|
|
|
1383 $bps{"$c_bp[0]:$c_bp[1]"} = "";
|
|
|
1384
|
|
|
1385 # EXTEND_STEM
|
|
|
1386 if (@c_stem) {
|
|
|
1387 $c_stem[0] = $c_bp[0];
|
|
|
1388 $c_stem[1] = $c_bp[1];
|
|
|
1389 } else {
|
|
|
1390 die "ERROR: in case '))' there has to be an open stem object!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1391 }
|
|
|
1392
|
|
|
1393 # case: '(x(x))' or '(x.)', '()' not allowed
|
|
|
1394 } else {
|
|
|
1395
|
|
|
1396 # case: '(x(x))' CLOSE_STEM close previous stem
|
|
|
1397 if (@c_stem) {
|
|
|
1398 print STDERR "TEST - closed stem:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_stem ) ) if ($i_debug);
|
|
|
1399 close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps );
|
|
|
1400
|
|
|
1401 # case '(x(x))' OPEN_STEM
|
|
|
1402 if (@open_blks) {
|
|
|
1403 my $tmp = pop(@open_blks);
|
|
|
1404 @p_open = @{$tmp};
|
|
|
1405
|
|
|
1406 # get current base-pair
|
|
|
1407 @c_bp = ( pop(@p_open), $idx );
|
|
|
1408
|
|
|
1409 # add it to the base-pair hash
|
|
|
1410 # as the stem is not finished yet, cannot add stem information
|
|
|
1411 $bps{"$c_bp[0]:$c_bp[1]"} = "";
|
|
|
1412
|
|
|
1413 # open stem
|
|
|
1414 $c_stem[0] = $c_bp[0];
|
|
|
1415 $c_stem[1] = $c_bp[1];
|
|
|
1416 $c_stem[2] = $c_bp[0];
|
|
|
1417 $c_stem[3] = $c_bp[1];
|
|
|
1418 } else {
|
|
|
1419 die "ERROR: there have to be open blocks to match " . "current closing bracket!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1420 }
|
|
|
1421
|
|
|
1422 # case: '(.(x))' => CLOSE_BL
|
|
|
1423 if ( @up && $up[-1]->[0] - 1 == $c_bp[0] && $up[-1]->[1] + 1 == $p_bp[0] ) {
|
|
|
1424 my @newBL = ( @c_bp, @p_bp );
|
|
|
1425 push( @bls, \@newBL );
|
|
|
1426 pop(@up);
|
|
|
1427 print STDERR "TEST - closed BL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newBL ) ) if ($i_debug);
|
|
|
1428
|
|
|
1429 # case: '(x(x).(x))' => CLOSE_ML
|
|
|
1430 } elsif ( @c_mls && $c_mls[-1]->[-1]->[-1] == $p_bp[0] ) {
|
|
|
1431 my @newml = @{ pop(@c_mls) };
|
|
|
1432 push( @{ $newml[-1] }, $p_bp[1] ); # close last BP
|
|
|
1433 $newml[0]->[1] = $c_bp[1]; # close 1st BP
|
|
|
1434 push( @mls, \@newml );
|
|
|
1435 print STDERR "TEST - closed ML1:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newml, 1 ) ) if ($i_debug);
|
|
|
1436 } else {
|
|
|
1437 die "ERROR: what is this case?\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1438 }
|
|
|
1439
|
|
|
1440 # case: (x.)
|
|
|
1441 } else {
|
|
|
1442
|
|
|
1443 # OPEN_STEM open new stem
|
|
|
1444 if (@open_blks) {
|
|
|
1445 @p_open = @{ pop(@open_blks) };
|
|
|
1446
|
|
|
1447 # get current base-pair
|
|
|
1448 @c_bp = ( pop(@p_open), $idx );
|
|
|
1449
|
|
|
1450 # add it to the base-pair hash
|
|
|
1451 # as the stem is not finished yet, cannot add stem information
|
|
|
1452 $bps{"$c_bp[0]:$c_bp[1]"} = "";
|
|
|
1453
|
|
|
1454 # open stem
|
|
|
1455 $c_stem[0] = $c_bp[0];
|
|
|
1456 $c_stem[1] = $c_bp[1];
|
|
|
1457 $c_stem[2] = $c_bp[0];
|
|
|
1458 $c_stem[3] = $c_bp[1];
|
|
|
1459 } else {
|
|
|
1460 die "ERROR: there have to be open blocks to match " . "current closing bracket!$curr_shrep->[0]\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1461 }
|
|
|
1462
|
|
|
1463 # CLOSE_C_UP
|
|
|
1464 if (@c_up) {
|
|
|
1465
|
|
|
1466 # case: (.) => close_HL
|
|
|
1467 if ( $c_up[0] - 1 == $c_bp[0] ) {
|
|
|
1468 my @newhl = @c_bp;
|
|
|
1469 push( @hls, \@newhl );
|
|
|
1470 @c_up = ();
|
|
|
1471 print STDERR "TEST - closed HL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newhl ) ) if ($i_debug);
|
|
|
1472
|
|
|
1473 # case: (x(x).)
|
|
|
1474 } elsif ( @p_bp && $c_up[0] - 1 == $p_bp[1] ) {
|
|
|
1475
|
|
|
1476 # case: ((x).) => CLOSE_BL
|
|
|
1477 if ( @stmbrks
|
|
|
1478 && $stmbrks[-1]->[0] == $c_bp[0]
|
|
|
1479 && $stmbrks[-1]->[1] == $p_bp[0]
|
|
|
1480 && $stmbrks[-1]->[2] == $p_bp[1] ) {
|
|
|
1481 my @newbl = ( @c_bp, @p_bp );
|
|
|
1482 push( @bls, \@newbl );
|
|
|
1483 pop(@stmbrks);
|
|
|
1484 @c_up = ();
|
|
|
1485 print STDERR "TEST - closed BL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newbl ) ) if ($i_debug);
|
|
|
1486
|
|
|
1487 # case: (x(x).), either CLOSE_ML or CLOSE_IL
|
|
|
1488 } else {
|
|
|
1489
|
|
|
1490 # case '(.(x).)' => CLOSE_IL
|
|
|
1491 if ( @up
|
|
|
1492 && ( $up[-1]->[0] - 1 == $c_bp[0] )
|
|
|
1493 && ( $up[-1]->[1] + 1 == $p_bp[0] ) ) {
|
|
|
1494 my @p_up = @{ pop(@up) };
|
|
|
1495 my @newil = ( @c_bp, @p_bp );
|
|
|
1496 push( @ils, \@newil );
|
|
|
1497 @c_up = ();
|
|
|
1498 print STDERR "TEST - closed IL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newil ) ) if ($i_debug);
|
|
|
1499
|
|
|
1500 # case (x(x)x(x).) => CLOSE_ML remember, array of tuples
|
|
|
1501 } elsif ( @c_mls && $c_mls[-1]->[-1]->[0] == $p_bp[0] ) {
|
|
|
1502 my @newml = @{ pop(@c_mls) };
|
|
|
1503 push( @{ $newml[-1] }, $p_bp[1] ); # close final BP
|
|
|
1504 $newml[0]->[1] = $c_bp[1]; # close 1st BP
|
|
|
1505 push( @mls, \@newml );
|
|
|
1506 @c_up = ();
|
|
|
1507 print STDERR "TEST - closed ML2:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newml, 1 ) ) if ($i_debug);
|
|
|
1508 } else {
|
|
|
1509 die "There should be no such case!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1510 }
|
|
|
1511 }
|
|
|
1512
|
|
|
1513 # there should be no more unknown unpaired regions left!
|
|
|
1514 } else {
|
|
|
1515 die "There shouldn't be any unknown unpaired " . "regions left\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1516 }
|
|
|
1517
|
|
|
1518 } else {
|
|
|
1519 die "ERROR: '()' is not allowed in $curr_shrep->[0]\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
|
|
|
1520 }
|
|
|
1521 }
|
|
|
1522 }
|
|
|
1523 }
|
|
|
1524 }
|
|
|
1525
|
|
|
1526 # EOS
|
|
|
1527
|
|
|
1528 # CLOSE_STEM
|
|
|
1529 print STDERR "TEST - closed stem EOS:\n" . ( test( $curr_shrep->[0], ( $win_size + $curr_gi ), $curr_gi, \@c_stem ) )
|
|
|
1530 if ( @c_stem && $i_debug );
|
|
|
1531 close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps ) if (@c_stem);
|
|
|
1532 die "ERROR: there should not be an open ML!\n" . $curr_shrep . "\n"
|
|
|
1533 if (@c_mls);
|
|
|
1534
|
|
|
1535 # case ().() CLOSE_EL
|
|
|
1536 if (@c_el) {
|
|
|
1537 if ( $c_el[-1] == $p_bp[0] ) {
|
|
|
1538 push( @c_el, $p_bp[1] );
|
|
|
1539 my @newel = @c_el;
|
|
|
1540 push( @els, \@newel );
|
|
|
1541 @c_el = ();
|
|
|
1542
|
|
|
1543 print STDERR "TEST - closed EL4 EOS:\n" . ( test( $curr_shrep->[0], ( $win_size + $curr_gi ), $curr_gi, \@newel ) )
|
|
|
1544 if ($i_debug);
|
|
|
1545 } else {
|
|
|
1546 die "ERROR: the base-pairs in the EL don't match" . $curr_shrep . "\n";
|
|
|
1547 }
|
|
|
1548 }
|
|
|
1549
|
|
|
1550 # case (). CLOSE_EL
|
|
|
1551 if (@c_up) {
|
|
|
1552 unless ($cue) {
|
|
|
1553 if (@p_bp) {
|
|
|
1554 push( @c_el, @p_bp );
|
|
|
1555 } else {
|
|
|
1556 push( @c_el, ( $curr_gi, $curr_gi ) );
|
|
|
1557 }
|
|
|
1558
|
|
|
1559 push( @c_el, ( ( $curr_gi + $win_size - 1 ), ( $curr_gi + $win_size - 1 ) ) );
|
|
|
1560 my @newel = @c_el;
|
|
|
1561 push( @els, \@newel );
|
|
|
1562 @c_el = ();
|
|
|
1563 print STDERR "TEST - closed EL5 EOS:\n" . ( test( $curr_shrep->[0], ( $win_size + $curr_gi ), $curr_gi, \@newel ) )
|
|
|
1564 if ($i_debug);
|
|
|
1565 ## check case where cue is given, but window is completely unstructured
|
|
|
1566 } else {
|
|
|
1567 unless (@p_bp) {
|
|
|
1568 my @newel = ( $curr_gi, $curr_gi, ( $curr_gi + $win_size - 1 ), ( $curr_gi + $win_size - 1 ) );
|
|
|
1569 push( @els, \@newel );
|
|
|
1570 }
|
|
|
1571 }
|
|
|
1572 }
|
|
|
1573
|
|
|
1574 # # update current graph index (since all shrep indices have been read)
|
|
|
1575 $curr_gi += $win_size;
|
|
|
1576
|
|
|
1577 # ===================================================
|
|
|
1578 # build graph lines
|
|
|
1579 my @graph_lines = ();
|
|
|
1580
|
|
|
1581 # ===================================================
|
|
|
1582
|
|
|
1583 # graph line objects
|
|
|
1584 my @edg = ();
|
|
|
1585 my @stackgraph = ();
|
|
|
1586 my @abstrgraph = ();
|
|
|
1587
|
|
|
1588 #-------------------> stems,stacks
|
|
|
1589 # build basic graph edges, stacks, and vertices for stems in the abstract graph
|
|
|
1590 my @tmp = keys %stems;
|
|
|
1591 foreach my $sk (@tmp) {
|
|
|
1592 my @stem = @{ $stems{$sk} }; # stem (i,j,k,l)=> (i,j)= (0,1) and (k,l)=(2,3)
|
|
|
1593 die "ERROR: There are not enough elements in this stem: ", ( join( ",", @stem ) ), "\n" unless ( scalar(@stem) == 4 );
|
|
|
1594
|
|
|
1595 # infer base-pairs from stem
|
|
|
1596 @p_bp = ();
|
|
|
1597 @c_bp = ();
|
|
|
1598 my @stem_relation_idxs = ();
|
|
|
1599 my @stack_relation_idxs = ();
|
|
|
1600 my $stacksize = $stem[2] - $stem[0] + 1;
|
|
|
1601 for ( my $x = 0 ; $x < $stacksize ; $x++ ) {
|
|
|
1602
|
|
|
1603 # for (my $i = $stem[0]; $i <= $stem[2] ; $i++){
|
|
|
1604 # for(my $j = $stem[1]; $j >= $stem[3]; $j--){
|
|
|
1605 my $i = $stem[0] + $x;
|
|
|
1606 my $j = $stem[1] - $x;
|
|
|
1607 @p_bp = @c_bp;
|
|
|
1608 @c_bp = ( $i, $j );
|
|
|
1609 if ( defined $bps{"$i:$j"} ) {
|
|
|
1610
|
|
|
1611 # add indices to stem relation object
|
|
|
1612 push( @stem_relation_idxs, $i );
|
|
|
1613 push( @stem_relation_idxs, $j );
|
|
|
1614
|
|
|
1615 # add edge to graph edges
|
|
|
1616 push( @edg, "e $i $j s" );
|
|
|
1617
|
|
|
1618 # add stack if option given
|
|
|
1619 if ($stacks) {
|
|
|
1620 if (@p_bp) {
|
|
|
1621 push( @stackgraph, "v $curr_gi P" );
|
|
|
1622 push( @stackgraph, "e $curr_gi $p_bp[0] p" );
|
|
|
1623 push( @stackgraph, "e $curr_gi $p_bp[1] p" );
|
|
|
1624 push( @stackgraph, "e $c_bp[0] $curr_gi p" );
|
|
|
1625 push( @stackgraph, "e $c_bp[1] $curr_gi p" );
|
|
|
1626 push( @stack_relation_idxs, $curr_gi );
|
|
|
1627
|
|
|
1628 # update graph index
|
|
|
1629 ++$curr_gi;
|
|
|
1630 }
|
|
|
1631 }
|
|
|
1632 } else {
|
|
|
1633 die "ERROR: There is no such base-pair in the base-pair hash:" . " ($i, $j)\n";
|
|
|
1634 }
|
|
|
1635 }
|
|
|
1636
|
|
|
1637 # }
|
|
|
1638 # }
|
|
|
1639 # add stem to abstract graph
|
|
|
1640 # (graph indices become a bit jumbled, but saves time)
|
|
|
1641 push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#S" );
|
|
|
1642
|
|
|
1643 # stem object now saves the curr_gi index for that stem
|
|
|
1644 $stems{$sk} = $curr_gi;
|
|
|
1645 ++$curr_gi;
|
|
|
1646
|
|
|
1647 # add relations
|
|
|
1648 push( @abstrgraph, "v $curr_gi ^R" );
|
|
|
1649 push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
|
|
|
1650 foreach my $r (@stem_relation_idxs) {
|
|
|
1651 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1652 }
|
|
|
1653 if ($stacks) {
|
|
|
1654 foreach my $r (@stack_relation_idxs) {
|
|
|
1655 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1656 }
|
|
|
1657 }
|
|
|
1658 ++$curr_gi;
|
|
|
1659 }
|
|
|
1660
|
|
|
1661 #-------------------> add HLs
|
|
|
1662 foreach my $hl (@hls) {
|
|
|
1663
|
|
|
1664 # check size
|
|
|
1665 die "ERROR: the HL object is the incorrect size\n" unless ( @{$hl} == 2 );
|
|
|
1666
|
|
|
1667 # add vertex
|
|
|
1668 push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#HL" );
|
|
|
1669
|
|
|
1670 # add edge from stem to hl
|
|
|
1671 if ( defined( $stems{ $bps{"$hl->[0]:$hl->[1]"} } ) ) {
|
|
|
1672 push( @abstrgraph, "e " . $stems{ $bps{"$hl->[0]:$hl->[1]"} } . " $curr_gi n" );
|
|
|
1673 } else {
|
|
|
1674 die "ERROR: stem not defined in base-pair hash for $hl->[0]:$hl->[1]\n";
|
|
|
1675 }
|
|
|
1676 ++$curr_gi;
|
|
|
1677
|
|
|
1678 # add relations
|
|
|
1679 push( @abstrgraph, "v $curr_gi ^R" );
|
|
|
1680 push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
|
|
|
1681
|
|
|
1682 # add all nodes between i to j, inclusively
|
|
|
1683 for ( my $r = $hl->[0] ; $r <= $hl->[1] ; $r++ ) {
|
|
|
1684 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1685 }
|
|
|
1686 ++$curr_gi;
|
|
|
1687 }
|
|
|
1688
|
|
|
1689 #------------------> add BLs
|
|
|
1690 foreach my $bl (@bls) {
|
|
|
1691
|
|
|
1692 # check size
|
|
|
1693 die "ERROR: the BL object is the incorrect size\n" unless ( scalar( @{$bl} ) == 4 );
|
|
|
1694
|
|
|
1695 # add vertex
|
|
|
1696 push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#BL" );
|
|
|
1697
|
|
|
1698 # add edges from bl to both adjacent stems
|
|
|
1699 if ( defined( $stems{ $bps{"$bl->[0]:$bl->[1]"} } )
|
|
|
1700 && defined( $stems{ $bps{"$bl->[2]:$bl->[3]"} } ) ) {
|
|
|
1701
|
|
|
1702 # add outer stem: S->BL
|
|
|
1703 push( @abstrgraph, "e " . $stems{ $bps{"$bl->[0]:$bl->[1]"} } . " $curr_gi n" );
|
|
|
1704
|
|
|
1705 # add inner stem: BL->S
|
|
|
1706 push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$bl->[2]:$bl->[3]"} } . " n" );
|
|
|
1707 } else {
|
|
|
1708 die "ERROR: stems not defined for BL\n";
|
|
|
1709 }
|
|
|
1710 ++$curr_gi;
|
|
|
1711
|
|
|
1712 # add relations
|
|
|
1713 push( @abstrgraph, "v $curr_gi ^R" );
|
|
|
1714 push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
|
|
|
1715
|
|
|
1716 # add all nodes between i to k and l to j, inclusively
|
|
|
1717 for ( my $r = $bl->[0] ; $r <= $bl->[2] ; $r++ ) {
|
|
|
1718 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1719 }
|
|
|
1720 for ( my $r = $bl->[3] ; $r <= $bl->[1] ; $r++ ) {
|
|
|
1721 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1722 }
|
|
|
1723 ++$curr_gi;
|
|
|
1724 }
|
|
|
1725
|
|
|
1726 #------------------> add ILs
|
|
|
1727 foreach my $il (@ils) {
|
|
|
1728
|
|
|
1729 # check size
|
|
|
1730 die "ERROR: the IL object is the incorrect size\n" unless ( scalar( @{$il} ) == 4 );
|
|
|
1731
|
|
|
1732 # add vertex
|
|
|
1733 push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#IL" );
|
|
|
1734
|
|
|
1735 # add edges from il to both adjacent stems
|
|
|
1736 if ( defined( $stems{ $bps{"$il->[0]:$il->[1]"} } )
|
|
|
1737 && defined( $stems{ $bps{"$il->[2]:$il->[3]"} } ) ) {
|
|
|
1738
|
|
|
1739 # add outer stem: S->IL
|
|
|
1740 push( @abstrgraph, "e " . $stems{ $bps{"$il->[0]:$il->[1]"} } . " $curr_gi n" );
|
|
|
1741
|
|
|
1742 # add inner stem: IL->S
|
|
|
1743 push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$il->[2]:$il->[3]"} } . " n" );
|
|
|
1744 } else {
|
|
|
1745 die "ERROR: stems not defined for IL\n";
|
|
|
1746 }
|
|
|
1747 ++$curr_gi;
|
|
|
1748
|
|
|
1749 # add relations
|
|
|
1750 push( @abstrgraph, "v $curr_gi ^R" );
|
|
|
1751 push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
|
|
|
1752
|
|
|
1753 # add all nodes between i to k and l to j, inclusively
|
|
|
1754 for ( my $r = $il->[0] ; $r <= $il->[2] ; $r++ ) {
|
|
|
1755 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1756 }
|
|
|
1757 for ( my $r = $il->[3] ; $r <= $il->[1] ; $r++ ) {
|
|
|
1758 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1759 }
|
|
|
1760 ++$curr_gi;
|
|
|
1761 }
|
|
|
1762
|
|
|
1763 #------------------> add ELs
|
|
|
1764 foreach my $el (@els) {
|
|
|
1765
|
|
|
1766 #check size
|
|
|
1767 die "ERROR: the EL object is the incorrect size\n" unless ( scalar( @{$el} ) == 4 );
|
|
|
1768
|
|
|
1769 # add vertex
|
|
|
1770 push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#EL" );
|
|
|
1771
|
|
|
1772 # add edges from el to both adjacent stems, if available
|
|
|
1773 unless ( $el->[0] == $el->[1] ) {
|
|
|
1774 if ( defined( $stems{ $bps{"$el->[0]:$el->[1]"} } ) ) {
|
|
|
1775 push( @abstrgraph, "e " . $stems{ $bps{"$el->[0]:$el->[1]"} } . " $curr_gi n" );
|
|
|
1776 } else {
|
|
|
1777 die "ERROR: stems not defined for left BP of EL\n";
|
|
|
1778 }
|
|
|
1779 }
|
|
|
1780 unless ( $el->[2] == $el->[3] ) {
|
|
|
1781 if ( defined( $stems{ $bps{"$el->[2]:$el->[3]"} } ) ) {
|
|
|
1782 push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$el->[2]:$el->[3]"} } . " n" );
|
|
|
1783 } else {
|
|
|
1784 die "ERROR: stems not defined for right BP of EL\n";
|
|
|
1785 }
|
|
|
1786 }
|
|
|
1787 ++$curr_gi;
|
|
|
1788
|
|
|
1789 # add relations
|
|
|
1790 push( @abstrgraph, "v $curr_gi ^R" );
|
|
|
1791 push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
|
|
|
1792
|
|
|
1793 # @r relation edges
|
|
|
1794 for ( my $r = $el->[1] ; $r <= $el->[2] ; $r++ ) {
|
|
|
1795 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1796 }
|
|
|
1797
|
|
|
1798 # if BP, not EOS
|
|
|
1799 unless ( $el->[0] == $el->[1] ) {
|
|
|
1800 push( @abstrgraph, "e $curr_gi $el->[0] \@r" );
|
|
|
1801 }
|
|
|
1802
|
|
|
1803 # if BP, not EOS
|
|
|
1804 unless ( $el->[2] == $el->[3] ) {
|
|
|
1805 push( @abstrgraph, "e $curr_gi $el->[3] \@r" );
|
|
|
1806 }
|
|
|
1807 ++$curr_gi;
|
|
|
1808 }
|
|
|
1809
|
|
|
1810 #------------------> add MLs
|
|
|
1811 foreach my $ml (@mls) {
|
|
|
1812
|
|
|
1813 #check size
|
|
|
1814 die "ERROR: the ML object is not big enough!\n" unless ( scalar( @{$ml} ) >= 3 );
|
|
|
1815
|
|
|
1816 # add vertex
|
|
|
1817 push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#ML" );
|
|
|
1818
|
|
|
1819 # add edge from outer stem to first BP
|
|
|
1820 my $closing_bp = shift( @{$ml} );
|
|
|
1821 if ( defined( $stems{ $bps{"$closing_bp->[0]:$closing_bp->[1]"} } ) ) {
|
|
|
1822 push( @abstrgraph, "e " . $stems{ $bps{"$closing_bp->[0]:$closing_bp->[1]"} } . " $curr_gi n" );
|
|
|
1823 } else {
|
|
|
1824 die "ERROR: ML is not defined for BP $closing_bp->[0]:$closing_bp->[1]\n";
|
|
|
1825 }
|
|
|
1826
|
|
|
1827 # add remaining stems
|
|
|
1828 foreach my $bp ( @{$ml} ) {
|
|
|
1829 if ( defined( $stems{ $bps{"$bp->[0]:$bp->[1]"} } ) ) {
|
|
|
1830 push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$bp->[0]:$bp->[1]"} } . " n" );
|
|
|
1831 } else {
|
|
|
1832 die "ERROR: ML is not defined for BP $bp->[0]:$bp->[1]\n";
|
|
|
1833 }
|
|
|
1834 }
|
|
|
1835 ++$curr_gi;
|
|
|
1836
|
|
|
1837 # add relations
|
|
|
1838 push( @abstrgraph, "v $curr_gi ^R" );
|
|
|
1839 push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
|
|
|
1840
|
|
|
1841 # add @r relation edges between i and k, where (i,j) is the first BP and (k,l)
|
|
|
1842 # the second BP
|
|
|
1843 my $pbp_aref = $closing_bp;
|
|
|
1844 my $cbp_aref = shift( @{$ml} );
|
|
|
1845 for ( my $r = $pbp_aref->[0] ; $r <= $cbp_aref->[0] ; $r++ ) {
|
|
|
1846 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1847 }
|
|
|
1848
|
|
|
1849 # rest of @r relation edges
|
|
|
1850 foreach my $bp ( @{$ml} ) {
|
|
|
1851 $pbp_aref = $cbp_aref;
|
|
|
1852 $cbp_aref = $bp;
|
|
|
1853 for ( my $r = $pbp_aref->[1] ; $r <= $cbp_aref->[0] ; $r++ ) {
|
|
|
1854 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1855 }
|
|
|
1856 }
|
|
|
1857
|
|
|
1858 # add final stretch between last BP (m,n) and first closing base-pair (i,j)
|
|
|
1859 for ( my $r = $cbp_aref->[1] ; $r <= $closing_bp->[1] ; $r++ ) {
|
|
|
1860 push( @abstrgraph, "e $curr_gi $r \@r" );
|
|
|
1861 }
|
|
|
1862 ++$curr_gi;
|
|
|
1863 }
|
|
|
1864
|
|
|
1865 @graph_lines = ( @edg, @stackgraph, @abstrgraph );
|
|
|
1866 return ( \@graph_lines, $curr_gi );
|
|
|
1867 }
|
|
|
1868
|
|
|
1869 # just to test the current shrep dot-bracket structure parsing
|
|
|
1870 sub test {
|
|
|
1871 my ( $shrep, $idx, $gi, $mark_aref, $isML ) = @_;
|
|
|
1872
|
|
|
1873 my $result = "";
|
|
|
1874
|
|
|
1875 $result .= "$shrep\n";
|
|
|
1876 my @a = split( "", $shrep );
|
|
|
1877 for ( my $i = 0 ; $i <= @a ; $i++ ) {
|
|
|
1878 my $mark = 0;
|
|
|
1879 foreach my $j ( @{$mark_aref} ) {
|
|
|
1880 if ($isML) {
|
|
|
1881 $mark = 1 if ( $i == $j->[0] - $gi );
|
|
|
1882 $mark = 1 if ( scalar( @{$j} == 2 && $i == $j->[1] - $gi ) );
|
|
|
1883 } else {
|
|
|
1884 $mark = 1 if ( $i == $j - $gi );
|
|
|
1885 }
|
|
|
1886
|
|
|
1887 }
|
|
|
1888 if ( $i == $idx - $gi ) {
|
|
|
1889 $result .= "^";
|
|
|
1890 } elsif ($mark) {
|
|
|
1891 $result .= "'";
|
|
|
1892 } else {
|
|
|
1893 $result .= " ";
|
|
|
1894 }
|
|
|
1895 }
|
|
|
1896 $result .= "\n";
|
|
|
1897
|
|
|
1898 return $result;
|
|
|
1899 }
|
|
|
1900
|
|
|
1901 ###########################################################################
|
|
|
1902 # Extends or closes a multi-loop. This means a ML already exists.
|
|
|
1903 # We can have the case ((x)( or ((x).(
|
|
|
1904 #
|
|
|
1905 # Input:
|
|
|
1906 # curr_shrep The current shrep object, with shrep string at pos 0
|
|
|
1907 # c_ml_aref The current open MLs array reference
|
|
|
1908 # p_bp_aref The previous base-pair array reference
|
|
|
1909 # stmbrks_aref The stem-breaks array reference
|
|
|
1910 # open_blks_aref The open bracket blocks array reference
|
|
|
1911 # idx The current index
|
|
|
1912 # #curr_gi The current graph index
|
|
|
1913 # Output:
|
|
|
1914 # None, it should just modify the ML accordingly
|
|
|
1915 ############################################################################
|
|
|
1916 sub extend_or_open_nested_ML {
|
|
|
1917 my ( $curr_shrep, $c_mls_aref, $p_bp_aref, $up_aref, $stmbrks_aref, $open_blks_aref, $idx, $curr_gi ) = @_;
|
|
|
1918
|
|
|
1919 # case '(x(x)x(x)x(' extend current ML
|
|
|
1920 if ( $c_mls_aref->[-1]->[-1]->[-1] == $p_bp_aref->[0] ) {
|
|
|
1921 push( @{ $c_mls_aref->[-1]->[-1] }, $p_bp_aref->[1] ); #close previous BP
|
|
|
1922 my @a = ($idx);
|
|
|
1923 push( @{ $c_mls_aref->[-1] }, \@a ); # add new open base-pair
|
|
|
1924 print STDERR "TEST - extended ML:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, $c_mls_aref->[-1], 1 ) ) if ($i_debug);
|
|
|
1925
|
|
|
1926 # case '(x(x)x((x)(' nested ML, => OPEN_ML with stem-break
|
|
|
1927 } elsif ( @$stmbrks_aref
|
|
|
1928 && $stmbrks_aref->[-1]->[1] == $p_bp_aref->[0]
|
|
|
1929 && $stmbrks_aref->[-1]->[2] == $p_bp_aref->[1] ) {
|
|
|
1930
|
|
|
1931 my @c_stmbrk = @{ pop( @{$stmbrks_aref} ) };
|
|
|
1932 my @a = ( $c_stmbrk[0], -1 );
|
|
|
1933 my @ml = ();
|
|
|
1934 push( @ml, \@a ); # 1st base-pair is in-waiting
|
|
|
1935 my @a2 = ( $c_stmbrk[1], $c_stmbrk[2] );
|
|
|
1936 push( @ml, \@a2 ); #2nd BP
|
|
|
1937 my @a3 = ($idx);
|
|
|
1938 push( @ml, \@a3 ); #first of current BP
|
|
|
1939 push( @{$c_mls_aref}, \@ml );
|
|
|
1940 print STDERR "TEST - opened nested ML with SB:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
|
|
|
1941
|
|
|
1942 # case '(x(x)x(.(x)(' nested ML, => OPEN_ML with UP
|
|
|
1943 } elsif ( @$up_aref
|
|
|
1944 && $up_aref->[-1]->[0] - 1 == $open_blks_aref->[-1]->[-1]
|
|
|
1945 && $up_aref->[-1]->[1] + 1 == $p_bp_aref->[0] ) {
|
|
|
1946 my @tmp_up = @{ pop( @{$up_aref} ) };
|
|
|
1947 my @a = ( $open_blks_aref->[-1]->[-1], -1 );
|
|
|
1948 my @ml = ();
|
|
|
1949 push( @ml, \@a ); # 1st BP awaits closing
|
|
|
1950 my @newbp = @{$p_bp_aref};
|
|
|
1951 push( @ml, \@newbp ); #2nd BP
|
|
|
1952 my @a2 = ($idx);
|
|
|
1953 push( @ml, \@a2 ); # current open BP, awaits closing
|
|
|
1954 push( @{$c_mls_aref}, \@ml );
|
|
|
1955 print STDERR "TEST - opened nested ML with prev UP\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
|
|
|
1956 } else {
|
|
|
1957 die "ERROR: in the case extend or open nested ML, " . "but the base-pairs don't fit for either!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) ) if ($i_debug);
|
|
|
1958 }
|
|
|
1959 }
|
|
|
1960
|
|
|
1961 ###########################################################################
|
|
|
1962 # Closes a stem-object for the abstract graph. If necessary, it also
|
|
|
1963 # creates a stem-break because the current open block is still open,
|
|
|
1964 # which can only mean that there is a stem-break between two consecutive
|
|
|
1965 # opening brackets.
|
|
|
1966 #
|
|
|
1967 # Input:
|
|
|
1968 # c_stem_aref : current stem array reference
|
|
|
1969 # stems_href : all stems hash reference
|
|
|
1970 # p_bp_aref : previous base-pair array reference
|
|
|
1971 # c_cp_aref : current base-pair array reference
|
|
|
1972 # p_open_aref : previously opened open-bracket block
|
|
|
1973 # p_stmbrk_aref : the previous stem-break
|
|
|
1974 # open_blks_aref : all old open blocks array reference
|
|
|
1975 # bps_href : all base-pairs hash reference
|
|
|
1976 #
|
|
|
1977 # Output:
|
|
|
1978 # None, it should just modify the input variables accordingly
|
|
|
1979 ############################################################################
|
|
|
1980 sub close_stem {
|
|
|
1981 my ( $c_stem_aref, $stems_href, $p_bp_aref, $c_bp_aref, $p_open_aref, $stmbrks_aref, $open_blks_aref, $bps_href ) = @_;
|
|
|
1982
|
|
|
1983 my @newstem = @$c_stem_aref;
|
|
|
1984 my $stemkey = "$newstem[0]:$newstem[1],$newstem[2]:$newstem[3]";
|
|
|
1985 $stems_href->{$stemkey} = \@newstem;
|
|
|
1986
|
|
|
1987 # reset and clean running variables
|
|
|
1988 @$c_stem_aref = ();
|
|
|
1989 @$p_bp_aref = @$c_bp_aref;
|
|
|
1990 @$c_bp_aref = ();
|
|
|
1991
|
|
|
1992 # add stem to base-pair hash
|
|
|
1993 $bps_href->{"$newstem[0]:$newstem[1]"} = $stemkey;
|
|
|
1994 $bps_href->{"$newstem[2]:$newstem[3]"} = $stemkey;
|
|
|
1995
|
|
|
1996 # there are still bases open in the previous open block, '[[]x',
|
|
|
1997 # where [] is a closed stem and [ is an open block,
|
|
|
1998 # must create STEM-BREAK
|
|
|
1999 # this cannot be the case ((x)), because there p_open is empty
|
|
|
2000 if (@$p_open_aref) {
|
|
|
2001
|
|
|
2002 # check: the next open base is adjacent to the previous BP
|
|
|
2003 if ( $p_open_aref->[-1] == $p_bp_aref->[0] - 1 ) {
|
|
|
2004 my @stmbrk = ( $p_open_aref->[-1], @$p_bp_aref );
|
|
|
2005 push( @$stmbrks_aref, \@stmbrk );
|
|
|
2006
|
|
|
2007 # push the previous open block back onto stack
|
|
|
2008 my @newopen = @$p_open_aref;
|
|
|
2009 push( @$open_blks_aref, \@newopen );
|
|
|
2010 @$p_open_aref = ();
|
|
|
2011 } else {
|
|
|
2012
|
|
|
2013 # error, this should not be possible
|
|
|
2014 die("ERROR in close_stem(), [.[]x should not occur here\n");
|
|
|
2015 }
|
|
|
2016 }
|
|
|
2017 }
|
|
|
2018
|
|
|
2019 ###########################################################################
|
|
|
2020 # Here the information about the structure is added to the graph. In a
|
|
|
2021 # first step each base-pair is added by adding edges between the vertices
|
|
|
2022 # of the respective nucleotides. Then, if two base-pairs are stacked,
|
|
|
2023 # we add an extra vertex and connect the four nucleotide vertices that
|
|
|
2024 # are involved with extra edges (symbols P used for vertex and p for edges).
|
|
|
2025 # While adding new vertices we keep track of the current graph index.
|
|
|
2026 #
|
|
|
2027 # Input:
|
|
|
2028 # curr_shrep : the current shrep as dot-bracket structure
|
|
|
2029 # curr_gi : the current graph index
|
|
|
2030 # win_size : the window size
|
|
|
2031 # stacks : the input parameter -stacks (whether to add stacking information)
|
|
|
2032 #
|
|
|
2033 # Output:
|
|
|
2034 # The structure graph information line-by-line as an array reference and
|
|
|
2035 # the current graph index.
|
|
|
2036 ############################################################################
|
|
|
2037 sub getStructGraph {
|
|
|
2038 my ( $curr_shrep, $curr_gi, $win_size, $stacks ) = @_;
|
|
|
2039
|
|
|
2040 my @struct = split( "", $curr_shrep->[0] );
|
|
|
2041
|
|
|
2042 my @edg = ();
|
|
|
2043 my @starts = ();
|
|
|
2044 my @pairs = ();
|
|
|
2045
|
|
|
2046 foreach my $idx ( 0 .. @struct - 1 ) {
|
|
|
2047
|
|
|
2048 push( @starts, $idx ) if ( $struct[$idx] eq "(" );
|
|
|
2049
|
|
|
2050 if ( $struct[$idx] eq ")" ) {
|
|
|
2051 my $start = $curr_gi + pop(@starts);
|
|
|
2052 my $end = $curr_gi + $idx;
|
|
|
2053 push( @edg, "e " . $start . " " . $end . " s " );
|
|
|
2054 my @pair = ( $start, $end );
|
|
|
2055 push( @pairs, \@pair );
|
|
|
2056 }
|
|
|
2057 }
|
|
|
2058
|
|
|
2059 # update current graph index
|
|
|
2060 $curr_gi = $curr_gi + $win_size;
|
|
|
2061
|
|
|
2062 # initiate structure lines array
|
|
|
2063 my @stacking_info = ();
|
|
|
2064
|
|
|
2065 # add stacking information to graph unless input option tells us not to
|
|
|
2066 if ($stacks) {
|
|
|
2067
|
|
|
2068 my $stacked_pairs = 0;
|
|
|
2069
|
|
|
2070 # get stacked base-pairs (they are ordered according to
|
|
|
2071 # position of closing nucleotide)
|
|
|
2072 for ( my $i = 1 ; $i < @pairs ; $i++ ) {
|
|
|
2073
|
|
|
2074 # add stacked base-pairs (vertices+edges)
|
|
|
2075 my $curr_bp_aref = $pairs[$i];
|
|
|
2076 my $prev_bp_aref = $pairs[ $i - 1 ];
|
|
|
2077
|
|
|
2078 # if the current base-pair is stacked on the previous base-pair,
|
|
|
2079 # when curr_start = prev_start - 1 AND curr_end = prev_end + 1
|
|
|
2080 if ( $curr_bp_aref->[0] == $prev_bp_aref->[0] - 1
|
|
|
2081 && $curr_bp_aref->[1] == $prev_bp_aref->[1] + 1 ) {
|
|
|
2082
|
|
|
2083 # add stacking vertex P
|
|
|
2084 push( @stacking_info, "v $curr_gi P" );
|
|
|
2085
|
|
|
2086 # add four edges from involved nucleotids
|
|
|
2087 push( @stacking_info, "e $curr_gi $prev_bp_aref->[0] p" );
|
|
|
2088 push( @stacking_info, "e $curr_gi $prev_bp_aref->[1] p" );
|
|
|
2089 push( @stacking_info, "e $curr_bp_aref->[0] $curr_gi p" );
|
|
|
2090 push( @stacking_info, "e $curr_bp_aref->[1] $curr_gi p" );
|
|
|
2091 ++$curr_gi; # add one to the index, ready for next vertex
|
|
|
2092
|
|
|
2093 }
|
|
|
2094 }
|
|
|
2095 }
|
|
|
2096
|
|
|
2097 my @str_graphlines = ( @edg, @stacking_info );
|
|
|
2098 return ( \@str_graphlines, $curr_gi );
|
|
|
2099 }
|
|
|
2100
|
|
|
2101 ###########################################################################
|
|
|
2102 # This method gives us the overall graph header for this sequence. This
|
|
|
2103 # graph includes all window calculations and all shreps for all windows.
|
|
|
2104 # Input:
|
|
|
2105 # seq_id : the sequence ID
|
|
|
2106 # seq_head : the header information for the sequence
|
|
|
2107 # wins_aref : the windows sizes for folding
|
|
|
2108 # h_shift : input parameter -shift
|
|
|
2109 # h_e : input parameter -e
|
|
|
2110 # h_c : input parameter -c
|
|
|
2111 # h_t : input parameter -t
|
|
|
2112 # h_u : input parameter -u
|
|
|
2113 # h_r : input parameter -r
|
|
|
2114 # h_M : input parameter -M
|
|
|
2115 # h_crop_unpaired_ends : input parameter -cue
|
|
|
2116 #
|
|
|
2117 # Output:
|
|
|
2118 # The header line as a string
|
|
|
2119 ############################################################################
|
|
|
2120 sub getGraphHeader {
|
|
|
2121 my ( $seq_id, $seq_head, $wins_aref, $h_shift, $h_e, $h_c, $h_t, $h_u, $h_r, $h_M, $h_crop_unpaired_ends, $h_i, $h_sample_len, $h_q, $h_Tp, $h_seqlen, $h_noStr ) = @_;
|
|
|
2122
|
|
|
2123 my $ret;
|
|
|
2124
|
|
|
2125 $ret = "t # SEQID $seq_id ";
|
|
|
2126 $ret .= "$seq_head " if ( defined $seq_head );
|
|
|
2127 $ret .= "SEQLEN $h_seqlen ";
|
|
|
2128
|
|
|
2129 return $ret . "\n" if ($h_noStr);
|
|
|
2130
|
|
|
2131 $ret .= "MAXWINSHREPS $h_M ";
|
|
|
2132 $ret .= "CUE 1 " if ($h_crop_unpaired_ends);
|
|
|
2133 $ret .= "RNASHAPES -w " . ( join( ",", @{$wins_aref} ) ) . " ";
|
|
|
2134 $ret .= "-SHIFT%w $h_shift " if ( not $globalFolding and $h_shift );
|
|
|
2135 $ret .= "-SHIFT 1 " if ( not $globalFolding and not $h_shift );
|
|
|
2136 $ret .= "-e $h_e " if ( defined $h_e );
|
|
|
2137 $ret .= "-c $h_c " if ( defined $h_c );
|
|
|
2138 $ret .= "-t $h_t " if ( defined $h_t );
|
|
|
2139 $ret .= "-u 1 " if ($h_u);
|
|
|
2140 $ret .= "-r 1 " if ($h_r);
|
|
|
2141 $ret .= "-i $h_i SAMPLE_LEN $h_sample_len" if ( not $h_q and $h_i );
|
|
|
2142 $ret .= "-q 1 " if ($h_q);
|
|
|
2143 $ret .= "-T $h_Tp" if ( $h_q and defined $h_Tp );
|
|
|
2144 $ret .= "\n";
|
|
|
2145
|
|
|
2146 return $ret;
|
|
|
2147 }
|
|
|
2148
|
|
|
2149 ###########################################################################
|
|
|
2150 # This method gives us the window header for the subgraph that includes
|
|
|
2151 # all shreps for the current window.
|
|
|
2152 # Input:
|
|
|
2153 # graphHead : the header line for the entire sequence graph
|
|
|
2154 # win_size : the size of the current window
|
|
|
2155 # win_start : the starting position of the current window
|
|
|
2156 # win_end : the end position of the current window
|
|
|
2157 # win_center : the centre position of the current window
|
|
|
2158 # win_global_fold : (boolean) whether the actual number of nucleotides
|
|
|
2159 # is smaller than the given window size (end of sequence?)
|
|
|
2160 #
|
|
|
2161 # Output:
|
|
|
2162 # The header line as a string
|
|
|
2163 ############################################################################
|
|
|
2164 sub getWindowHeader {
|
|
|
2165 my ( $graphHead, $win_size, $win_start, $win_end, $win_center, $win_global_fold, $win_sample, $used_t ) = @_;
|
|
|
2166
|
|
|
2167 chomp $graphHead;
|
|
|
2168
|
|
|
2169 $graphHead =~ s/t # //;
|
|
|
2170
|
|
|
2171 my $ret = "w # $graphHead ";
|
|
|
2172 $ret .= "SHAPE_TYPE $used_t ";
|
|
|
2173 $ret .= "GLOBALFOLD $win_global_fold ";
|
|
|
2174 $ret .= "WSIZE $win_size ";
|
|
|
2175 $ret .= "WSTART $win_start ";
|
|
|
2176 $ret .= "WEND $win_end ";
|
|
|
2177 $ret .= "WCENT $win_center ";
|
|
|
2178 $ret .= "WSAMPLE $win_sample ";
|
|
|
2179 $ret .= "\n";
|
|
|
2180
|
|
|
2181 return $ret;
|
|
|
2182 }
|
|
|
2183
|
|
|
2184 ###########################################################################
|
|
|
2185 # This method gives us the shape header that includes the single shrep
|
|
|
2186 # connected components.
|
|
|
2187 # Input:
|
|
|
2188 # winHead : the header for the current window
|
|
|
2189 # shrep : information about the current shrep
|
|
|
2190 #
|
|
|
2191 # Output:
|
|
|
2192 # The header line as a string
|
|
|
2193 ############################################################################
|
|
|
2194 sub getShapeHeader {
|
|
|
2195 my ( $winHead, $shrep ) = @_;
|
|
|
2196
|
|
|
2197 chomp $winHead;
|
|
|
2198
|
|
|
2199 $winHead =~ s/w # //;
|
|
|
2200 $winHead =~ s/t # //;
|
|
|
2201
|
|
|
2202 my $ret = "s # $winHead ";
|
|
|
2203
|
|
|
2204 my @info = @{$shrep};
|
|
|
2205 $ret .= join( " ", @info[ 1 .. $#info ] );
|
|
|
2206 $ret .= "\n";
|
|
|
2207
|
|
|
2208 return $ret;
|
|
|
2209 }
|
|
|
2210
|
|
|
2211 ###########################################################################
|
|
|
2212 # This method gives us the sequence header that includes the unstructured
|
|
|
2213 # sequence connected components.
|
|
|
2214 # Input:
|
|
|
2215 # winHead : the header for the current window
|
|
|
2216 # shrep : information about the current shrep
|
|
|
2217 #
|
|
|
2218 # Output:
|
|
|
2219 # The header line as a string
|
|
|
2220 ############################################################################
|
|
|
2221 sub getSeqHeader {
|
|
|
2222 my ( $winHead, $shrep ) = @_;
|
|
|
2223
|
|
|
2224 chomp $winHead;
|
|
|
2225
|
|
|
2226 $winHead =~ s/w # //;
|
|
|
2227 $winHead =~ s/t # //;
|
|
|
2228
|
|
|
2229 my $ret = "u # $winHead ";
|
|
|
2230
|
|
|
2231 my @info = @{$shrep};
|
|
|
2232 $ret .= join( " ", @info[ 1 .. $#info ] );
|
|
|
2233 $ret .= "\n";
|
|
|
2234
|
|
|
2235 return $ret;
|
|
|
2236 }
|
|
|
2237
|
|
|
2238
|
|
|
2239 ##################################################################################
|
|
|
2240 # This method parses a fasta file and is useful if the header ID is not-unique!!
|
|
|
2241 # It returns the header lines in an array and the sequence lines in the same
|
|
|
2242 # order. It is then up to the user to extract the parts from the header that is
|
|
|
2243 # necessary for the script.
|
|
|
2244 # Furthermore, the method deals with multiple lines, and returns a single sequence
|
|
|
2245 # without line breaks.
|
|
|
2246 # Input:
|
|
|
2247 # file The name of the fasta file
|
|
|
2248 # Output:
|
|
|
2249 # (1) An array reference for each header line
|
|
|
2250 # (2) An array reference for each sequence in same order as the headers
|
|
|
2251 ##################################################################################
|
|
|
2252 sub read_fasta_with_nonunique_headers_meta {
|
|
|
2253 my ($file) = @_;
|
|
|
2254 my $FUNCTION = "read_fasta_file in Sequences.pm";
|
|
|
2255
|
|
|
2256 my $header = "";
|
|
|
2257 my $seqstring = "";
|
|
|
2258 my @headers = ();
|
|
|
2259 my @sequences = ();
|
|
|
2260 my @meta = ();
|
|
|
2261 my %seq_meta = ();
|
|
|
2262 open( IN_HANDLE, "<$file" ) || die "ERROR in $FUNCTION:\n" . "Couldn't open the following file in package Tool," . " sub read_fasta_file: $file\n";
|
|
|
2263
|
|
|
2264 while ( my $line = <IN_HANDLE> ) {
|
|
|
2265 chomp($line);
|
|
|
2266
|
|
|
2267 # header (can contain one space after > symbol)
|
|
|
2268 if ( $line =~ /^\>(.*)/ ) {
|
|
|
2269 if ($header) {
|
|
|
2270 $seqstring =~ s/\s*//g; ## do not allow spaces in sequence
|
|
|
2271 push( @headers, $header );
|
|
|
2272 push( @sequences, $seqstring );
|
|
|
2273 push( @meta, +{%seq_meta} ); ## anonymous hash reference
|
|
|
2274 #print keys %seq_meta;
|
|
|
2275 $seqstring = "";
|
|
|
2276 undef(%seq_meta);
|
|
|
2277 }
|
|
|
2278 $header = $1;
|
|
|
2279
|
|
|
2280 } elsif ( $line =~ /(.+)\s+(#\S+)$/ && $header ) {
|
|
|
2281
|
|
|
2282 if ( exists $seq_meta{$2} ) {
|
|
|
2283 $seq_meta{$2} .= $1;
|
|
|
2284 } else {
|
|
|
2285 $seq_meta{$2} = $1;
|
|
|
2286 }
|
|
|
2287
|
|
|
2288 } elsif ($header) {
|
|
|
2289 $seqstring .= $line
|
|
|
2290 }
|
|
|
2291 }
|
|
|
2292
|
|
|
2293 if ($header) {
|
|
|
2294 $seqstring =~ s/\s*//g; ## do not allow spaces in sequence
|
|
|
2295 push( @headers, $header );
|
|
|
2296 push( @sequences, $seqstring );
|
|
|
2297 push( @meta, +{%seq_meta} );
|
|
|
2298 $seqstring = "";
|
|
|
2299 %seq_meta = ();
|
|
|
2300 }
|
|
|
2301 return ( \@headers, \@sequences, \@meta );
|
|
|
2302 }
|
|
|
2303
|
|
|
2304
|
|
|
2305 1;
|
|
|
2306
|
|
|
2307 #############################################################################
|
|
|
2308 # Programming description
|
|
|
2309 #############################################################################
|
|
|
2310 # Substructure graphs for machine learning with Fabrizio
|
|
|
2311 # -------------------------------------------------------
|
|
|
2312 #
|
|
|
2313 # (1) Parameters (RNAshapes parameter):
|
|
|
2314 # - Window sizes [] (-w)
|
|
|
2315 # - window shift size (-W)
|
|
|
2316 # - calculate structure probabilities (-r)
|
|
|
2317 # - energy range kcal/mol (-e) OR energy relative percentage (%) to MFE (-c)
|
|
|
2318 # - shape type 1 to 5 (-t)
|
|
|
2319 # - ignore unstable substructures (-u)
|
|
|
2320 # - max shreps
|
|
|
2321 #
|
|
|
2322 #
|
|
|
2323 # (2) For each sequence, generate one graph/file that consists of all windows. The general format for one graph is as follows:
|
|
|
2324 #
|
|
|
2325 # t # seq_id parameters
|
|
|
2326 # v graph_index nt_type window_size window_centre abs_seq_pos
|
|
|
2327 # ...
|
|
|
2328 # e m m+1 > 1 (backbone)
|
|
|
2329 # ...
|
|
|
2330 # e base_i_graph_index base_j_graph_index s shrep_e shrep_p ...
|
|
|
2331 #
|
|
|
2332 # For each window (subgraph) we create a subgraph (of the subgraph) for each substructure.
|
|
|
2333 # We have a running index (gi) for each subgraph. All vertex and edge indices of the subgraph add
|
|
|
2334 # the running graph index to the actual window position. For example
|
|
|
2335 #
|
|
|
2336 #
|
|
|
2337 # Sequence: AAACC CUUUG GG
|
|
|
2338 # 01234 56789 01
|
|
|
2339 #
|
|
|
2340 # Window=10 substructure1 = (((...))). centre 5.5
|
|
|
2341 #
|
|
|
2342 # v 0 A 10 5.5 0
|
|
|
2343 # v 1 A 10 5.5 1
|
|
|
2344 # v 2 A 10 5.5 2
|
|
|
2345 # v 3 C 10 5.5 3
|
|
|
2346 # v 4 C 10 5.5 4
|
|
|
2347 # v 5 C 10 5.5 5
|
|
|
2348 # v 6 U 10 5.5 6
|
|
|
2349 # v 7 U 10 5.5 7
|
|
|
2350 # v 8 U 10 5.5 8
|
|
|
2351 # v 9 G 10 5.5 9
|
|
|
2352 # e 0 1 > 1
|
|
|
2353 # e 1 2 > 1
|
|
|
2354 # e 2 3 > 1
|
|
|
2355 # e 3 4 > 1
|
|
|
2356 # e 4 5 > 1
|
|
|
2357 # e 5 6 > 1
|
|
|
2358 # e 6 7 > 1
|
|
|
2359 # e 7 8 > 1
|
|
|
2360 # e 8 9 > 1
|
|
|
2361 # e 0 8 s -15.0 0.1223
|
|
|
2362 # e 1 7 s -15.0 0.1223
|
|
|
2363 # e 2 6 s -15.0 0.1223
|
|
|
2364 #
|
|
|
2365 # gi = 9+1 = 10
|
|
|
2366 #
|
|
|
2367 # Window = 10 substructure2 = .(((...))) centre 6.5
|
|
|
2368 # v 10 A 10 6.5 2
|
|
|
2369 # v 11 C 10 6.5 3
|
|
|
2370 # v 12 C 10 6.5 4
|
|
|
2371 # v 13 C 10 6.5 5
|
|
|
2372 # v 14 U 10 6.5 6
|
|
|
2373 # v 15 U 10 6.5 7
|
|
|
2374 # v 16 U 10 6.5 8
|
|
|
2375 # v 17 G 10 6.5 9
|
|
|
2376 # v 18 G 10 6.5 10
|
|
|
2377 # v 19 G 10 6.5 11
|
|
|
2378 # e 10 11 > 1
|
|
|
2379 # e 11 12 > 1
|
|
|
2380 # e 12 13 > 1
|
|
|
2381 # e 13 14 > 1
|
|
|
2382 # e 14 15 > 1
|
|
|
2383 # e 15 16 > 1
|
|
|
2384 # e 16 17 > 1
|
|
|
2385 # e 17 18 > 1
|
|
|
2386 # e 18 19 > 1
|
|
|
2387 # e 11 19 s -17.0 0.156
|
|
|
2388 # e 12 18 s -17.0 0.156
|
|
|
2389 # e 13 17 s -17.0 0.156
|
|
|
2390 #
|
|
|
2391 # gi = 19+1 = 20
|
|
|
2392 #
|
|
|
2393 #
|
|
|
2394 #
|
|
|
2395 # Write one perl script to create graphs for a set of sequences, fasta2shrep_gspan.pl.
|
|
|
2396 #
|
|
|
2397 # INPUT:
|
|
|
2398 # -f fasta file with all sequences to compute
|
|
|
2399 # parameters as above
|
|
|
2400 #
|
|
|
2401 # OUTPUT:
|
|
|
2402 # one file per sequence that contains graph called seq_id.gspan
|
|
|
2403 #
|
|
|
2404 # (1) for each window size call RNAshapes and write to a tmp file
|
|
|
2405 # (2) parse result of RNAshapes (catch RNAshapes error - sequence too long?) check for max shreps.
|
|
|
2406 # (3) convert RNAshapes result to subgraph -> write to file (readpipe) look at efficiency and errors
|
|
|
2407 # (4) repeat (1) to (3) for each sequence
|
|
|
2408
|
|
|
2409 # ABSTRACT STRUCTURE GRAPH
|
|
|
2410 # To each shrep graph (s #) add an abstract structure graph with special
|
|
|
2411 # abstract relations. For example, to add the abstract structure graph to
|
|
|
2412 # the previous shrep structure with gi numbers from 10 to 19, we first
|
|
|
2413 # identify the abstract shape, i.e. EL-S-HL and then create nodes (labelled
|
|
|
2414 # with a given name-space followed the actual node name and separated by a
|
|
|
2415 # hash) and edges for this graph as follows:
|
|
|
2416 # NOTE: at the moment we add the adjacent base-pairs to the loop definitions
|
|
|
2417 # but these can be removed if necessary in the future
|
|
|
2418 # v 20 abstruct#EL
|
|
|
2419 # v 21 abstruct#S
|
|
|
2420 # v 22 abstruct#HL
|
|
|
2421 # e 20 21 n
|
|
|
2422 # e 21 22 n
|
|
|
2423 # v 23 ^R
|
|
|
2424 # e 20 23 ^r
|
|
|
2425 # e 23 10 @r
|
|
|
2426 # e 23 11 @r
|
|
|
2427 # e 23 19 @r
|
|
|
2428 # v 24 ^R
|
|
|
2429 # e 21 24 ^r
|
|
|
2430 # e 24 11 @r
|
|
|
2431 # e 24 12 @r
|
|
|
2432 # e 24 13 @r
|
|
|
2433 # e 24 17 @r
|
|
|
2434 # e 24 18 @r
|
|
|
2435 # e 24 19 @r
|
|
|
2436 # v 25 ^R
|
|
|
2437 # e 22 25 ^r
|
|
|
2438 # e 25 13 @r
|
|
|
2439 # e 25 14 @r
|
|
|
2440 # e 25 15 @r
|
|
|
2441 # e 25 16 @r
|
|
|
2442 # e 25 17 @r
|
|
|
2443
|
|
|
2444 # gi = 25+1 = 26
|
|
|
2445
|
|
|
2446 # ANOTATION FILE
|
|
|
2447 # This is a file that labels a sequence region with a given annotation. In one
|
|
|
2448 # file we can have annotations within different name-spaces, for example target
|
|
|
2449 # sites predicted with different tools.
|
|
|
2450 # File format is a tab-delimited file with the following columns:
|
|
|
2451 # SEQID (same ID as in the fasta file header)
|
|
|
2452 # i (left positon of region)
|
|
|
2453 # j (right position of region, if one position i=j)
|
|
|
2454 # NAMESPACE#LABEL (give each annotation type one name-space and choose a label depending on the task)
|
|
|
2455
|
|
|
2456 # E.g we have 2 different miRNAs and we predict the target-sites with (1) IntaRNA and (2) TargetSearch,
|
|
|
2457 # then this could be IntaRNA#miR1, IntaRNA#miR2, TargetSearch#miR1, and TargetSearch#miR2.
|
|
|
2458 # All labels can be used more than once for one or more sequences.
|
|
|
2459 # all labels with the same namespace and the same sequence ID are grouped into
|
|
|
2460 # one abstract graph, according to the order of i.
|
|
|
2461 # Example
|
|
|
2462 # SEQID i j NAMESPACE#LABEL
|
|
|
2463 # s1 10 20 IntaRNA#miR1
|
|
|
2464 # s1 54 60 IntaRNA#miR2
|
|
|
2465 # s1 15 25 TargetSearch#miR1
|
|
|
2466 # s1 54 60 TargetSearch#miR2
|
|
|
2467 # s2 ...
|
|
|
2468
|
|
|
2469 # We create connected component per sequence ID per name-space,
|
|
|
2470 # so that all labels with
|
|
|
2471 # the same namespace are grouped together (for one sequence) as follows:
|
|
|
2472 # v 26 IntaRNA#miR1
|
|
|
2473 # v 27 IntaRNA#miR2
|
|
|
2474 # e 26 27 n
|
|
|
2475 # then create R nodes and ^r and @r relations as before (from i to j)
|
|
|
2476
|
|
|
2477 # If a sequence graph option is given, then these graphs (u #) have to
|
|
|
2478 # be labelled with this annotation
|
|
|
2479
|
