# HG changeset patch
# User bgruening
# Date 1383059419 14400
# Node ID b01beb1702902e13630f114998922a54170abbf6

Uploaded

diff -r 000000000000 -r b01beb170290 fasta2shrep_gspan.pl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fasta2shrep_gspan.pl	Tue Oct 29 11:10:19 2013 -0400
@@ -0,0 +1,2479 @@
+#!/usr/bin/perl
+#use feature ':5.10';
+use File::Basename;
+use lib dirname($0);    # search skript directory for modules
+use strict 'vars';
+use warnings;
+use Getopt::Long;
+use Pod::Usage;
+#use List::Util qw/ min max /;
+
+#use StructureLibrary::RNAtoolsConfig;
+#use StructureLibrary::Sequence;
+
+use Cwd qw(getcwd abs_path);
+use File::Temp qw(tempdir);
+use File::Copy;
+use POSIX qw(ceil);
+
+use FindBin;
+use lib "$FindBin::Bin";
+#use GraphClust;
+
+#use vars qw(%CONFIG);
+#*CONFIG = \%GraphClust::CONFIG;
+
+=head1 NAME
+
+
+fasta2shrep_gspan.pl -fasta mysequences.fasta -wins "50,100,150,200" -shift 5 -M 8
+
+=head1 SYNOPSIS
+
+Options:
+
+		HELP
+        -help   brief help message
+        -man    full documentation
+
+        COMPULSORY
+        -fasta	<STRING> e.g. "sequence.fasta"
+        		All sequences in fasta format.
+
+       	OPTIONS
+        -wins		[INTEGER] e.g. "50,100,200"
+        		A list of window sizes to use.
+        		If none are given (empty string ''), then the entire sequence is
+        		taken with no windows. Each window > 1 required!
+        -shift		<INTEGER> e.g. 20
+        		The shift of the window, relative to the window size given in
+        		percent. So you give which percent of the window size shall be
+        		used for the shift. Of course the shift is rounded down to the
+        		nearest whole number.
+        		Example 20 % of a window 150 would result in a step size of 30 nt.
+        		It is a relative parameter, as you can give different window sizes.
+        		If you do not give this parameter there is a default shift of 1 nt.
+        -cue		Crop unpaired ends.
+        		If you give this flag, then the unpaired ends of each
+        		single structure are ignored. E.g. the structure
+        		...(((...))).. becomes just (((...)))
+        -stack		Adds stacking information to graphs. This adds an additional
+        		vertex (type P) for each pair of stacked base-pairs and four edges
+        		(type p) from each of the involved bases to the new vertex.
+        -e		<FLOAT> e.g. 5.0
+        		Energy range in kcal/mol (RNAshapes)
+        		Use only one of -e and -c!
+        -c		<INTEGER> e.g. 10
+        		Relative energy range, i.e. percentage (%) of MFE energy (RNAshapes)
+        		Use only one of -e and -c!
+        -t		<INTEGER> [1-5] e.g. 3 OR "3=0,4=100,5=200"
+        		The shape type (RNAshapes). Default is 3.
+        		With the list format, the shape level can be changed for different window length
+        		"4=100" means that shape level 4 is used from length 100nt (window length)
+        		The first given length has to be 0! Not continuous given levels are allowed!
+        -M		<INTEGER> e.g. 10
+        		Max number of shreps that should be taken per window.
+        -u 		Ignore unstable structures (RNAshapes).
+        		This option filters out closed structures with positive free energy.
+        -r		Calculate structure probabilities for shreps (RNAshapes)
+        -i		<INT> e.g. 10
+        		Turn on structure sampling and gives number of sampling iterations.
+        		Default no sampling (i=0)
+        -sample-len	<INT> e.g. 100
+        		Only in sampling mode: Sampling is only used for seqs/windows >= given length,
+        		Default: sample all lengths (0), if -i > 0
+        -q		Turn on shape probabilities for RNAshapes, no sampling mode allowed
+        -Tp		<FLOAT> e.g 0.001
+        		Filter cutoff for shape probabilities, applied before -M filter!
+        -seq-graph-win	add for each window a graph which contains no structure
+        -seq-graph-t	add for each 't #' a graph which contains no structure
+        -seq-graph-alph change the alphabet of unstructured graphs
+        -annotate		<STRING> annotation.tab
+        				A file with annotations to be added as abstract graphs
+        				on the sequence leven (if given) and on the structure
+        				(SHREP) level. The format is has the following TAB-delimited
+        				columns: SEQID, START, END, NAMESPACE#LABEL.
+        				Labels with the same name-space and SEQID form connected
+        				components, which is a sequence of label vertices ordered
+        				by the START position in the sequence.
+        -abstr			Add abstract structure graphs to the single shrep graph
+        				instances.
+        -nostr			Calculate no structures, only add sequence information,
+        				if this is given, then -seq-graph-win AND/OR -seq-graph-t
+        				are required.
+        -match-shape    <SHAPE>
+                all seqs/windows will be constraint folded into that shape via 
+                RNAshapes (if structure is given in another way this struct will be kept), 
+                if this shape is not possible within given energy range, produce a 
+                specific t graph with only one vertex 'X'. By this the instance 
+                becomes very unsimilar to all other graphs (for knn)
+        -vp     enable graph computation with viewpoints:
+                svmsgdnspdk will center on those nucleotides that are given
+                via capital letters and ignore those given as lowercase letters
+        -tmp		<STRING> e.g. "/scratch/1/sita/tmp"
+        		A directory for writing temporary files
+        -o		<STRING> e.g. "ProjectX/MySequences/GSPAN/"
+        		Output directory for gspan files containing graphs.
+        -group		<INTEGER> e.g. 5
+                        Combine/group that number of input seqs into 1 gspan file
+                        output name is then '<INT>.group.gspan.bz2'
+        -sge            Use SGE cluster for each sequence separately
+        -sge-logdir     stdout directory for SGE call
+        -sge-errdir     sdterr directory for SGE call
+        -stdout         send graphs to stdout instead of writing to files
+        -ignore-header  don't write fasta id part after first space to gspan
+        -debug          additional debug output
+
+
+        DEFAULT VALUES
+        -wins	""
+        -shift	1 nt
+        -c		10
+        -t		3
+        -M		0 # selects all shreps
+        -tmp    "/var/tmp/fasta2shrep"
+        -o		"CURRENT_DIR/GSPAN/"
+
+        SGE mode
+        -sge-logdir "CURRENT_DIR/GSPAN/SGE_log"
+        -sge-errdir "CURRENT_DIR/GSPAN/SGE_log"
+        -task-id <NUM>
+
+=head1 DESCRIPTION
+
+=cut
+
+###############################################################################
+# end handler for temporary directory
+# adds an error handler that deletes the directory in case of error
+# SIGUSR{1/2} are sent by the sge prior to the uncatchable SIGKILL if the
+# option -notify was set
+###############################################################################
+$SIG{'INT'}  = 'end_handler';
+$SIG{'TERM'} = 'end_handler';
+$SIG{'ABRT'} = 'end_handler';
+$SIG{'USR1'} = 'end_handler';
+$SIG{'USR2'} = 'end_handler';
+
+sub end_handler {
+  print STDERR "signal ", $_[0], " caught, cleaning up temporary files\n";
+
+  # change into home directory. deletion of the temporary directory will
+  # fail if it is the current working directory
+  chdir();
+  File::Temp::cleanup();
+  die();
+}
+
+###############################################################################
+# PARSE COMMAND LINE OPTIONS
+###############################################################################
+
+# command line options
+my ( $i_help, $i_man, $i_debug, $i_fas, $i_wins, $i_shift, $i_crop_unpaired_ends,
+  $i_r, $i_e, $i_c, $i_t, $i_u, $i_M, $i_o, $i_sge, $i_jobid, $i_tmp,
+  $i_ignore_seq_header, $i_stacks, $i_stdout, $i_q, $i_T, $i_i,
+  $i_sample_min_length, $i_sge_logDir, $i_sge_errDir, $i_groupsize, $i_annotate,
+  $i_abstr, $i_no_structure, $i_vp, $i_matchShape, $i_rnashapes_binpath );
+
+my ( $i_add_seq_graph_win, $i_add_seq_graph_t, $i_change_seq_graph_alph );
+
+my $options = GetOptions(
+  "help"            => \$i_help,
+  "man"             => \$i_man,
+  "debug"           => \$i_debug,
+  "fasta=s"         => \$i_fas,
+  "wins=s"          => \$i_wins,
+  "shift=f"         => \$i_shift,
+  "cue"             => \$i_crop_unpaired_ends,
+  "stack"           => \$i_stacks,
+  "r"               => \$i_r,
+  "e=f"             => \$i_e,
+  "c=i"             => \$i_c,
+  "t=s"             => \$i_t,
+  "u"               => \$i_u,
+  "M=i"             => \$i_M,
+  "tmp=s"           => \$i_tmp,
+  "o=s"             => \$i_o,
+  "i=i"             => \$i_i,
+  "sample-len=i"    => \$i_sample_min_length,
+  "q"               => \$i_q,
+  "Tp=f"            => \$i_T,
+  "seq-graph-win"   => \$i_add_seq_graph_win,
+  "seq-graph-t"     => \$i_add_seq_graph_t,
+  "-seq-graph-alph" => \$i_change_seq_graph_alph,
+  "sge"             => \$i_sge,
+  "task-id=i"       => \$i_jobid,
+  "ignore-header"   => \$i_ignore_seq_header,
+  "stdout"          => \$i_stdout,
+  "sge-logdir=s"    => \$i_sge_logDir,
+  "sge-errdir=s"    => \$i_sge_errDir,
+  "group=i"         => \$i_groupsize,
+  "annotate=s"      => \$i_annotate,
+  "abstr"           => \$i_abstr,
+  "nostr"           => \$i_no_structure,
+  "vp"              => \$i_vp,
+  "match-shape=s"   => \$i_matchShape,
+  "rnashapes-bin=s" => \$i_rnashapes_binpath,
+);
+
+pod2usage( -exitstatus => 1, -verbose => 1 ) if $i_help;
+pod2usage( -exitstatus => 0, -verbose => 2 ) if $i_man;
+($options) or pod2usage(2);
+
+# check compulsory options
+($i_fas) or pod2usage("Error: the option -fasta is compulsory!\n");
+( -e $i_fas ) or pod2usage("Error: no such file - $i_fas!\n");
+$i_fas = abs_path($i_fas);
+
+# check other options and set default values
+pod2usage("Error: either -e OR -c can be given, but NOT both!\n") if ( $i_e && $i_c );
+( $i_e or $i_c ) or $i_c = 10;    # set -c 10, if neither -e or -c are given
+($i_M) or $i_M = 0;    # max number of shreps is 0 (=means take all computed)
+($i_i) or $i_i = 0;    # default no sampling else sampling iterations
+($i_sample_min_length) or $i_sample_min_length = 0;
+pod2usage("\nError: use --sample-len <INT> only with -i <1..INT> !\n") if ( not $i_i and $i_sample_min_length );
+($i_T) or $i_T = 0;
+($i_q) or $i_q = 0;
+pod2usage("\nError: Sampling (-i) not possible with shape probabilities (-q)!\n") if ( $i_i and $i_q );
+($i_add_seq_graph_win)     or $i_add_seq_graph_win     = 0;
+($i_add_seq_graph_t)       or $i_add_seq_graph_t       = 0;
+($i_change_seq_graph_alph) or $i_change_seq_graph_alph = 0;
+($i_no_structure)          or $i_no_structure          = 0;
+
+if ($i_change_seq_graph_alph) {
+  ($i_add_seq_graph_t) or ($i_add_seq_graph_win) or pod2usage( "Error: " .
+"When giving the parameter -seq-graph-alph, then either -seq-graph-t or -seq-graph-win" .
+      " must also be given!\n" );
+}
+( -e $i_annotate ) or pod2usage("Error: no such file - $i_annotate!\n")
+  if ($i_annotate);
+
+$i_add_seq_graph_t = 1 if ($i_no_structure);
+
+($i_t) or $i_t = 3;    # default abstraction type is 3
+my $change_shape_level = 0;
+my @level_lens = ( -1, -1, -1, -1, -1 ); ## array_idx-1=shape level, value=start length of this level
+
+if ( $i_t !~ /^\d+$/ ) {
+
+  my @t_minlens = split( ",", $i_t );    ## -t "3=100,4=200"
+
+  foreach my $idx ( 1 .. @t_minlens ) {
+
+    my $level = $t_minlens[ $idx - 1 ];    ##  $level = "3=100"
+
+    die "$level Wrong -t format! Example: -t 3=0,4=100,5=200\n" if ( $level !~ /^\d+\=\d+$/ );
+    my @lev_len = split( "=", $level );    ##  $level = "3=100"
+
+    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 );
+    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 );
+    die "Wrong -t format! Length >= 0 expected! Example: -t 3=0,4=100,5=200\n" if ( $lev_len[1] < 0 );
+    $change_shape_level = 1;
+    $level_lens[ $lev_len[0] - 1 ] = $lev_len[1];
+  }
+
+  ($i_debug) and print STDERR "i_t = $i_t - change $change_shape_level -  shape level lengths#" . join( ":", @level_lens ) . "#" . join( ":", @t_minlens ) . "#\n";
+}
+
+## checks for match-shape
+if ($i_matchShape){
+  die "Please provide correct match shape string like '[]'! Exit...\n\n"
+    if ($i_matchShape !~ /^[\[\]_]+$/);
+  ## RNAshapes prodices anyway only 1 structure, no suboptimal structs in match-shape folding
+  $i_M = 1;
+}
+
+my $CURRDIR = getcwd;
+
+# set up tmp directory
+# default tmp is /var/tmp, usually not NFS mounted!
+( defined $i_tmp ) or $i_tmp = '/var/tmp/';
+my $tmp_template = 'fasta2shrep-XXXXXX';
+
+# CLEANUP => 1 : automatically delete at exit
+$i_tmp = tempdir( $tmp_template, DIR => $i_tmp, CLEANUP => 1 );
+
+# create GSPAN directory when not printing to stdout
+if ( not $i_stdout ) {
+  if ($i_o) {
+    ( -e $i_o ) or system("mkdir -p $i_o");
+  } else {
+    system("mkdir -p GSPAN");
+    $i_o = $CURRDIR . "/GSPAN/";
+  }
+}
+
+###############################################################################
+# GLOBAL VARIABLES
+###############################################################################
+$i_rnashapes_binpath or $i_rnashapes_binpath = "RNAshapes";
+die "Please provide full RNAshapes path (e.g. /usr/bin/RNAshapes) !\n" if (!-f "$i_rnashapes_binpath");
+my $rnashapes_loc = $i_rnashapes_binpath;
+if ( !$rnashapes_loc || !-e $rnashapes_loc ) {
+  my $loc = `which RNAshapes`;
+  chomp($loc);
+  die "\nCannot find RNAshapes binary! Exit...\n\n" if ( !$loc );
+  die "\nCannot find RNAshapes binary! Exit...\n\n" if ( !-e $loc );
+  $rnashapes_loc = $loc;
+}
+
+my @WINDOWS = ();
+@WINDOWS = split( ",", $i_wins ) if ($i_wins);
+my $globalFolding;
+$globalFolding = 1 unless @WINDOWS;
+my $CURRUSER = getlogin;
+my $SEQNO    = 0;          # used to id the sequences
+my $GSPANNO  = 0;          # used for gspan filenames
+# minimum length of sequence that will work with RNAshapes
+# sequences with one or two nucleotides will be restricted to sequence-only graphs
+# this ensures that no sequences are skipped and external info kept synchronized
+my $GSPAN_SEQ_MINLEN = 3;
+
+# name spaces
+my $ABSTRUCT = "AS";
+###############################################################################
+# EXECUTION CODE
+###############################################################################
+
+# read fasta file into hash
+my ( $headers_aref, $sequences_aref, $metainfo_aref ) = read_fasta_with_nonunique_headers_meta($i_fas);
+
+# call script again on the sge cluster in a batch job
+if ($i_sge) {
+  my $SCRIPTNAME    = $0;
+  my $CLUSTERSUBMIT = $FindBin::Bin."/fasta2shrep_gspan.sge";
+
+  die "Cannot find SGE submit script $CLUSTERSUBMIT! Exit...\n\n" if ( !-e $CLUSTERSUBMIT );
+
+  my $SGE_jobs = @{$sequences_aref};
+  $SGE_jobs = ceil( @{$sequences_aref} / $i_groupsize ) if ($i_groupsize);
+
+  my $params = "-fasta $i_fas ";
+  $params .= "-wins $i_wins "                    if ($i_wins);
+  $params .= "-shift $i_shift "                  if ($i_shift);
+  $params .= "-e $i_e "                          if ($i_e);
+  $params .= "-c $i_c "                          if ($i_c);
+  $params .= "-t $i_t "                          if ($i_t);
+  $params .= "-u "                               if ($i_u);
+  $params .= "-r "                               if ($i_r);
+  $params .= "-M $i_M "                          if ($i_M);
+  $params .= "-o $i_o "                          if ($i_o);
+  $params .= "-i $i_i "                          if ($i_i);
+  $params .= "--sample-len $i_sample_min_length" if ($i_sample_min_length);
+  $params .= "-q "                               if ($i_q);
+  $params .= "-Tp $i_T "                         if ($i_T);
+  $params .= "--seq-grap-win "                   if ($i_add_seq_graph_win);
+  $params .= "--seq-graph-t "                    if ($i_add_seq_graph_t);
+  $params .= "--seq-graph-alph "                 if ($i_change_seq_graph_alph);
+  $params .= "-ignore-header "                   if ($i_ignore_seq_header);
+  $params .= "-cue "                             if ($i_crop_unpaired_ends);
+  $params .= "-stack "                           if ($i_stacks);
+  $params .= "--group $i_groupsize "             if ($i_groupsize);
+  $params .= "-annotate $i_annotate"             if ($i_annotate);
+  $params .= "-abstr"                            if ($i_abstr);
+  $params .= "-nostr "                           if ($i_no_structure);
+  $params .= "--debug "                          if ($i_debug);
+  $params .= "-vp "                              if ($i_vp);
+  $params .= "--match-shape $i_matchShape "      if ($i_matchShape);
+
+  print "used script:" . $SCRIPTNAME . "\n";
+  print "used submit script:" . $CLUSTERSUBMIT . "\n";
+
+  $i_sge_logDir = "$i_o/SGE_log" if ( !$i_sge_logDir );
+  mkdir($i_sge_logDir);
+
+  $i_sge_errDir = $i_sge_logDir if ( !$i_sge_errDir );
+  mkdir($i_sge_errDir);
+
+  my $ssh = 1;    ## can be used for debug shell script call
+  if ($ssh) {
+    system( "ssh $CURRUSER\@biui.informatik.uni-freiburg.de "
+        . "'export SGE_ROOT=/opt/sge-6.0/; cd $CURRDIR; "
+        . "/opt/sge-6.0/bin/lx24-amd64/qsub -t 1-$SGE_jobs -o $i_sge_errDir/ -e $i_sge_errDir/ "
+        . "$CLUSTERSUBMIT $CURRDIR $SCRIPTNAME \"$params\" ' " );
+  } else {
+    system("$CLUSTERSUBMIT $CURRDIR $SCRIPTNAME '$params' ");
+  }
+  exit;
+}
+
+## compute shreps and gspan, either local or after SGE submission
+
+# TODO read and process annotations
+
+my @used_seq_headers;
+my @used_seqs;
+my @used_meta;
+my $group_idx;
+
+if ($i_jobid) {
+  ## just process the one sequence as given by the jobid number
+  my $used_grouping = 1;   ## if no group is given, make 1 seq per job = group=1
+  $used_grouping = $i_groupsize if ($i_groupsize);
+  my $st  = ( $i_jobid - 1 ) * $used_grouping + 1;
+  my $end = $st + $used_grouping - 1;
+  $end = @{$sequences_aref} if ( $end > @{$sequences_aref} );
+
+  foreach my $idx ( $st .. $end ) {
+    push( @used_seq_headers, $headers_aref->[ $idx - 1 ] );
+    push( @used_seqs,        $sequences_aref->[ $idx - 1 ] );
+    push( @used_meta,        $metainfo_aref->[ $idx - 1 ] );
+  }
+
+  $group_idx = $i_jobid;
+
+  #($i_debug) and print STDERR "st $st end $end gr $group_idx job $i_jobid gs $i_groupsize\n";
+
+  $GSPANNO = $i_jobid-1 if (!$i_groupsize);
+
+} else {
+  ## process all sequences at once
+  @used_seq_headers = @{$headers_aref};
+  @used_seqs        = @{$sequences_aref};
+  @used_meta        = @{$metainfo_aref};
+}
+
+my $out;
+my $gspanfile;
+my $out_no_match_shape;
+
+if ($i_matchShape && !$i_stdout && !$i_groupsize){
+  open($out_no_match_shape,">$i_o/fasta2shrep.no_match");
+}
+
+# for each sequence in the fasta file
+while ( my $seq = shift @used_seqs ) {
+  my $tmp_header = shift @used_seq_headers;
+  my $tmp_meta   = shift @used_meta;
+
+  my ( $seq_id, $seq_header ) = ( $tmp_header =~ /(\S+)\s*([\S*\s*]*)/ );
+  $i_ignore_seq_header and $seq_header = '';
+
+  my $seq_fasta = generate_single_fasta_from_sequence_X( $seq_id, $seq );
+  my $seq_len = length($seq);
+
+  # only print sequence graphs for sequences below this threshold
+  my $no_structure_override = ($seq_len < $GSPAN_SEQ_MINLEN) ? 1 : 0;
+
+  $GSPANNO++;
+
+  # set outstream for gspan output to correct file/STDOUT
+  if ($i_stdout) {
+    $out = \*STDOUT;
+  } elsif ( !$i_groupsize ) {
+    $gspanfile = $i_tmp . '/' . $GSPANNO . '.gspan';
+    open( $out, "| bzip2 -f > $gspanfile.bz2" );
+  } elsif ( ( $GSPANNO - 1 ) % $i_groupsize == 0 ) {
+
+    if ( $GSPANNO > 1 ) {
+      close($out);
+      system("mv $gspanfile.bz2 $i_o/$group_idx.group.gspan.bz2");
+      if ($out_no_match_shape){
+        close($out_no_match_shape);
+        system("mv $gspanfile.no_match $i_o/$group_idx.group.gspan.no_match");
+      }
+    }
+
+    if ( !$i_jobid ) {
+      $group_idx = int( ( $GSPANNO - 1 ) / $i_groupsize ) + 1;
+    }
+
+    $gspanfile = "$i_tmp/$group_idx.group.gspan";
+    open( $out, "| bzip2 -f > $gspanfile.bz2" );
+    open( $out_no_match_shape, ">$gspanfile.no_match" ) if ($out_no_match_shape);
+  }
+
+  ## do not use folding windows in special cases
+  if ( $globalFolding || $i_no_structure || $no_structure_override) {
+    @WINDOWS = ();
+    push( @WINDOWS, $seq_len );
+  }
+
+  ##check win sizes for global folding
+  my @WINDOWS_used = sort { $a <=> $b } @WINDOWS;
+  foreach my $w_idx ( 0 .. ( @WINDOWS_used - 1 ) ) {
+    if ( $WINDOWS_used[$w_idx] >= $seq_len && $WINDOWS_used[$w_idx] > 1 ) {
+      if ( $w_idx < $#WINDOWS_used ) {
+        @WINDOWS_used = @WINDOWS_used[ 0 .. $w_idx ];
+        last;
+      }
+    }
+  }
+
+  ## use seq graph only if no shape folding wanted
+  $i_add_seq_graph_t = 1 if ($i_no_structure || $no_structure_override);
+
+  ## no shape info in graphheader if we have fixed structure (according LocaRNA handling)
+  ## use tags #FS or #S for provided structure
+  my $graph_header;
+  my @struct_meta = grep { $_ =~ /#FS/ || $_ =~ /#S/ } keys %{$tmp_meta};
+
+  if (@struct_meta) {
+    $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 );
+  } else {
+    $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 );
+  }
+
+  print $out $graph_header;
+
+  my $gi = 0;    # current graph index
+
+  ## add graph with no structure at all depending on $i_add_seq_graph_t
+  if ($i_add_seq_graph_t) {
+
+    ($gi) = convertSeqWindow( $seq, $seq_len, 1, $gi, $graph_header, $out, $i_annotate, $seq );
+
+  }
+
+  ## encode fixed structure only if provided and structures wanted in general
+  if ( @struct_meta && !$i_no_structure && !$no_structure_override) {
+
+    my $struct_meta = "#FS";
+    $struct_meta = "#S" if ( !exists $tmp_meta->{$struct_meta} );
+
+    my $seq_shrep = [ $tmp_meta->{$struct_meta}, "ENERGY", "0.00", "SHAPE", $struct_meta ];
+    $gi = convertShapeWindow( [$seq_shrep], $seq, $seq_len, 1, $gi, $out,
+      $graph_header, $i_annotate, $i_abstr, $i_crop_unpaired_ends, $i_stacks, $seq );
+    @WINDOWS_used = ();    ## no shape folding if we have a fixed structure
+
+  }
+
+  ## ignore RNAshapes folding if wanted, but do correct file move afterwards
+  ## (just "next" does not work due to output)
+  @WINDOWS_used = () if ($i_no_structure or $no_structure_override);
+
+  #for each window size in list
+  foreach my $win_size (@WINDOWS_used) {
+
+    # calculate shift size from percentage
+    my $curr_shift = 1;
+    if ($i_shift) {
+      $curr_shift = ( $i_shift / 100 ) * $win_size;
+      $curr_shift = int($curr_shift);                 #round down
+      $curr_shift = 1 unless ($curr_shift);           # just in case it is 0
+    }
+    ($i_debug) and print STDERR "winsize: $win_size curr_shift: $curr_shift\n";
+    ($i_debug) and print STDERR "\nNext: $seq_id\t winsize:$win_size \n";
+
+    # choose current shape level, depending on $i_t
+    my $curr_t = 0;
+    if ($change_shape_level) {
+      for ( my $i = 0 ; $i < @level_lens ; $i++ ) {
+        $curr_t = $i + 1 if ( $level_lens[$i] != -1 && ( $level_lens[$i] <= $win_size ) );
+      }
+      ($i_debug) and print STDERR "$win_size curr type $curr_t\n";
+    } else {
+      $curr_t = $i_t;
+    }
+
+    my $rnashapesoutput_fh;
+
+    # call RNAshapes and write to $rnashapesoutput_fh
+    $rnashapesoutput_fh = call_RNAshapes( $seq_fasta, $rnashapes_loc, $win_size,
+      $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 );
+
+    # read RNAshapes output from $rnashapesoutput_fh and write subgraph
+    # to gspan file
+    my $gi_old = $gi;
+    $gi = convert_RNAshapes_output( $rnashapesoutput_fh, $gi, $i_M, $out, $graph_header,
+      $win_size, $seq_len, $curr_t, $i_annotate, $i_abstr, $i_crop_unpaired_ends, $i_stacks, $seq );
+
+    ## no (match) shape found at all for this seq
+    if ($gi == $gi_old+1){
+      $gi = convertSeqWindow( "X", 1, 1, $gi, $graph_header, $out, $i_annotate, $seq );
+      print $out_no_match_shape $seq_id."\n" if (!$i_stdout && $out_no_match_shape);
+    }
+  }    ## foreach WINDOW_used
+
+  if ( !$i_stdout && !$i_groupsize ) {
+    close($out);
+    move "$gspanfile.bz2", "$i_o/$GSPANNO.gspan.bz2";
+  }
+
+  system("rm $seq_fasta");
+
+}    ## while @used_seqs
+
+if ($i_groupsize) {
+  close($out);
+  system("mv $gspanfile.bz2 $i_o/$group_idx.group.gspan.bz2");
+  close($out_no_match_shape) if ($out_no_match_shape);
+  move "$gspanfile.no_match", "$i_o/$group_idx.group.gspan.no_match" if ($out_no_match_shape);
+} elsif ( $out_no_match_shape && !$i_stdout && !$i_groupsize ) {
+  close($out_no_match_shape);
+}
+
+###############################################################################
+# METHODS
+###############################################################################
+
+############################################################################
+# Generates fasta file for a single sequence (one-lined-fasta). This
+# fasta is stored in the temp directory and should be deleted at the end.
+# Input:
+# seq_id : the sequence ID
+# seq : the sequence
+#
+# Output:
+# The fasta file name
+############################################################################
+sub generate_single_fasta_from_sequence_X {
+  my ( $seq_id, $seq ) = @_;
+
+  $seq = uc($seq);
+  $seq =~ tr/T/U/;
+  $seq =~ s/[^AUCGN]/N/g;
+
+  my $outfas = $i_tmp . "/seq_" . $SEQNO++ . ".fasta";
+  my $host   = readpipe("hostname");
+  open( FAS, ">$outfas" ) or die "$host Cannot open file $outfas! Exit...\n\n";
+  print FAS ">$seq_id\n$seq";
+  close(FAS);
+
+  return $outfas;
+}
+
+############################################################################
+# RNAshapes is called with the given input or default parameters.
+# Input:
+# seq_fasta : the sequence fasta file
+# rnashapes_location : the location of the installation files for RNAshapes
+# win_size : the current window size
+# shift : the input parameter -shift
+# e : the input parameter -e
+# c : the input parameter -c
+# t : the input parameter -t
+# u : the input parameter -u
+# r : the input parameter -r
+#
+# Output: none
+############################################################################
+sub call_RNAshapes {
+  my ( $seq_fasta, $rnashapes_location, $win_size, $shift, $e, $c, $t, $u, $r, $q, $T, $i, $sample_length, $seqLen, $matchShape ) = @_;
+  my $FUNCTION = "call_RNAshapes in fasta2shrep_gspan.pl";
+
+  ($seq_fasta) or die("INPUT ERROR in $FUNCTION: the fasta file is compulsory!\n");
+  ($rnashapes_location) or die( "INPUT ERROR in $FUNCTION: the RNAshapes location" . " is compulsory!\n" );
+  die "$rnashapes_location does not exist! Exit...\n\n" if ( !-e $rnashapes_location );
+  my $call = $rnashapes_location . " -o 1 ";    # the output format is of type 1
+  $call .= "-q " if ($q);
+  $call .= "-T $T " if ( $q and $T );
+  $call .= "-w $win_size ";
+  $call .= "-W $shift " if ($i_shift);
+  die("ERROR in $FUNCTION: Give only one of the options -c or -e (RNAshapes)!\n")
+
+    if ( $e && $c );
+  $call .= "-e $e " if ($e);
+  $call .= "-c $c " if ($c);
+  $call .= "-t $t " if ($t);
+  $call .= "-u "    if ( $u and not $i );       ## not possible in sampling mode
+  $call .= "-r "    if ($r);
+  $call .= "-m $matchShape " if ($matchShape);
+
+  ## 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
+  ## necessary to do sampling a large window is given, sample_length is shorter than window, but seq is longer than sample_len
+  $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
+  $call .= " < $seq_fasta";
+
+  ($i_debug) and print STDERR "$seqLen $sample_length $win_size $call\n";
+
+  open my $rnashapesoutput, "$call |" or die( "ERROR in $FUNCTION: The following call " . "could not be carried out!\n$call\n" );
+
+  return $rnashapesoutput;
+}
+
+############################################################################
+# The output of RNAshapes for one sequence and one window size is read
+# and converted into graph format.
+# Input:
+# rnashapeoutput : filehandle for RNAshapes output in format -o 1
+# curr_gi : current graph index (for vertices)
+# maxShreps : max number of shreps to convert to graphs
+# graph_file_hdl : the output handler for the graph file
+# graphHead : the header line for the complete graph (for sequence)
+# winSize : current window size in input for RNAshapes
+# seqLen  : full input seq length
+# used_t  :
+# annotate:
+# abstr   :
+# cue     :
+# stacks  :
+# orig_seq : the nucleotide sequence as read from fasta
+#
+# Output:
+# The current graph index
+############################################################################
+sub convert_RNAshapes_output {
+  my ( $rnashapesoutput, $curr_gi, $maxShreps, $graph_file_hdl, $graphHead, $winSize, $seqLen, $used_t, $annotate, $abstr, $cue, $stacks, $orig_seq ) = @_;
+
+  ## omit first line in output, contains fasta header line
+  my $line       = <$rnashapesoutput>;
+  my $win_shreps = [];
+  my $win_start;
+  my $win_end;
+  my $win_seq;
+  my $win_shrep_count = 0;
+  my $win_sample      = 0;
+  my $winHead;
+  my $win_globalFolding;
+  my $win_size_real;
+
+  # reading RNAshapes output
+  while ( $line = <$rnashapesoutput> ) {
+
+    if ( $line =~ /^(\d+)\s+(\d+)$/ ) {
+      ## line: "<start>    <end>"
+
+      if ( @{$win_shreps} > 0 ) {
+
+        print $graph_file_hdl $winHead;
+
+        ## remove SHAPE="_" shrep if it exists in $win_shreps
+        ## do this only when have a sequence graph already
+        if ($i_add_seq_graph_t) {
+          my @new_win_shreps = ();
+          map { push( @new_win_shreps, $_ ) if ( $_->[ $#{$_} ] ne "_" ) } @{$win_shreps};
+          $win_shreps = \@new_win_shreps;
+        }
+
+        ## add graph with no structure depending on $i_add_seq_graph_win to win
+        if ($i_add_seq_graph_win ) {
+          ( $curr_gi ) = convertSeqWindow( $win_seq, $win_size_real,
+            $win_start, $curr_gi, $winHead, $graph_file_hdl, $annotate, $orig_seq );
+        }
+
+        $curr_gi = convertShapeWindow( $win_shreps, $win_seq, $win_size_real,
+          $win_start, $curr_gi, $graph_file_hdl, $winHead, $annotate, $abstr,
+          $cue, $stacks, $orig_seq );
+      }
+
+      ## set new window params
+      $win_shreps      = [];
+      $win_start       = $1;
+      $win_end         = $2;
+      $win_shrep_count = 0;
+      $win_size_real   = $win_end - $win_start + 1;
+
+      if ( ($win_size_real) >= $seqLen ) {
+        $win_globalFolding = 1;
+      } else {
+        $win_globalFolding = 0;
+      }
+      my $win_center = $win_start + ( ( $win_size_real + 1 ) / 2 );
+
+      $winHead = getWindowHeader( $graphHead, $winSize, $win_start, $win_end, $win_center, $win_globalFolding, $win_sample, $used_t );
+
+    } elsif ( $line =~ /^(\S+)$/ ) {
+      ## line: "CUUAUGAGUAAGGAAAAUAACGAUUCGGGGUGACGCCCGAAUCCUCACUG"
+      $win_seq = uc($1); ## to be 101% sure that we have by default uppercase chars
+
+    } elsif ( $line =~ /^Results for (\d+) iterations:$/ ) {
+      ## line: "Results for 10 iterations:"
+      $win_sample = $1;
+
+    } elsif ( $line =~ /^Shape\s+\S+\s+not found within energy range.*$/ ) {
+      ## line: "Shape [] not found within energy range (-24.75 to -27.50). Try -c or -e to increase range."
+      $win_shreps      = [];
+
+    } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+(\S+)$/ ) {
+      ## line:"...((((..(((....)))))))...........(((((......)))))  (-10.10)  [[]][]"
+      ## take only $maxShreps shreps per window if set
+      next if ( $maxShreps && $win_shrep_count >= $maxShreps );
+      push( @{$win_shreps}, [ $1, "ENERGY", $2, "SHAPE", $3 ] );
+      $win_shrep_count++;
+
+    } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+\((\S+)\)\s+(\S+)$/ ) {
+      ## line:"((((..((((...((.((.((.....)).)).))...))))..))))...  (-10.60)  (0.7795360)  [[[[[]]]]]"
+      ## take only $maxShreps shreps per window if set
+      next if ( $maxShreps && $win_shrep_count >= $maxShreps );
+      push( @{$win_shreps}, [ $1, "ENERGY", $2, "PROB", $3, "SHAPE", $4 ] );
+      $win_shrep_count++;
+
+    } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+(\S+)\s+(\S+)$/ ) {
+      ## line:"((((..((((...((.((.((.....)).)).))...))))..))))...  (-10.60) 0.3000000 [[[[[]]]]]"
+      ## take only $maxShreps shreps per window if set
+      next if ( $maxShreps && $win_shrep_count >= $maxShreps );
+      push( @{$win_shreps}, [ $1, "ENERGY", $2, "SHAPEPROB", $3, "SHAPE", $4 ] );
+      $win_shrep_count++;
+
+    } elsif ( $line =~ /^([\(\)\.]+)\s+\((\S+)\)\s+\((\S+)\)\s+(\S+)\s+(\S+)$/ ) {
+      ## line:"((((..((((...((.((.((.....)).)).))...))))..))))...  (-10.60)  (0.7795360) 0.3000000 [[[[[]]]]]"
+      ## take only $maxShreps shreps per window if set
+      next if ( $maxShreps && $win_shrep_count >= $maxShreps );
+      push( @{$win_shreps}, [ $1, "ENERGY", $2, "PROB", $3, "SHAPEPROB", $4, "SHAPE", $5 ] );
+      $win_shrep_count++;
+
+    } else {
+      next if ( $line =~ /^$/ );
+      die "Unexpected shape output format!\nline=$line\n\nExit...\n\n";
+    }
+  }
+
+  ## convert last windows
+  if ( @{$win_shreps} > 0 ) {
+
+    print $graph_file_hdl $winHead;
+
+    ## remove SHAPE="_" shrep if it exists in $win_shreps
+    ## do this only when have a sequence graph already
+    if ($i_add_seq_graph_t) {
+      my @new_win_shreps = ();
+      map { push( @new_win_shreps, $_ ) if ( $_->[ $#{$_} ] ne "_" ) } @{$win_shreps};
+      $win_shreps = \@new_win_shreps;
+    }
+
+    ## add graph with no structure depending on $i_add_seq_graph_win to win
+    if ($i_add_seq_graph_win) {
+      ($curr_gi) = convertSeqWindow( $win_seq, $win_size_real, $win_start,
+        $curr_gi, $winHead, $graph_file_hdl, $annotate, $orig_seq );
+    }
+
+    $curr_gi = convertShapeWindow( $win_shreps, $win_seq, $win_size_real,
+      $win_start, $curr_gi, $graph_file_hdl, $winHead, $annotate, $abstr, $cue,
+      $stacks, $orig_seq );
+  }
+
+  close($rnashapesoutput);
+
+  return $curr_gi + 1;    # return the gi (graph index) for the next subgraph
+}
+
+## Sub to create graph for a complete unstructured sequence
+# TODO: document function
+sub convertSeqWindow {
+  my ( $win_seq, $win_size_real, $win_start, $curr_gi, $winHead, $graph_file_hdl, $annotate, $orig_seq ) = @_;
+
+  my $seq_shrep;
+
+  # use a different alphabet (lowercase) if option seq-graph-alph set
+  my $seq_graph_sequence = $i_change_seq_graph_alph ? lc($win_seq) : uc($win_seq);
+  $seq_shrep = [ "." x $win_size_real, "ENERGY", "0.00", "SHAPE", "_", "STRUCT", "." x $win_size_real, "SEQ", $seq_graph_sequence ];
+  my $backboneGraph_ref = getBackboneGraph( $seq_graph_sequence, $curr_gi, $win_start, 0, ( $win_size_real - 1 ), $orig_seq );
+
+  print $graph_file_hdl getSeqHeader( $winHead, $seq_shrep );
+  print $graph_file_hdl join( "\n", @{$backboneGraph_ref} ) . "\n";
+
+  $curr_gi += $win_size_real;
+
+  if ($annotate) {
+    ## TODO add annotation graphs
+  }
+
+  return ($curr_gi);
+}
+
+############################################################################
+# The results for one window are converted in this method to GSPAN graphs.
+# Input:
+# win_shreps_aref : the array ref of shreps (dot-bracket structures) for the
+# 					current window
+# win_seq : the nucleotide sequence for the current window
+# win_size_real : the current (true) window size
+# win_start : the starting position of win_seq in the original seq (1-n)
+# curr_gi : the current graph index
+# graph_file_hdl : the graph file handler
+# winHead   : the header line for the current sequence window
+# TODO annotate
+# abstr : input parameter -abstr
+# cue   : input parameter -cue
+# stacks : input parameter -stacks
+# orig_seq : the nucleotide sequence as read from fasta
+#
+# Output:
+# The current graph index
+############################################################################
+sub convertShapeWindow {
+  my ( $win_shreps_aref, $win_seq, $win_size_real, $win_start, $curr_gi,
+    $graph_file_hdl, $winHead, $annotate, $abstr, $cue, $stacks, $orig_seq ) = @_;
+
+  ## generate for each shrep a connected component in gspan file
+  foreach my $shrep ( @{$win_shreps_aref} ) {
+
+    # get the current gi as it was at the beginning of this shrep
+    my $shrep_gi = $curr_gi;
+
+    # cut off unpaired ends, if option is given
+    my $shrep_struct     = $shrep->[0];
+    my $crop_index_left  = 0;
+    my $crop_index_right = length($shrep_struct) - 1;
+
+    # find croping indices with numbering from 0 to n-1
+    if ($cue) {
+      $crop_index_left = index( $shrep_struct, "(" );    # find 1st occ of "("
+      $crop_index_right = rindex( $shrep_struct, ")" );  # find last occ of ")"
+
+      # if the complete window is unpaired, then don't crop
+      if ( $crop_index_left == -1 ) {
+        $crop_index_left  = 0;
+        $crop_index_right = length($shrep_struct) - 1;
+      }
+    }
+
+    # create structure graph
+    my $backboneGraph_ref = getBackboneGraph( $win_seq, $curr_gi, $win_start, $crop_index_left, $crop_index_right, $orig_seq );
+    my $structGraph_ref;
+
+    # add both basic edges, stacks and abstract graph
+    if ($abstr) {
+      ( $structGraph_ref, $curr_gi ) = getStructPlusAbstractStructGraph( $shrep, $curr_gi, $win_size_real, $cue, $stacks );
+
+      # just add basic structure graph edges, and stacks
+    } else {
+      ( $structGraph_ref, $curr_gi ) = getStructGraph( $shrep, $curr_gi, $win_size_real, $stacks );
+    }
+
+    # additional information for shrep header: sequence and dot-bracket
+    my $crop_length = $crop_index_right - $crop_index_left + 1;
+    my $shrepheader_struct = substr( $shrep_struct, $crop_index_left, $crop_length );
+    my $shrepheader_seq = substr( $win_seq, $crop_index_left, $crop_length );
+    push( @{$shrep}, 'STRUCT', $shrepheader_struct, 'SEQ', $shrepheader_seq );
+
+    # print structure graph to file
+    print $graph_file_hdl getShapeHeader( $winHead, $shrep );
+    print $graph_file_hdl join( "\n", @{$backboneGraph_ref} ) . "\n";
+
+    # don't print empty array; sgdnspdk breaks with empty lines
+    if ( @{$structGraph_ref} > 0 ) {
+      print $graph_file_hdl join( "\n", @{$structGraph_ref} ) . "\n";
+
+      if ($annotate) {
+
+        # TODO add annotations to shrep structure
+        # $win_start, $win_end (1-n)
+      }
+
+    }
+  }
+
+  return $curr_gi;
+}
+
+###########################################################################
+# Here the vertices for the nucleotides are generated and also the backbone
+# edges, which connect the nucleotides in their correct order (i.e order
+# of the given sequence)
+# Input:
+# win_seq : the sequence for the current window
+# win_start : the starting position for the current window
+# curr_crop_i_left : the cropping index for the left end of sequence (-cue)
+# curr_crop_i_right : the cropping index for the right end of sequence (-cue)
+# orig_seq : the nucleotide sequence as read from fasta
+#
+# Output:
+# The lines of the graph in an array reference to be printed to a file with
+# one element per line.
+############################################################################
+sub getBackboneGraph {
+  my ( $win_seq, $curr_gi, $win_start, $curr_crop_i_left, $curr_crop_i_right, $orig_seq ) = @_;
+
+  # RNAshapes substitutes T -> U, sequence only does not
+  # thus, just transform all t/T -> u/U to remain consistent
+  $orig_seq =~ tr /tT/uU/;
+  $win_seq  =~ tr /tT/uU/;
+
+  my @seq;
+  @seq = split( "", $win_seq );
+
+  # if the vp option is set we need to obtain the original capitalization from
+  # $orig_seq and extract the windowed sequence manually
+  my @seq_vp;
+  if ( defined $i_vp ) {
+    my $win_len = length($win_seq);
+    my $capitalized_win_seq = substr( $orig_seq, $win_start - 1, $win_len );
+    ( uc($capitalized_win_seq) eq uc($win_seq) ) or
+      die( "error: windowed sequence generated due to vp option not equal " .
+        "to sequence reported by RNASHAPES.\n" .
+        "'${capitalized_win_seq}' != '${win_seq}'" );
+    @seq_vp = split( "", $capitalized_win_seq );
+  }
+
+  my @vert = ();
+  my @edg  = ();
+
+  my $curr_abs_pos = $win_start + $curr_crop_i_left;
+  $curr_gi += $curr_crop_i_left;
+
+  # set vertice labeled with 'v' or 'V' according to sequence and vp option
+  # when vp set,
+  # uppercase nucleotides are annotated with 'v'
+  # lowercase nucleotides are annotated with 'V'
+  my $vertice_label = ( defined $i_vp and ( $seq_vp[$curr_crop_i_left] =~ /[a-z]/ ) ) ? 'V' : 'v';
+
+  # create backbone vertice of first nucleotide
+  push( @vert, join( ' ', $vertice_label, $curr_gi, $seq[$curr_crop_i_left], $curr_abs_pos ) );
+
+  foreach my $idx ( ( $curr_crop_i_left + 1 ) .. $curr_crop_i_right ) {
+    $curr_abs_pos++;
+    $curr_gi++;
+
+    # set vertice label as described above
+    $vertice_label = ( defined $i_vp and ( $seq_vp[$idx] =~ /[a-z]/ ) ) ? 'V' : 'v';
+    push( @vert, join( ' ', $vertice_label, $curr_gi, $seq[$idx], $curr_abs_pos ) );
+    push( @edg, join( ' ', "e", $curr_gi - 1, $curr_gi, '> 1' ) );
+  }
+
+  my @ret = ( @vert, @edg );
+
+  return \@ret;
+}
+
+###########################################################################
+# This method does the same as getStructGraph, but is extended to identify
+# the abstract parts of the structure and also add this to the graph
+# via abstract relations, see information below.
+# We already have the backbone graph, now the base-pair edges need to be
+# added to this graph and the abstract graph is added after this. Finally
+# the abstract relations, pointing from the abstract level to the basic
+# structure level is added. All the while, we keep track of the current
+# index graph. Furthermore, in the basic structure graph, we have vertices
+# stacks that connects the four stacked base-pairs involved. (symbols P
+# for vertex and p for edges). The abstract vertices are labelled according
+# to the given name-space, then #, then the abstract structure type, e.g.
+# HL for hairpin-loop. The edges between the structure types are denoted
+# with n, and the relation vertex is ^R, with ^r going from the abstract
+# structure to the relation and from the relation to the basic structure
+# is @r.
+#
+# Input:
+# curr_shrep : the current shrep as dot-bracket structure
+# curr_gi : the current graph index
+# win_size : the window size
+# cue : the input parameter -cue (whether to crop unpaired ends)
+# stacks : the input parameter -stacks (whether to add stack information)
+#
+# Output:
+# The structure graph information line-by-line as an array reference and
+# the current graph index.
+#
+############################################################################
+sub getStructPlusAbstractStructGraph {
+  my ( $curr_shrep, $curr_gi, $win_size, $cue, $stacks ) = @_;
+
+  #  $win_size =  length($curr_shrep->[0]);
+
+  my @struct = split( "", $curr_shrep->[0] );
+
+#  print "=================testing: new shrep : gi=$curr_gi ===============================\n" if ($i_debug);
+
+  # OBJECTS and VARIABLES
+  # all indices are saved according to the current graph index, $curr_gi,
+  # which is at the beginning of the current window, so that the nucleotide
+  # index can be inferred using $idx + $curr_gi
+  # opening brackets
+  my @open_blks = (); # open blocks of consecutive opening brackets (array of arrays)
+  my @p_open = ();  # currently open block of the last identified complete block
+  my @c_open = ();  # current "incomplete" block of consecutive opening brackets
+                    # base-pair objects
+  my %bps    = ();  # hash of all base-pairs with stem for a value, if it is at
+                    # the BP is opening or closing a stem object
+  my @c_bp   = ();  # current BP, just being closed at current index pos (i,j)
+  my @p_bp   = ();  # previously closed BP at position before current index
+                    # stem objects
+  my %stems  = ();  # all stems, key="i:j,k:l", value=stem array, (i,j) outer BP
+  my @c_stem = ();  # current incomplete stem object
+  my @stmbrks = (); # all stem-breaks in the form of (i,k,l), where i is the remaining
+                    # open base and (k,l) is the subsequent closing base-pair
+
+#  my @p_stmbrk	= (); # previous stem-break in the form of (i,k,l), where i is the remaining
+#  					  # open base and (k,l) is the subsequent closing base-pair
+# loop objects
+  my @up    = (); # all unpaired regions (i,j) that are of unkown type
+  my @c_up  = (); # current incomplete unpaired object
+  my @hls   = (); # all hairpin loop objects
+  my @bls   = (); # all bulge-loop objects
+  my @ils   = (); # all internal loops
+  my @mls   = (); # all multi-loops
+  my @c_mls = (); # current incomplete multi-loops! (there can be more than one)
+  my @els   = (); # all external loops
+  my @c_el  = (); # current incomplete external loop
+
+  # iterate through shrep structure an identify abstract parts
+  foreach my $idx ( 0 .. @struct - 1 ) {
+
+    #===================================================
+    # current char is the opening bracket of a base-pair
+    if ( $struct[$idx] eq "(" ) {
+
+      #===================================================
+
+      # update current window index to current graph index
+      $idx += $curr_gi;
+
+      # case: '((' currently there is an open block, extend it
+      if (@c_open) {
+        push( @c_open, $idx );
+
+# case: '.(' or ')('  or BOS'(' there is no open block, create a new one ".(" or ")("
+      } else {
+
+        # begin a new current open block
+        push( @c_open, $idx );
+
+        # case: ')('
+        if (@c_stem) {
+
+          # CLOSE STEM
+          print STDERR "TEST - closed stem:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_stem ) ) if ($i_debug);
+          close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps );
+
+          # case (x(x)(
+          if (@open_blks) {
+
+            # EXTEND_ML, if one exists
+            if (@c_mls) {
+              extend_or_open_nested_ML( $curr_shrep, \@c_mls, \@p_bp, \@up, \@stmbrks, \@open_blks, $idx, $curr_gi );
+
+              # OPEN_ML
+            } else {
+
+              # case 'x((x)(' the opening of the ML was a stem-break
+              if ( @stmbrks
+                && $stmbrks[-1]->[1] == $p_bp[0]
+                && $stmbrks[-1]->[2] == $p_bp[1] ) {
+                my @c_stmbrk = @{ pop(@stmbrks) };
+                my @a        = ( $c_stmbrk[0], -1 );
+                my @ml       = ();
+                push( @ml, \@a );    # 1st base-pair is in-waiting
+                my @a2 = ( $c_stmbrk[1], $c_stmbrk[2] );
+                push( @ml, \@a2 );    #2nd BP
+                my @a3 = ($idx);
+                push( @ml,    \@a3 );    #first of current BP
+                push( @c_mls, \@ml );
+                print STDERR "TEST - opened ML with SB:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
+
+                # case 'x(.(x)('
+              } elsif (@up) {
+                my @tmp_up = @{ pop(@up) };
+
+                # compare adjacent BP to UP to see if they fit
+                if ( $tmp_up[1] + 1 == $p_bp[0]
+                  && $open_blks[-1]->[-1] == $tmp_up[0] - 1 ) {
+                  my @a = ( $open_blks[-1]->[-1], -1 );
+                  my @ml = ();
+                  push( @ml, \@a );    # 1st BP awaits closing
+                  my @newbp = @p_bp;
+                  push( @ml, \@newbp );    # 2nd BP
+                  my @a2 = ($idx);
+                  push( @ml,    \@a2 );    # current open BP, awaits closing
+                  push( @c_mls, \@ml );
+                  print STDERR "TEST - opened ML with prev UP1:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) )
+
+                    if ($i_debug);
+                } else {
+                  die "ERROR: the base-pairs to not match the " . "previous unpaired region?!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+                }
+              } else {
+                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 ) );
+              }
+            }
+
+            # case: 'x(x)('
+          } else {
+
+            # case: '.(x)(' => CLOSE_EL
+            if (@c_el) {
+              if ( $c_el[-1] == $p_bp[0] ) {
+                push( @c_el, $p_bp[1] );
+                my @newel = @c_el;
+                push( @els, \@newel );
+                @c_el = ();
+                print STDERR "TEST - closed EL1:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newel ) ) if ($i_debug);
+              } else {
+                die "ERROR: the previous BP must match the current EL\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+              }
+            }
+
+            # in both cases OPEN_EL
+            push( @c_el, @p_bp );
+            push( @c_el, $idx );
+            print STDERR "TEST - opened EL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_el ) ) if ($i_debug);
+          }
+
+          # case '.('
+        } elsif (@c_up) {
+
+          # case (x(x).( OR 'x(.('
+          if (@open_blks) {
+
+            # case 'x(x).(' ML
+            if ( @p_bp && ( $p_bp[1] == $c_up[0] - 1 ) ) {
+
+              # case '((x)x(x).(' EXTEND_ML, if one exists
+              if (@c_mls) {
+
+                extend_or_open_nested_ML( $curr_shrep, \@c_mls, \@p_bp, \@up, \@stmbrks, \@open_blks, $idx, $curr_gi );
+
+                # case 'x((x).(' or 'x(.(x).(' OPEN_ML
+              } else {
+
+                # case 'x((x).(' the opening of the ML was a stem-break
+                if ( @stmbrks
+                  && $stmbrks[-1]->[1] == $p_bp[0]
+                  && $stmbrks[-1]->[2] == $p_bp[1] ) {
+                  my @c_stmbrk = @{ pop(@stmbrks) };
+                  my @a        = ( $c_stmbrk[0], -1 );
+                  my @ml       = ();
+                  push( @ml, \@a );    # 1st base-pair is in-waiting
+                  my @a2 = ( $c_stmbrk[1], $c_stmbrk[2] );
+                  push( @ml, \@a2 );    #2nd BP
+                  my @a3 = ($idx);
+                  push( @ml,    \@a3 );    #first of current BP
+                  push( @c_mls, \@ml );
+                  print STDERR "TEST - open ML with SB:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
+
+                  # case 'x(.(x).(' opening the ML with initial UP region
+                } elsif (@up) {
+
+                  # get previous unpaired region
+                  my $tmp  = pop(@up);
+                  my @p_up = @{$tmp};
+
+          # case: '(.(x).(' compare adjacent BP to UP to see if they create a ML
+                  if ( $p_up[1] + 1 == $p_bp[0]
+                    && $open_blks[-1]->[-1] == $p_up[0] - 1 ) {
+                    my @a = ( $open_blks[-1]->[-1], -1 );
+                    my @ml = ();
+                    push( @ml, \@a );    # 1st BP awaits closing
+                    my @newbp = @p_bp;
+                    push( @ml, \@newbp );    # 2nd BP
+                    my @a2 = ($idx);
+                    push( @ml,    \@a2 );    # current open BP, awaits closing
+                    push( @c_mls, \@ml );
+                    print STDERR "TEST - opened ML with prev UP2:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) )
+
+                      if ($i_debug);
+                  }
+
+                  # has to be an ML, but can't find first part
+                } else {
+                  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 ) );
+                }
+              }
+
+              # case '(.(
+            } else {
+              my @newup = @c_up;
+              push( @up, \@newup );
+              print STDERR "TEST - closed UP:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newup ) ) if ($i_debug);
+            }
+
+            # case: 'x(x).(' OR 'EOS.('
+          } else {
+
+            # case: '.(x)(' => CLOSE_EL
+            if (@c_el) {
+              if ( $c_el[-1] == $c_bp[0] ) {
+                push( @c_el, $p_bp[1] );
+                my @newel = @c_el;
+                push( @els, \@newel );
+                @c_el = ();
+                print STDERR "TEST - closed EL2:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newel ) ) if ($i_debug);
+              } else {
+                die "ERROR: the previous BP must match the current EL\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+              }
+            }
+
+            # in both cases OPEN_EL
+            if (@p_bp) {
+              push( @c_el, @p_bp );
+              push( @c_el, $idx );
+            } else {
+
+              # ignore this if -cue is given
+              # TODO check the gi index
+              if ($cue) {
+                @c_el = ();
+              } else {
+                push( @c_el, ( $curr_gi, $curr_gi ) );
+                push( @c_el, $idx );
+                print STDERR "TEST - opened EL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_el ) ) if ($i_debug);
+              }
+            }
+          }
+          @c_up = ();
+
+        } else {
+
+          # no more cases except opening bracket at beginning of sequence
+          # in this case do nothing, as the index has already been added
+          die "ERROR: '((', ')(', '.(' have all been covered\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) )
+            unless ( $idx == $curr_gi );
+
+        }
+      }
+
+      #===================================================
+      # current char is the unpaired base
+    } elsif ( $struct[$idx] eq "." ) {
+
+      #===================================================
+
+      # update current window index to current graph index
+      $idx += $curr_gi;
+
+      # case ').' => CLOSE_STEM
+      if (@c_stem) {
+        print STDERR "TEST - closed stem:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_stem ) ) if ($i_debug);
+        close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps );
+      }
+
+      # case '..' extend UP, @c_up
+      if (@c_up) {
+        $c_up[1] = $idx;
+
+        # case 'x.' open new UP
+      } else {
+
+        # OPEN_UP, @c_up
+        @c_up = ( $idx, $idx );
+
+        # case '(.' just come to the end of an open block, push it onto stack
+        if (@c_open) {
+          my @newopen = @c_open;
+          push( @open_blks, \@newopen );
+          @c_open = ();
+
+          # case ').' or '.'
+          # (have already closed stem and created stem-break and emptied p_open)
+        } else {
+
+          # case 'x(x).' or '.' => CLOSE_EL
+          if ( !@open_blks ) {
+
+            # case 'x().().'
+            if (@c_el) {
+
+              # double check
+              if ( $c_el[-1] == $p_bp[0] ) {
+                push( @c_el, $p_bp[1] );
+                my @newel = @c_el;
+                push( @els, \@newel );
+                @c_el = ();
+                print STDERR "TEST - closed EL3:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newel ) ) if ($i_debug);
+              } else {
+                die "This case should not occur!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+              }
+            }
+          }
+        }
+      }
+
+      #===================================================
+      # current char is the closing bracket of a base-pair
+    } elsif ( $struct[$idx] eq ")" ) {
+
+      #===================================================
+
+      # update current window index to current graph index
+      $idx += $curr_gi;
+
+      # save previous base-pair
+      if (@c_bp) {
+        @p_bp = @c_bp;
+      }
+
+      # case: '((x))', extend stem
+      if (@p_open) {
+
+        # get current base-pair
+        @c_bp = ( pop(@p_open), $idx );
+
+        # add it to the base-pair hash
+        $bps{"$c_bp[0]:$c_bp[1]"} = "";
+
+        # EXTEND_STEM
+        if (@c_stem) {
+          $c_stem[0] = $c_bp[0];
+          $c_stem[1] = $c_bp[1];
+        } else {
+          die "ERROR: in case '))' there has to be an open stem object!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+        }
+
+        # case: '(x(x))' or '(x.)', '()' not allowed
+      } else {
+
+        # case: '(x(x))'  CLOSE_STEM close previous stem
+        if (@c_stem) {
+          print STDERR "TEST - closed stem:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@c_stem ) ) if ($i_debug);
+          close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps );
+
+          # case '(x(x))' OPEN_STEM
+          if (@open_blks) {
+            my $tmp = pop(@open_blks);
+            @p_open = @{$tmp};
+
+            # get current base-pair
+            @c_bp = ( pop(@p_open), $idx );
+
+            # add it to the base-pair hash
+            # as the stem is not finished yet, cannot add stem information
+            $bps{"$c_bp[0]:$c_bp[1]"} = "";
+
+            # open stem
+            $c_stem[0] = $c_bp[0];
+            $c_stem[1] = $c_bp[1];
+            $c_stem[2] = $c_bp[0];
+            $c_stem[3] = $c_bp[1];
+          } else {
+            die "ERROR: there have to be open blocks to match " . "current closing bracket!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+          }
+
+          # case: '(.(x))' => CLOSE_BL
+          if ( @up && $up[-1]->[0] - 1 == $c_bp[0] && $up[-1]->[1] + 1 == $p_bp[0] ) {
+            my @newBL = ( @c_bp, @p_bp );
+            push( @bls, \@newBL );
+            pop(@up);
+            print STDERR "TEST - closed BL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newBL ) ) if ($i_debug);
+
+            # case: '(x(x).(x))' => CLOSE_ML
+          } elsif ( @c_mls && $c_mls[-1]->[-1]->[-1] == $p_bp[0] ) {
+            my @newml = @{ pop(@c_mls) };
+            push( @{ $newml[-1] }, $p_bp[1] );    # close last BP
+            $newml[0]->[1] = $c_bp[1];            # close 1st BP
+            push( @mls, \@newml );
+            print STDERR "TEST - closed ML1:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newml, 1 ) ) if ($i_debug);
+          } else {
+            die "ERROR: what is this case?\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+          }
+
+          # case: (x.)
+        } else {
+
+          # OPEN_STEM open new stem
+          if (@open_blks) {
+            @p_open = @{ pop(@open_blks) };
+
+            # get current base-pair
+            @c_bp = ( pop(@p_open), $idx );
+
+            # add it to the base-pair hash
+            # as the stem is not finished yet, cannot add stem information
+            $bps{"$c_bp[0]:$c_bp[1]"} = "";
+
+            # open stem
+            $c_stem[0] = $c_bp[0];
+            $c_stem[1] = $c_bp[1];
+            $c_stem[2] = $c_bp[0];
+            $c_stem[3] = $c_bp[1];
+          } else {
+            die "ERROR: there have to be open blocks to match " . "current closing bracket!$curr_shrep->[0]\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+          }
+
+          # CLOSE_C_UP
+          if (@c_up) {
+
+            # case: (.) => close_HL
+            if ( $c_up[0] - 1 == $c_bp[0] ) {
+              my @newhl = @c_bp;
+              push( @hls, \@newhl );
+              @c_up = ();
+              print STDERR "TEST - closed HL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newhl ) ) if ($i_debug);
+
+              # case: (x(x).)
+            } elsif ( @p_bp && $c_up[0] - 1 == $p_bp[1] ) {
+
+              # case: ((x).) => CLOSE_BL
+              if ( @stmbrks
+                && $stmbrks[-1]->[0] == $c_bp[0]
+                && $stmbrks[-1]->[1] == $p_bp[0]
+                && $stmbrks[-1]->[2] == $p_bp[1] ) {
+                my @newbl = ( @c_bp, @p_bp );
+                push( @bls, \@newbl );
+                pop(@stmbrks);
+                @c_up = ();
+                print STDERR "TEST - closed BL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newbl ) ) if ($i_debug);
+
+                # case: (x(x).), either CLOSE_ML or CLOSE_IL
+              } else {
+
+                # case '(.(x).)' => CLOSE_IL
+                if ( @up
+                  && ( $up[-1]->[0] - 1 == $c_bp[0] )
+                  && ( $up[-1]->[1] + 1 == $p_bp[0] ) ) {
+                  my @p_up = @{ pop(@up) };
+                  my @newil = ( @c_bp, @p_bp );
+                  push( @ils, \@newil );
+                  @c_up = ();
+                  print STDERR "TEST - closed IL:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newil ) ) if ($i_debug);
+
+                  # case (x(x)x(x).) => CLOSE_ML remember, array of tuples
+                } elsif ( @c_mls && $c_mls[-1]->[-1]->[0] == $p_bp[0] ) {
+                  my @newml = @{ pop(@c_mls) };
+                  push( @{ $newml[-1] }, $p_bp[1] );    # close final BP
+                  $newml[0]->[1] = $c_bp[1];            # close 1st BP
+                  push( @mls, \@newml );
+                  @c_up = ();
+                  print STDERR "TEST - closed ML2:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@newml, 1 ) ) if ($i_debug);
+                } else {
+                  die "There should be no such case!\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+                }
+              }
+
+              # there should be no more unknown unpaired regions left!
+            } else {
+              die "There shouldn't be any unknown unpaired " . "regions left\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+            }
+
+          } else {
+            die "ERROR: '()' is not allowed in $curr_shrep->[0]\n" . ( test( $curr_shrep->[0], $idx, $curr_gi ) );
+          }
+        }
+      }
+    }
+  }
+
+  # EOS
+
+  # CLOSE_STEM
+  print STDERR "TEST - closed stem EOS:\n" . ( test( $curr_shrep->[0], ( $win_size + $curr_gi ), $curr_gi, \@c_stem ) )
+    if ( @c_stem && $i_debug );
+  close_stem( \@c_stem, \%stems, \@p_bp, \@c_bp, \@p_open, \@stmbrks, \@open_blks, \%bps ) if (@c_stem);
+  die "ERROR: there should not be an open ML!\n" . $curr_shrep . "\n"
+    if (@c_mls);
+
+  # case ().() CLOSE_EL
+  if (@c_el) {
+    if ( $c_el[-1] == $p_bp[0] ) {
+      push( @c_el, $p_bp[1] );
+      my @newel = @c_el;
+      push( @els, \@newel );
+      @c_el = ();
+
+      print STDERR "TEST - closed EL4 EOS:\n" . ( test( $curr_shrep->[0], ( $win_size + $curr_gi ), $curr_gi, \@newel ) )
+        if ($i_debug);
+    } else {
+      die "ERROR: the base-pairs in the EL don't match" . $curr_shrep . "\n";
+    }
+  }
+
+  # case (). CLOSE_EL
+  if (@c_up) {
+    unless ($cue) {
+      if (@p_bp) {
+        push( @c_el, @p_bp );
+      } else {
+        push( @c_el, ( $curr_gi, $curr_gi ) );
+      }
+
+      push( @c_el, ( ( $curr_gi + $win_size - 1 ), ( $curr_gi + $win_size - 1 ) ) );
+      my @newel = @c_el;
+      push( @els, \@newel );
+      @c_el = ();
+      print STDERR "TEST - closed EL5 EOS:\n" . ( test( $curr_shrep->[0], ( $win_size + $curr_gi ), $curr_gi, \@newel ) )
+        if ($i_debug);
+      ## check case where cue is given, but window is completely unstructured
+    } else {
+      unless (@p_bp) {
+        my @newel = ( $curr_gi, $curr_gi, ( $curr_gi + $win_size - 1 ), ( $curr_gi + $win_size - 1 ) );
+        push( @els, \@newel );
+      }
+    }
+  }
+
+  #  # update current graph index (since all shrep indices have been read)
+  $curr_gi += $win_size;
+
+  # ===================================================
+  # build graph lines
+  my @graph_lines = ();
+
+  # ===================================================
+
+  # graph line objects
+  my @edg        = ();
+  my @stackgraph = ();
+  my @abstrgraph = ();
+
+ #-------------------> stems,stacks
+ # build basic graph edges, stacks, and vertices for stems in the abstract graph
+  my @tmp = keys %stems;
+  foreach my $sk (@tmp) {
+    my @stem = @{ $stems{$sk} }; # stem (i,j,k,l)=> (i,j)= (0,1) and (k,l)=(2,3)
+    die "ERROR: There are not enough elements in this stem: ", ( join( ",", @stem ) ), "\n" unless ( scalar(@stem) == 4 );
+
+    # infer base-pairs from stem
+    @p_bp = ();
+    @c_bp = ();
+    my @stem_relation_idxs  = ();
+    my @stack_relation_idxs = ();
+    my $stacksize           = $stem[2] - $stem[0] + 1;
+    for ( my $x = 0 ; $x < $stacksize ; $x++ ) {
+
+      #    for (my $i = $stem[0]; $i <= $stem[2] ; $i++){
+      #      for(my $j = $stem[1]; $j >= $stem[3]; $j--){
+      my $i = $stem[0] + $x;
+      my $j = $stem[1] - $x;
+      @p_bp = @c_bp;
+      @c_bp = ( $i, $j );
+      if ( defined $bps{"$i:$j"} ) {
+
+        # add indices to stem relation object
+        push( @stem_relation_idxs, $i );
+        push( @stem_relation_idxs, $j );
+
+        # add edge to graph edges
+        push( @edg, "e $i $j s" );
+
+        # add stack if option given
+        if ($stacks) {
+          if (@p_bp) {
+            push( @stackgraph,          "v $curr_gi P" );
+            push( @stackgraph,          "e $curr_gi $p_bp[0] p" );
+            push( @stackgraph,          "e $curr_gi $p_bp[1] p" );
+            push( @stackgraph,          "e $c_bp[0] $curr_gi p" );
+            push( @stackgraph,          "e $c_bp[1] $curr_gi p" );
+            push( @stack_relation_idxs, $curr_gi );
+
+            # update graph index
+            ++$curr_gi;
+          }
+        }
+      } else {
+        die "ERROR: There is no such base-pair in the base-pair hash:" . " ($i, $j)\n";
+      }
+    }
+
+    #  	  }
+    #    }
+    # add stem to abstract graph
+    # (graph indices become a bit jumbled, but saves time)
+    push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#S" );
+
+    # stem object now saves the curr_gi index for that stem
+    $stems{$sk} = $curr_gi;
+    ++$curr_gi;
+
+    # add relations
+    push( @abstrgraph, "v $curr_gi ^R" );
+    push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
+    foreach my $r (@stem_relation_idxs) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    if ($stacks) {
+      foreach my $r (@stack_relation_idxs) {
+        push( @abstrgraph, "e $curr_gi $r \@r" );
+      }
+    }
+    ++$curr_gi;
+  }
+
+  #-------------------> add HLs
+  foreach my $hl (@hls) {
+
+    # check size
+    die "ERROR: the HL object is the incorrect size\n" unless ( @{$hl} == 2 );
+
+    # add vertex
+    push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#HL" );
+
+    # add edge from stem to hl
+    if ( defined( $stems{ $bps{"$hl->[0]:$hl->[1]"} } ) ) {
+      push( @abstrgraph, "e " . $stems{ $bps{"$hl->[0]:$hl->[1]"} } . " $curr_gi n" );
+    } else {
+      die "ERROR: stem not defined in base-pair hash for $hl->[0]:$hl->[1]\n";
+    }
+    ++$curr_gi;
+
+    # add relations
+    push( @abstrgraph, "v $curr_gi ^R" );
+    push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
+
+    # add all nodes between i to j, inclusively
+    for ( my $r = $hl->[0] ; $r <= $hl->[1] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    ++$curr_gi;
+  }
+
+  #------------------> add BLs
+  foreach my $bl (@bls) {
+
+    # check size
+    die "ERROR: the BL object is the incorrect size\n" unless ( scalar( @{$bl} ) == 4 );
+
+    # add vertex
+    push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#BL" );
+
+    # add edges from bl to both adjacent stems
+    if ( defined( $stems{ $bps{"$bl->[0]:$bl->[1]"} } )
+      && defined( $stems{ $bps{"$bl->[2]:$bl->[3]"} } ) ) {
+
+      # add outer stem: S->BL
+      push( @abstrgraph, "e " . $stems{ $bps{"$bl->[0]:$bl->[1]"} } . " $curr_gi n" );
+
+      # add inner stem: BL->S
+      push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$bl->[2]:$bl->[3]"} } . " n" );
+    } else {
+      die "ERROR: stems not defined for BL\n";
+    }
+    ++$curr_gi;
+
+    # add relations
+    push( @abstrgraph, "v $curr_gi ^R" );
+    push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
+
+    # add all nodes between i to k and l to j, inclusively
+    for ( my $r = $bl->[0] ; $r <= $bl->[2] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    for ( my $r = $bl->[3] ; $r <= $bl->[1] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    ++$curr_gi;
+  }
+
+  #------------------> add ILs
+  foreach my $il (@ils) {
+
+    # check size
+    die "ERROR: the IL object is the incorrect size\n" unless ( scalar( @{$il} ) == 4 );
+
+    # add vertex
+    push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#IL" );
+
+    # add edges from il to both adjacent stems
+    if ( defined( $stems{ $bps{"$il->[0]:$il->[1]"} } )
+      && defined( $stems{ $bps{"$il->[2]:$il->[3]"} } ) ) {
+
+      # add outer stem: S->IL
+      push( @abstrgraph, "e " . $stems{ $bps{"$il->[0]:$il->[1]"} } . " $curr_gi n" );
+
+      # add inner stem: IL->S
+      push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$il->[2]:$il->[3]"} } . " n" );
+    } else {
+      die "ERROR: stems not defined for IL\n";
+    }
+    ++$curr_gi;
+
+    # add relations
+    push( @abstrgraph, "v $curr_gi ^R" );
+    push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
+
+    # add all nodes between i to k and l to j, inclusively
+    for ( my $r = $il->[0] ; $r <= $il->[2] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    for ( my $r = $il->[3] ; $r <= $il->[1] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    ++$curr_gi;
+  }
+
+  #------------------> add ELs
+  foreach my $el (@els) {
+
+    #check size
+    die "ERROR: the EL object is the incorrect size\n" unless ( scalar( @{$el} ) == 4 );
+
+    # add vertex
+    push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#EL" );
+
+    # add edges from el to both adjacent stems, if available
+    unless ( $el->[0] == $el->[1] ) {
+      if ( defined( $stems{ $bps{"$el->[0]:$el->[1]"} } ) ) {
+        push( @abstrgraph, "e " . $stems{ $bps{"$el->[0]:$el->[1]"} } . " $curr_gi n" );
+      } else {
+        die "ERROR: stems not defined for left BP of EL\n";
+      }
+    }
+    unless ( $el->[2] == $el->[3] ) {
+      if ( defined( $stems{ $bps{"$el->[2]:$el->[3]"} } ) ) {
+        push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$el->[2]:$el->[3]"} } . " n" );
+      } else {
+        die "ERROR: stems not defined for right BP of EL\n";
+      }
+    }
+    ++$curr_gi;
+
+    # add relations
+    push( @abstrgraph, "v $curr_gi ^R" );
+    push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
+
+    # @r relation edges
+    for ( my $r = $el->[1] ; $r <= $el->[2] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+
+    # if BP, not EOS
+    unless ( $el->[0] == $el->[1] ) {
+      push( @abstrgraph, "e $curr_gi $el->[0] \@r" );
+    }
+
+    # if BP, not EOS
+    unless ( $el->[2] == $el->[3] ) {
+      push( @abstrgraph, "e $curr_gi $el->[3] \@r" );
+    }
+    ++$curr_gi;
+  }
+
+  #------------------> add MLs
+  foreach my $ml (@mls) {
+
+    #check size
+    die "ERROR: the ML object is not big enough!\n" unless ( scalar( @{$ml} ) >= 3 );
+
+    # add vertex
+    push( @abstrgraph, "v $curr_gi ${ABSTRUCT}#ML" );
+
+    # add edge from outer stem to first BP
+    my $closing_bp = shift( @{$ml} );
+    if ( defined( $stems{ $bps{"$closing_bp->[0]:$closing_bp->[1]"} } ) ) {
+      push( @abstrgraph, "e " . $stems{ $bps{"$closing_bp->[0]:$closing_bp->[1]"} } . " $curr_gi n" );
+    } else {
+      die "ERROR: ML is not defined for BP $closing_bp->[0]:$closing_bp->[1]\n";
+    }
+
+    # add remaining stems
+    foreach my $bp ( @{$ml} ) {
+      if ( defined( $stems{ $bps{"$bp->[0]:$bp->[1]"} } ) ) {
+        push( @abstrgraph, "e $curr_gi " . $stems{ $bps{"$bp->[0]:$bp->[1]"} } . " n" );
+      } else {
+        die "ERROR: ML is not defined for BP $bp->[0]:$bp->[1]\n";
+      }
+    }
+    ++$curr_gi;
+
+    # add relations
+    push( @abstrgraph, "v $curr_gi ^R" );
+    push( @abstrgraph, "e " . ( $curr_gi - 1 ) . " $curr_gi ^r" );
+
+  # add @r relation edges between i and k, where (i,j) is the first BP and (k,l)
+  # the second BP
+    my $pbp_aref = $closing_bp;
+    my $cbp_aref = shift( @{$ml} );
+    for ( my $r = $pbp_aref->[0] ; $r <= $cbp_aref->[0] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+
+    # rest of @r relation edges
+    foreach my $bp ( @{$ml} ) {
+      $pbp_aref = $cbp_aref;
+      $cbp_aref = $bp;
+      for ( my $r = $pbp_aref->[1] ; $r <= $cbp_aref->[0] ; $r++ ) {
+        push( @abstrgraph, "e $curr_gi $r \@r" );
+      }
+    }
+
+    # add final stretch between last BP (m,n) and first closing base-pair (i,j)
+    for ( my $r = $cbp_aref->[1] ; $r <= $closing_bp->[1] ; $r++ ) {
+      push( @abstrgraph, "e $curr_gi $r \@r" );
+    }
+    ++$curr_gi;
+  }
+
+  @graph_lines = ( @edg, @stackgraph, @abstrgraph );
+  return ( \@graph_lines, $curr_gi );
+}
+
+# just to test the current shrep dot-bracket structure parsing
+sub test {
+  my ( $shrep, $idx, $gi, $mark_aref, $isML ) = @_;
+
+  my $result = "";
+
+  $result .= "$shrep\n";
+  my @a = split( "", $shrep );
+  for ( my $i = 0 ; $i <= @a ; $i++ ) {
+    my $mark = 0;
+    foreach my $j ( @{$mark_aref} ) {
+      if ($isML) {
+        $mark = 1 if ( $i == $j->[0] - $gi );
+        $mark = 1 if ( scalar( @{$j} == 2 && $i == $j->[1] - $gi ) );
+      } else {
+        $mark = 1 if ( $i == $j - $gi );
+      }
+
+    }
+    if ( $i == $idx - $gi ) {
+      $result .= "^";
+    } elsif ($mark) {
+      $result .= "'";
+    } else {
+      $result .= " ";
+    }
+  }
+  $result .= "\n";
+
+  return $result;
+}
+
+###########################################################################
+# Extends or closes a multi-loop. This means a ML already exists.
+# We can have the case ((x)( or ((x).(
+#
+# Input:
+# curr_shrep	The current shrep object, with shrep string at pos 0
+# c_ml_aref		The current open MLs array reference
+# p_bp_aref		The previous base-pair array reference
+# stmbrks_aref	The stem-breaks array reference
+# open_blks_aref	The open bracket blocks array reference
+# idx			The current index
+# #curr_gi		The current graph index
+# Output:
+# None, it should just modify the ML accordingly
+############################################################################
+sub extend_or_open_nested_ML {
+  my ( $curr_shrep, $c_mls_aref, $p_bp_aref, $up_aref, $stmbrks_aref, $open_blks_aref, $idx, $curr_gi ) = @_;
+
+  # case '(x(x)x(x)x(' extend current ML
+  if ( $c_mls_aref->[-1]->[-1]->[-1] == $p_bp_aref->[0] ) {
+    push( @{ $c_mls_aref->[-1]->[-1] }, $p_bp_aref->[1] );    #close previous BP
+    my @a = ($idx);
+    push( @{ $c_mls_aref->[-1] }, \@a );    # add new open base-pair
+    print STDERR "TEST - extended ML:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, $c_mls_aref->[-1], 1 ) ) if ($i_debug);
+
+    # case '(x(x)x((x)(' nested ML, => OPEN_ML with stem-break
+  } elsif ( @$stmbrks_aref
+    && $stmbrks_aref->[-1]->[1] == $p_bp_aref->[0]
+    && $stmbrks_aref->[-1]->[2] == $p_bp_aref->[1] ) {
+
+    my @c_stmbrk = @{ pop( @{$stmbrks_aref} ) };
+    my @a        = ( $c_stmbrk[0], -1 );
+    my @ml       = ();
+    push( @ml, \@a );                       # 1st base-pair is in-waiting
+    my @a2 = ( $c_stmbrk[1], $c_stmbrk[2] );
+    push( @ml, \@a2 );                      #2nd BP
+    my @a3 = ($idx);
+    push( @ml,            \@a3 );           #first of current BP
+    push( @{$c_mls_aref}, \@ml );
+    print STDERR "TEST - opened nested ML with SB:\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
+
+    # case '(x(x)x(.(x)(' nested ML, => OPEN_ML with UP
+  } elsif ( @$up_aref
+    && $up_aref->[-1]->[0] - 1 == $open_blks_aref->[-1]->[-1]
+    && $up_aref->[-1]->[1] + 1 == $p_bp_aref->[0] ) {
+    my @tmp_up = @{ pop( @{$up_aref} ) };
+    my @a      = ( $open_blks_aref->[-1]->[-1], -1 );
+    my @ml     = ();
+    push( @ml, \@a );                       # 1st BP awaits closing
+    my @newbp = @{$p_bp_aref};
+    push( @ml, \@newbp );                   #2nd BP
+    my @a2 = ($idx);
+    push( @ml,            \@a2 );           # current open BP, awaits closing
+    push( @{$c_mls_aref}, \@ml );
+    print STDERR "TEST - opened nested ML with prev UP\n" . ( test( $curr_shrep->[0], $idx, $curr_gi, \@ml, 1 ) ) if ($i_debug);
+  } else {
+    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);
+  }
+}
+
+###########################################################################
+# Closes a stem-object for the abstract graph. If necessary, it also
+# creates a stem-break because the current open block is still open,
+# which can only mean that there is a stem-break between two consecutive
+# opening brackets.
+#
+# Input:
+# c_stem_aref : current stem array reference
+# stems_href : all stems hash reference
+# p_bp_aref : previous base-pair array reference
+# c_cp_aref : current base-pair array reference
+# p_open_aref : previously opened open-bracket block
+# p_stmbrk_aref : the previous stem-break
+# open_blks_aref : all old open blocks array reference
+# bps_href : all base-pairs hash reference
+#
+# Output:
+# None, it should just modify the input variables accordingly
+############################################################################
+sub close_stem {
+  my ( $c_stem_aref, $stems_href, $p_bp_aref, $c_bp_aref, $p_open_aref, $stmbrks_aref, $open_blks_aref, $bps_href ) = @_;
+
+  my @newstem = @$c_stem_aref;
+  my $stemkey = "$newstem[0]:$newstem[1],$newstem[2]:$newstem[3]";
+  $stems_href->{$stemkey} = \@newstem;
+
+  # reset and clean running variables
+  @$c_stem_aref = ();
+  @$p_bp_aref   = @$c_bp_aref;
+  @$c_bp_aref   = ();
+
+  # add stem to base-pair hash
+  $bps_href->{"$newstem[0]:$newstem[1]"} = $stemkey;
+  $bps_href->{"$newstem[2]:$newstem[3]"} = $stemkey;
+
+  # there are still bases open in the previous open block, '[[]x',
+  # where [] is a closed stem and [ is an open block,
+  # must create STEM-BREAK
+  # this cannot be the case ((x)), because there p_open is empty
+  if (@$p_open_aref) {
+
+    # check: the next open base is adjacent to the previous BP
+    if ( $p_open_aref->[-1] == $p_bp_aref->[0] - 1 ) {
+      my @stmbrk = ( $p_open_aref->[-1], @$p_bp_aref );
+      push( @$stmbrks_aref, \@stmbrk );
+
+      # push the previous open block back onto stack
+      my @newopen = @$p_open_aref;
+      push( @$open_blks_aref, \@newopen );
+      @$p_open_aref = ();
+    } else {
+
+      # error, this should not be possible
+      die("ERROR in close_stem(), [.[]x should not occur here\n");
+    }
+  }
+}
+
+###########################################################################
+# Here the information about the structure is added to the graph. In a
+# first step each base-pair is added by adding edges between the vertices
+# of the respective nucleotides. Then, if two base-pairs are stacked,
+# we add an extra vertex and connect the four nucleotide vertices that
+# are involved with extra edges (symbols P used for vertex and p for edges).
+# While adding new vertices we keep track of the current graph index.
+#
+# Input:
+# curr_shrep : the current shrep as dot-bracket structure
+# curr_gi : the current graph index
+# win_size : the window size
+# stacks : the input parameter -stacks (whether to add stacking information)
+#
+# Output:
+# The structure graph information line-by-line as an array reference and
+# the current graph index.
+############################################################################
+sub getStructGraph {
+  my ( $curr_shrep, $curr_gi, $win_size, $stacks ) = @_;
+
+  my @struct = split( "", $curr_shrep->[0] );
+
+  my @edg    = ();
+  my @starts = ();
+  my @pairs  = ();
+
+  foreach my $idx ( 0 .. @struct - 1 ) {
+
+    push( @starts, $idx ) if ( $struct[$idx] eq "(" );
+
+    if ( $struct[$idx] eq ")" ) {
+      my $start = $curr_gi + pop(@starts);
+      my $end   = $curr_gi + $idx;
+      push( @edg, "e " . $start . " " . $end . " s " );
+      my @pair = ( $start, $end );
+      push( @pairs, \@pair );
+    }
+  }
+
+  # update current graph index
+  $curr_gi = $curr_gi + $win_size;
+
+  # initiate structure lines array
+  my @stacking_info = ();
+
+  # add stacking information to graph unless input option tells us not to
+  if ($stacks) {
+
+    my $stacked_pairs = 0;
+
+    # get stacked base-pairs (they are ordered according to
+    # position of closing nucleotide)
+    for ( my $i = 1 ; $i < @pairs ; $i++ ) {
+
+      # add stacked base-pairs (vertices+edges)
+      my $curr_bp_aref = $pairs[$i];
+      my $prev_bp_aref = $pairs[ $i - 1 ];
+
+      # if the current base-pair is stacked on the previous base-pair,
+      # when curr_start = prev_start - 1 AND curr_end = prev_end + 1
+      if ( $curr_bp_aref->[0] == $prev_bp_aref->[0] - 1
+        && $curr_bp_aref->[1] == $prev_bp_aref->[1] + 1 ) {
+
+        # add stacking vertex P
+        push( @stacking_info, "v $curr_gi P" );
+
+        # add four edges from involved nucleotids
+        push( @stacking_info, "e $curr_gi $prev_bp_aref->[0] p" );
+        push( @stacking_info, "e $curr_gi $prev_bp_aref->[1] p" );
+        push( @stacking_info, "e $curr_bp_aref->[0] $curr_gi p" );
+        push( @stacking_info, "e $curr_bp_aref->[1] $curr_gi p" );
+        ++$curr_gi;    # add one to the index, ready for next vertex
+
+      }
+    }
+  }
+
+  my @str_graphlines = ( @edg, @stacking_info );
+  return ( \@str_graphlines, $curr_gi );
+}
+
+###########################################################################
+# This method gives us the overall graph header for this sequence. This
+# graph includes all window calculations and all shreps for all windows.
+# Input:
+# seq_id : the sequence ID
+# seq_head : the header information for the sequence
+# wins_aref : the windows sizes for folding
+# h_shift : input parameter -shift
+# h_e : input parameter -e
+# h_c : input parameter -c
+# h_t : input parameter -t
+# h_u : input parameter -u
+# h_r : input parameter -r
+# h_M : input parameter -M
+# h_crop_unpaired_ends : input parameter -cue
+#
+# Output:
+# The header line as a string
+############################################################################
+sub getGraphHeader {
+  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 ) = @_;
+
+  my $ret;
+
+  $ret = "t # SEQID $seq_id ";
+  $ret .= "$seq_head " if ( defined $seq_head );
+  $ret .= "SEQLEN $h_seqlen ";
+
+  return $ret . "\n" if ($h_noStr);
+
+  $ret .= "MAXWINSHREPS $h_M ";
+  $ret .= "CUE 1 " if ($h_crop_unpaired_ends);
+  $ret .= "RNASHAPES -w " . ( join( ",", @{$wins_aref} ) ) . " ";
+  $ret .= "-SHIFT%w $h_shift " if ( not $globalFolding and $h_shift );
+  $ret .= "-SHIFT 1 " if ( not $globalFolding and not $h_shift );
+  $ret .= "-e $h_e "                         if ( defined $h_e );
+  $ret .= "-c $h_c "                         if ( defined $h_c );
+  $ret .= "-t $h_t "                         if ( defined $h_t );
+  $ret .= "-u 1 "                            if ($h_u);
+  $ret .= "-r 1 "                            if ($h_r);
+  $ret .= "-i $h_i SAMPLE_LEN $h_sample_len" if ( not $h_q and $h_i );
+  $ret .= "-q 1 "                            if ($h_q);
+  $ret .= "-T $h_Tp"                         if ( $h_q and defined $h_Tp );
+  $ret .= "\n";
+
+  return $ret;
+}
+
+###########################################################################
+# This method gives us the window header for the subgraph that includes
+# all shreps for the current window.
+# Input:
+# graphHead : the header line for the entire sequence graph
+# win_size : the size of the current window
+# win_start : the starting position of the current window
+# win_end : the end position of the current window
+# win_center : the centre position of the current window
+# win_global_fold : (boolean) whether the actual number of nucleotides
+# is smaller than the given window size (end of sequence?)
+#
+# Output:
+# The header line as a string
+############################################################################
+sub getWindowHeader {
+  my ( $graphHead, $win_size, $win_start, $win_end, $win_center, $win_global_fold, $win_sample, $used_t ) = @_;
+
+  chomp $graphHead;
+
+  $graphHead =~ s/t # //;
+
+  my $ret = "w # $graphHead ";
+  $ret .= "SHAPE_TYPE $used_t ";
+  $ret .= "GLOBALFOLD $win_global_fold ";
+  $ret .= "WSIZE $win_size ";
+  $ret .= "WSTART $win_start ";
+  $ret .= "WEND $win_end ";
+  $ret .= "WCENT $win_center ";
+  $ret .= "WSAMPLE $win_sample ";
+  $ret .= "\n";
+
+  return $ret;
+}
+
+###########################################################################
+# This method gives us the shape header that includes the single shrep
+# connected components.
+# Input:
+# winHead : the header for the current window
+# shrep : information about the current shrep
+#
+# Output:
+# The header line as a string
+############################################################################
+sub getShapeHeader {
+  my ( $winHead, $shrep ) = @_;
+
+  chomp $winHead;
+
+  $winHead =~ s/w # //;
+  $winHead =~ s/t # //;
+
+  my $ret = "s # $winHead ";
+
+  my @info = @{$shrep};
+  $ret .= join( " ", @info[ 1 .. $#info ] );
+  $ret .= "\n";
+
+  return $ret;
+}
+
+###########################################################################
+# This method gives us the sequence header that includes the unstructured
+# sequence connected components.
+# Input:
+# winHead : the header for the current window
+# shrep : information about the current shrep
+#
+# Output:
+# The header line as a string
+############################################################################
+sub getSeqHeader {
+  my ( $winHead, $shrep ) = @_;
+
+  chomp $winHead;
+
+  $winHead =~ s/w # //;
+  $winHead =~ s/t # //;
+
+  my $ret = "u # $winHead ";
+
+  my @info = @{$shrep};
+  $ret .= join( " ", @info[ 1 .. $#info ] );
+  $ret .= "\n";
+
+  return $ret;
+}
+
+
+##################################################################################
+# This method parses a fasta file and is useful if the header ID is not-unique!!
+# It returns the header lines in an array and the sequence lines in the same
+# order. It is then up to the user to extract the parts from the header that is
+# necessary for the script.
+# Furthermore, the method deals with multiple lines, and returns a single sequence
+# without line breaks.
+# Input:
+#		file		The name of the fasta file
+# Output:
+#	(1)	An array reference for each header line
+#	(2) An array reference for each sequence in same order as the headers
+##################################################################################
+sub read_fasta_with_nonunique_headers_meta {
+  my ($file) = @_;
+  my $FUNCTION = "read_fasta_file in Sequences.pm";
+
+  my $header    = "";
+  my $seqstring = "";
+  my @headers   = ();
+  my @sequences = ();
+  my @meta      = ();
+  my %seq_meta  = ();
+  open( IN_HANDLE, "<$file" ) || die "ERROR in $FUNCTION:\n" . "Couldn't open the following file in package Tool," . " sub read_fasta_file: $file\n";
+
+  while ( my $line = <IN_HANDLE> ) {
+    chomp($line);
+
+    # header (can contain one space after > symbol)
+    if ( $line =~ /^\>(.*)/ ) {
+      if ($header) {
+        $seqstring =~ s/\s*//g;    ## do not allow spaces in sequence
+        push( @headers,   $header );
+        push( @sequences, $seqstring );
+        push( @meta,      +{%seq_meta} );    ## anonymous hash reference
+                                             #print keys %seq_meta;
+        $seqstring = "";
+        undef(%seq_meta);
+      }
+      $header = $1;
+
+    } elsif ( $line =~ /(.+)\s+(#\S+)$/ && $header ) {
+
+      if ( exists $seq_meta{$2} ) {
+        $seq_meta{$2} .= $1;
+      } else {
+        $seq_meta{$2} = $1;
+      }
+
+    } elsif ($header) {
+      $seqstring .= $line
+    }
+  }
+
+  if ($header) {
+    $seqstring =~ s/\s*//g;    ## do not allow spaces in sequence
+    push( @headers,   $header );
+    push( @sequences, $seqstring );
+    push( @meta,      +{%seq_meta} );
+    $seqstring = "";
+    %seq_meta  = ();
+  }
+  return ( \@headers, \@sequences, \@meta );
+}
+
+
+1;
+
+#############################################################################
+# Programming description
+#############################################################################
+#	Substructure graphs for machine learning with Fabrizio
+#	-------------------------------------------------------
+#
+#	(1) Parameters (RNAshapes parameter):
+#		- Window sizes [] (-w)
+#		- window shift size (-W)
+#	- calculate structure probabilities (-r)
+#		- energy range kcal/mol (-e) OR energy relative percentage (%) to MFE (-c)
+#		- shape type 1 to 5 (-t)
+#		- ignore unstable substructures (-u)
+#		- max shreps
+#
+#
+#	(2) For each sequence, generate one graph/file that consists of all windows. The general format for one graph is as follows:
+#
+#	t # seq_id parameters
+#	v graph_index nt_type window_size window_centre abs_seq_pos
+#	...
+#	e m m+1 > 1 (backbone)
+#	...
+#	e base_i_graph_index base_j_graph_index s shrep_e shrep_p ...
+#
+#	For each window (subgraph) we create a subgraph (of the subgraph) for each substructure.
+#	We have a running index (gi) for each subgraph. All vertex and edge indices of the subgraph add
+#	the running graph index to the actual window position. For example
+#
+#
+#	Sequence: AAACC CUUUG GG
+#		  01234 56789 01
+#
+#	Window=10 substructure1 = (((...))). centre 5.5
+#
+#	v 0 A 10 5.5 0
+#	v 1 A 10 5.5 1
+#	v 2 A 10 5.5 2
+#	v 3 C 10 5.5 3
+#	v 4 C 10 5.5 4
+#	v 5 C 10 5.5 5
+#	v 6 U 10 5.5 6
+#	v 7 U 10 5.5 7
+#	v 8 U 10 5.5 8
+#	v 9 G 10 5.5 9
+#	e 0 1 > 1
+#	e 1 2 > 1
+#	e 2 3 > 1
+#	e 3 4 > 1
+#	e 4 5 > 1
+#	e 5 6 > 1
+#	e 6 7 > 1
+#	e 7 8 > 1
+#	e 8 9 > 1
+#	e 0 8 s -15.0 0.1223
+#	e 1 7 s -15.0 0.1223
+#	e 2 6 s -15.0 0.1223
+#
+#	gi = 9+1 = 10
+#
+#	Window = 10 substructure2 = .(((...))) centre 6.5
+#	v 10 A 10 6.5 2
+#	v 11 C 10 6.5 3
+#	v 12 C 10 6.5 4
+#	v 13 C 10 6.5 5
+#	v 14 U 10 6.5 6
+#	v 15 U 10 6.5 7
+#	v 16 U 10 6.5 8
+#	v 17 G 10 6.5 9
+#	v 18 G 10 6.5 10
+#	v 19 G 10 6.5 11
+#	e 10 11 > 1
+#	e 11 12 > 1
+#	e 12 13 > 1
+#	e 13 14 > 1
+#	e 14 15 > 1
+#	e 15 16 > 1
+#	e 16 17 > 1
+#	e 17 18 > 1
+#	e 18 19 > 1
+#	e 11 19 s -17.0 0.156
+#	e 12 18 s -17.0 0.156
+#	e 13 17 s -17.0 0.156
+#
+#	gi = 19+1 = 20
+#
+#
+#
+#	Write one perl script to create graphs for a set of sequences, fasta2shrep_gspan.pl.
+#
+#	INPUT:
+#		-f fasta file with all sequences to compute
+#		parameters as above
+#
+#	OUTPUT:
+#		one file per sequence that contains graph called seq_id.gspan
+#
+#	(1) for each window size call RNAshapes and write to a tmp file
+#	(2) parse result of RNAshapes (catch RNAshapes error - sequence too long?) check for max shreps.
+#	(3) convert RNAshapes result to subgraph -> write to file (readpipe) look at efficiency and errors
+#	(4) repeat (1) to (3) for each sequence
+
+# ABSTRACT STRUCTURE GRAPH
+# To each shrep graph (s #) add an abstract structure graph with special
+# abstract relations. For example, to add the abstract structure graph to
+# the previous shrep structure with gi numbers from 10 to 19, we first
+# identify the abstract shape, i.e. EL-S-HL and then create nodes (labelled
+# with a given name-space followed the actual node name and separated by a
+# hash) and edges for this graph as follows:
+# NOTE: at the moment we add the adjacent base-pairs to the loop definitions
+# but these can be removed if necessary in the future
+# v 20 abstruct#EL
+# v 21 abstruct#S
+# v 22 abstruct#HL
+# e 20 21 n
+# e 21 22 n
+# v 23 ^R
+# e 20 23 ^r
+# e 23 10 @r
+# e 23 11 @r
+# e 23 19 @r
+# v 24 ^R
+# e 21 24 ^r
+# e 24 11 @r
+# e 24 12 @r
+# e 24 13 @r
+# e 24 17 @r
+# e 24 18 @r
+# e 24 19 @r
+# v 25 ^R
+# e 22 25 ^r
+# e 25 13 @r
+# e 25 14 @r
+# e 25 15 @r
+# e 25 16 @r
+# e 25 17 @r
+
+# gi = 25+1 = 26
+
+# ANOTATION FILE
+# This is a file that labels a sequence region with a given annotation. In one
+# file we can have annotations within different name-spaces, for example target
+# sites predicted with different tools.
+# File format is a tab-delimited file with the following columns:
+# SEQID (same ID as in the fasta file header)
+# i (left positon of region)
+# j (right position of region, if one position i=j)
+# NAMESPACE#LABEL (give each annotation type one name-space and choose a label depending on the task)
+
+# E.g we have 2 different miRNAs and we predict the target-sites with (1) IntaRNA and (2) TargetSearch,
+# then this could be IntaRNA#miR1, IntaRNA#miR2, TargetSearch#miR1, and TargetSearch#miR2.
+# All labels can be used more than once for one or more sequences.
+# all labels with the same namespace and the same sequence ID are grouped into
+# one abstract graph, according to the order of i.
+# Example
+# SEQID 	i 	j 	NAMESPACE#LABEL
+# s1	10 	20	IntaRNA#miR1
+# s1	54	60	IntaRNA#miR2
+# s1	15	25	TargetSearch#miR1
+# s1	54	60 	TargetSearch#miR2
+# s2 ...
+
+# We create connected component per sequence ID per name-space,
+# so that all labels with
+# the same namespace are grouped together (for one sequence) as follows:
+# v 26 IntaRNA#miR1
+# v 27 IntaRNA#miR2
+# e 26 27 n
+# then create R nodes and ^r and @r relations as before (from i to j)
+
+# If a sequence graph option is given, then these graphs (u #) have to
+# be labelled with this annotation
+
diff -r 000000000000 -r b01beb170290 readme.rst
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/readme.rst	Tue Oct 29 11:10:19 2013 -0400
@@ -0,0 +1,48 @@
+Galaxy wrapper for fasta2shrep_gspan
+====================================
+
+This wrapper is copyright 2013 by Björn Grüning.
+
+This is a wrapper for the command line tool of fasta2shrep_gspan.pl.
+
+
+
+
+Installation
+============
+
+The recommended installation is by means of the toolshed_.
+If you need to install it manually here is a short introduction.
+
+.. _toolshed: http://toolshed.g2.bx.psu.edu/view/bgruening/sequence2gspan
+
+
+
+History
+=======
+
+- v0.1: Initial public release
+
+
+=============
+Licence (MIT)
+=============
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+
diff -r 000000000000 -r b01beb170290 repository_dependencies.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/repository_dependencies.xml	Tue Oct 29 11:10:19 2013 -0400
@@ -0,0 +1,4 @@
+<?xml version="1.0"?>
+<repositories description="This requires the EDeN datatypes.">
+    <repository changeset_revision="59b3b6ce10bb" name="eden_toolbox" owner="bgruening" toolshed="http://testtoolshed.g2.bx.psu.edu" />
+</repositories>
diff -r 000000000000 -r b01beb170290 sequence2gspan.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sequence2gspan.xml	Tue Oct 29 11:10:19 2013 -0400
@@ -0,0 +1,167 @@
+<tool id="sequence2gpsan_rnashapes" name="Sequence to gspan" version="0.3">
+    <description>converter using RNAshapes</description>
+    <parallelism method="multi" split_inputs="infile" split_mode="to_size" split_size="10" shared_inputs="" merge_outputs="outfile"></parallelism>
+    <requirements>
+        <requirement type="package" version="2.1.6">rnashapes</requirement>
+    </requirements>
+    <command interpreter="perl">
+
+        fasta2shrep_gspan.pl
+            --fasta $infile
+
+            #if $annotate_infile:
+                -annotate $annotate_infile
+            #end if
+
+            -shift $shift
+            $cue
+            $stack
+            #if $energy_opts.energy_opts_selector == 'energy':
+                -e $energy_opts.e
+            #elif $energy_opts.energy_opts_selector == 'rel_energy':
+                -c $energy_opts.c
+            #end if
+
+            #if $structure_sampling_opts.structure_sampling_opts_selector == 'on':
+                -i $structure_sampling_opts.i
+                -sample_len $structure_sampling_opts.sample_len
+            #end if
+
+            $q
+            -Tp $Tp
+            $seq_graph_win
+            $seq_graph_t
+            $seq_graph_alph
+            $abstr
+            $nostr
+            $vp
+            $ignore_header
+            -M $M
+            $u
+            $r
+
+            -stdout
+
+            > $output
+
+    </command>
+    <inputs>
+        <param name="ifile" type="data" format="fasta" label="Nucleotide Sequence in FASTA format"/>
+
+        <!-- todo repeat tag, normally only 3 windows are specified-->
+        <param name="wins" type="text" value="" label="A list of window sizes to use"
+            help="if none are given, then the entire sequence is taken with no windows. Each window > 1 required! Example: 50,100,200. (-wins)" />
+        <param name="t" type="text" value="3=0" label="The shape type (RNAshapes types from 1-5)"
+            help="With the list format, the shape level can be changed for different window length '4=100' means that shape level 4 is used from length 100nt (window length) The first given length has to be 0! Not continuous given levels are allowed! (-t)" />
+
+
+        <param name="shift" type="integer" value="1" label="The shift of the window, relative to the window size given in percent" 
+            help="Example 20 % of a window 150 would result in a step size of 30 nt. It is a relative parameter, as you can give different window sizes. Default shift of 1 nt. (-shift)">
+            <validator type="in_range" min="0" max="100" />
+        </param>
+        <param name="cue" type="boolean" label="Crop unpaired ends" truevalue="-cue" falsevalue="" checked="false" 
+            help="If you give this flag, then the unpaired ends of each single structure are ignored. E.g. the structure ...(((...))).. becomes just (((...))). (-cue)" />
+
+        <param name="stack" type="boolean" label="Adds stacking information to graphs" truevalue="-stack" falsevalue="" checked="false" 
+            help="This adds an additional vertex (type P) for each pair of stacked base-pairs and four edges (type p) from each of the involved bases to the new vertex. (-stack)" />
+
+        <conditional name="energy_opts">
+            <param name="energy_opts_selector" type="select" label="Energy calculation mode">
+                <option value="rel_energy">Relative energy range, i.e. percentage (%) of MFE energy (RNAshapes)</option>
+                <option value="energy">Energy range in kcal/mol (RNAshapes)</option>
+                <option value="">None</option>
+            </param>
+            <when value="" />
+            <when value="energy">
+                <param name="e" type="float" value="5.0" label="Energy range in kcal/mol (RNAshapes)" help="-e">
+                    <validator type="in_range" min="0.0" />
+                </param>
+            </when>
+            <when value="rel_energy">
+                <param name="c" type="float" value="10" label="Relative energy range, i.e. percentage (%) of MFE energy (RNAshapes)" help="(-c)">
+                    <validator type="in_range" min="0.0" />
+                </param>
+            </when>
+        </conditional>
+        <param name="M" type="integer" value="0" label="Max number of shreps that should be taken per window" help="0 means selects all shreps. (-M)">
+            <validator type="in_range" min="0" />
+        </param>
+        <param name="u" type="boolean" label="Ignore unstable RNAshapes" truevalue="-u" falsevalue="" checked="false" 
+            help="This option filters out closed structures with positive free energy. (-u)" />
+
+        <param name="r" type="boolean" label="Calculate structure probabilities for RNAshapes" help="(-r)" truevalue="-r" falsevalue="" checked="false" />
+
+
+        <conditional name="structure_sampling_opts">
+            <param name="structure_sampling_opts_selector" type="select" label="Structure Sampling">
+                <option value="on">On</option>
+                <option value="off">Off</option>
+            </param>
+            <when value="on">
+                <param name="i" type="integer" value="1" label="Number of sampling iterations" help="(-i)">
+                    <validator type="in_range" min="1" />
+                </param>
+                <param name="smaple_len" type="integer" value="0" label="Sampling is only used for seqs/windows >= given length" help="Default 0 (sample all length). (-sample-len)">
+                    <validator type="in_range" min="0" />
+                </param>
+            </when>
+            <when value="off" />
+        </conditional>
+
+        <param name="q" type="boolean" label="Turn on shape probabilities for RNAshapes, no sampling mode allowed" help="(-q)" truevalue="-q" falsevalue="" checked="false" />
+        <param name="Tp" type="float" value="0.001" label="Filter cutoff for shape probabilities, applied before -M filter!" help="(-Tp)">
+            <validator type="in_range" min="0.0" />
+        </param>
+
+        <param name="seq_graph_win" type="boolean" label="Add for each window a graph which contains no structure" 
+            help="(-seq-graph-win)" 
+            truevalue="-seq-graph-win" falsevalue="" checked="false" />
+        <param name="seq_graph_t" type="boolean" label="add for each 't #' a graph which contains no structure" 
+            help="(-seq-graph-t)" 
+            truevalue="-seq-graph-t" falsevalue="" checked="false" />
+        <param name="seq_graph_alph" type="boolean" label="Change the alphabet of unstructured graphs" 
+            help="(-seq-graph-alph)" 
+            truevalue="-seq-graph-alph" falsevalue="" checked="false" />
+
+        <param format="tabular" name="annotate_infile" type="data" optional="True" label="A file with annotations to be added as abstract graphs" 
+            help="The format is has the following TAB-delimited columns: SEQID, START, END, NAMESPACE#LABEL. Labels with the same name-space and SEQID form connected components, which is a sequence of label vertices ordered by the START position in the sequence."/>
+
+        <param name="abstr" type="boolean" label="Add abstract structure graphs to the single shrep graph instances" 
+            help="(-abstr)" 
+            truevalue="-abstr" falsevalue="" checked="false" />
+        <param name="nostr" type="boolean" label="Calculate no structures, only add sequence information" 
+            help="-seq-graph-win AND/OR -seq-graph-t are required (-nostr)" 
+            truevalue="-nostr" falsevalue="" checked="false" />
+
+        <!--
+            -match-shape    SHAPE
+            all seqs/windows will be constraint folded into that shape via 
+            RNAshapes (if structure is given in another way this struct will be kept), 
+            if this shape is not possible within given energy range, produce a 
+            specific t graph with only one vertex 'X'. By this the instance 
+            becomes very unsimilar to all other graphs (for knn)
+        -->
+        <param name="vp" type="boolean" label="Enable graph computation with viewpoints" 
+            help="svmsgdnspdk will center on those nucleotides that are given via capital letters and ignore those given as lowercase letters (-vp)" 
+            truevalue="-vp" falsevalue="" checked="false" />
+
+        <param name="ignore_header" type="boolean" label="Don't write fasta id part after first space to gspan" help="(-ignore-header)" truevalue="-ignore-header" falsevalue="" checked="false" />
+
+    </inputs>
+    <outputs>
+        <data format="gspan" name="output" label="${tool.name} on ${on_string}" />
+    </outputs>
+    <tests>
+    </tests>
+    <help>
+
+**What it does**
+
+
+
+**References**
+
+Sita et al...
+
+    </help>
+</tool>
diff -r 000000000000 -r b01beb170290 tool_dependencies.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tool_dependencies.xml	Tue Oct 29 11:10:19 2013 -0400
@@ -0,0 +1,6 @@
+<?xml version="1.0"?>
+<tool_dependency>
+    <package name="rnashapes" version="2.1.6">
+        <repository changeset_revision="e9a43f76591a" name="package_rnashapes_2_1_6" owner="rnateam" toolshed="http://testtoolshed.g2.bx.psu.edu" />
+    </package>
+</tool_dependency>