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