1#!/usr/bin/perl -I.libs
2
3######################### We start with some black magic to print on failure.
4# (It may become useful if the test is moved to ./t subdirectory.)
5use strict;
6use Test::More tests => 4;
7
8use RNA;
9use warnings;
10use RNAHelpers qw(:Paths :Messages);
11
12
13my $datadir = getDataDirPath();
14
15my $seq_con     = "CCCAAAAGGGCCCAAAAGGG";
16my $str_con     = "..........(((....)))";
17my $str_con_def = "(((....)))(((....)))";
18my $seq_long    = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUCAUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUCAGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUC";
19my $fc; #fold compound reference
20my $ss=""; #return string
21my $mfe=0;
22my $energy=0;
23my @data = ();
24my $d = 0;
25my $failed = 0;
26my $e   = 0;
27my $s = "";
28my $c = 0;
29
30sub mfe_window_callback {
31
32  my ($start, $end, $structure, $energy, $data) = @_;
33  my %Lfold_hit = ();
34  $Lfold_hit{'structure'} = $structure;
35  $Lfold_hit{'start'}     = $start;
36  $Lfold_hit{'end'}       = $end;
37  $Lfold_hit{'energy'}    = $energy;
38
39  push @{$data}, \%Lfold_hit;
40}
41
42
43
44##################################
45##test_sc_set_up
46
47MsgChecking("whether setting unpaired soft constraints to all nucleotides at once works");
48my $seq_sc  =      "CCCAAAAGGG";
49$fc         = new RNA::fold_compound($seq_sc);
50($ss, $mfe) = $fc->mfe();
51printf("%s [%6.2f] \n",$ss,$mfe);
52is($ss,"(((....)))");
53
54$fc->sc_init();
55
56my @m= (0.0,-5.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);  #E 1 0 1 -5 ,  position 1 gets -5 if unpaired , vector starts with 0 and not 1
57
58$fc->sc_set_up(\@m);
59($ss,$mfe) = $fc->mfe();
60printf("%s [%6.2f] \n", $ss, $mfe);
61is(sprintf ("%6.2f", $mfe), sprintf("%6.2f", -5.70));
62undef @m;
63##################################   We have no argin typemap for multidimensional lists yet, so this doesn't work
64#test_sc_set_bp
65#MsgChecking("whether setting base pair soft constraints to all nucleotides at once works");
66#
67#@m = ();
68#push @m, [(0) x (length($seq_sc) + 1)] for 0 .. length($seq_sc);
69#
70# add energy of -5 to basepair 1-9 if formed, prefered structure should now be ((.....))., with a energy of -4.90
71#$m[1][9] = -5.0; # base 1-9 should get -5.0 if basepair
72#$m[9][1] = -5.0; # base 1-9 should get -5.0 if basepair
73#
74#$fc = new RNA::fold_compound($seq_sc);
75#$fc->sc_set_bp(\@m);
76#($ss, $mfe) = $fc->mfe();
77#print("%s [%6.2f]\n", $ss, $mfe);
78#is(sprintf ("%6.2f", $mfe), sprintf("%6.2f", -4.90));
79
80
81##################################
82##test_sc_mfe_window_add_bp
83
84MsgChecking("whether adding base pair soft constraints works in conjunction with sliding-window MFE prediction");
85
86$fc = new RNA::fold_compound($seq_long, undef, RNA::OPTION_WINDOW);
87
88# add twice -10.0 kcal/mol on base pair (55,60)
89$fc->sc_add_bp(55, 60, -10.0, RNA::OPTION_WINDOW);
90$fc->sc_add_bp(55, 60, -10.0, RNA::OPTION_WINDOW);
91
92# allow pair (55,60) to actually form (might be non-canonical)
93$fc->hc_add_bp(55, 60, RNA::CONSTRAINT_CONTEXT_NO_REMOVE | RNA::CONSTRAINT_CONTEXT_ALL_LOOPS);
94
95@data = ();
96$mfe = $fc->mfe_window_cb(\&mfe_window_callback, \@data);
97
98$failed = 0;
99foreach my $hit (@data) {
100    if ($hit->{'start'} <= 55 && $hit->{'end'} >= 60) {
101        # compose actual dot bracket string (including potential 5' dangle nucleotide
102        if ($hit->{'start'} > 1) {
103          $d = 1;
104          $ss = ".".$hit->{'structure'};
105        } else {
106          $d = 0;
107          $ss = $hit->{'structure'};
108        }
109
110        # get corresponding subsequence
111        $s = substr($seq_long, $hit->{'start'} - 1 - $d, $hit->{'end'} - $hit->{'start'} + 2);
112
113        # re-evaluate free energy of subsequence/hit
114        $e = RNA::energy_of_struct($s, $ss);
115
116        # energy difference between both must be -20.0 kcal/mol (if the constrained base pair is present)
117        $failed = 1 if ! (sprintf("%6.2f", $hit->{'energy'}) eq sprintf("%6.2f", $e - 20.0));
118    }
119}
120ok($failed == 0);
121
122##################################
123##test_sc_mfe_window_add_bp
124
125MsgChecking("whether adding unpaired soft constraints works in conjunction with sliding-window MFE prediction");
126
127$fc = new RNA::fold_compound($seq_long, undef, RNA::OPTION_WINDOW);
128
129# add twice -5.0 kcal/mol per unpaired nucleotide in segment [55,60]
130for my $i (50..60) {
131    $fc->sc_add_up($i, -5.0, RNA::OPTION_WINDOW);
132}
133
134# force segment [50,60] to stay unpaired
135for my $i (50..60) {
136    $fc->hc_add_up($i, RNA::CONSTRAINT_CONTEXT_ALL_LOOPS);
137}
138
139@data = ();
140$mfe = $fc->mfe_window_cb(\&mfe_window_callback, \@data);
141
142$failed = 0;
143
144foreach my $hit (@data) {
145    if ($hit->{'start'} <= 55 && $hit->{'end'} >= 60) {
146        # count number of unpiared nucleotides with bonus
147        $c = 0;
148        foreach my $i (50..60) {
149            $c++ if substr($hit->{'structure'}, $i - $hit->{'start'}, 1) eq ".";
150        }
151
152        # compose actual dot bracket string (including potential 5' dangle nucleotide
153        if ($hit->{'start'} > 1) {
154          $d = 1;
155          $ss = ".".$hit->{'structure'};
156        } else {
157          $d = 0;
158          $ss = $hit->{'structure'};
159        }
160
161        # get corresponding subsequence
162        $s = substr($seq_long, $hit->{'start'} - 1 - $d, $hit->{'end'} - $hit->{'start'} + 1 + $d);
163
164        # re-evaluate free energy of subsequence/hit
165        $e = RNA::energy_of_struct($s, $ss);
166
167        # energy difference between both must be c * -5.0 kcal/mol
168        $failed = 1 if ! (sprintf("%6.2f", $hit->{'energy'}) eq sprintf("%6.2f", $e + ($c * -5.0)));
169    }
170}
171ok($failed == 0);
172
173undef $fc;
174