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;
7use Data::Dumper;
8use FileHandle;
9
10use RNA;
11use warnings;
12use RNAHelpers qw(:Paths);
13
14
15my $datadir = getDataDirPath();
16
17my $seq_con     = "CCCAAAAGGGCCCAAAAGGG";
18my $str_con     = "..........(((....)))";
19my $str_con_def = "(((....)))(((....)))";
20my $fc; #fold compound reference
21my $ss=""; #return string
22my $mfe=0;
23my $energy=0;
24
25
26
27sub getShapeDataFromFile
28{
29    my $filepath = shift(@_);
30    my @retVec;
31    push(@retVec,-999.0);# data list is 1-based, so we add smth. at pos 0
32    my $count=1;
33    open(my $fh, "<", $filepath) || die "Couldn't open '".$filepath."' for reading because: ".$!;
34
35    foreach my $line (<$fh>)
36    {
37
38        my @sp= split(/ /,$line);
39        my $pos = $sp[0];
40        my $value = $sp[1];
41
42        if($pos == $count)
43        {
44            push(@retVec,$value +0.00); # make float from string
45        }else
46        {
47            for(my $i=$pos;$i<=$count;$i++)
48            {
49                push(@retVec,-999.0);
50            }
51            push(@retVec,$value);
52            $count=$pos;
53        }
54        $count+=1;
55
56    }
57    return @retVec;
58}
59
60sub getShapeSequenceFromFile
61{
62    my $f = shift(@_);
63    my $retSeq ="";
64    open(my $fh, "<", $f) || die "Couldn't open '".$f."' for reading because: ".$!;
65    foreach my $line (<$fh>)
66    {
67        chomp $line;
68        return $line; #return the first line
69    }
70}
71
72##################################
73##test_sc_add_deigan
74print "test_sc_add_deigan\n";
75my $seq  =  getShapeSequenceFromFile($datadir . "Lysine_riboswitch_T._martima.db");
76my @reactivities = getShapeDataFromFile($datadir . "Lysine_riboswitch_T._martima.shape_2rows");
77
78$fc= new RNA::fold_compound($seq);
79print "@reactivities \n";
80
81$fc->sc_add_SHAPE_deigan(\@reactivities,1.9,-0.7,RNA::OPTION_DEFAULT);
82($ss,$mfe) = $fc->mfe();
83printf("%s [%6.2f] \n",$ss,$mfe);
84is(sprintf ("%6.2f",$mfe), sprintf("%6.2f",-121.55));
85##################################
86##test_sc_add_SHAPE_deigan2
87print "test_sc_add_SHAPE_deigan2";
88my $seq_ribo  =  getShapeSequenceFromFile($datadir . "TPP_riboswitch_E.coli.db");
89my @reactivities_ribo = getShapeDataFromFile($datadir . "TPP_riboswitch_E.coli.shape_2rows");
90
91$fc= new RNA::fold_compound($seq_ribo);
92print "@reactivities_ribo \n";
93
94$fc->sc_add_SHAPE_deigan(\@reactivities_ribo,1.9,-0.7,RNA::OPTION_DEFAULT);
95($ss,$mfe) = $fc->mfe();
96printf("%s [%6.2f] \n",$ss,$mfe);
97is(sprintf ("%6.2f",$mfe), sprintf("%6.2f",-52.61));
98##################################
99# I just added completely randomly choosen sequences and shape data, this test only checks if the function can be called, not if it results correct results.
100##test_sc_add_SHAPE_deigan_ali
101print "test_sc_add_SHAPE_deigan_ali";
102# shape data from
103my $shapeSeq1 = "CCCAAAAGGG";
104my $shapeSeq2 = "AAUAAAAAUU";
105
106#my @shapeData1 = (-999.0,0.04,1.12,1.1,-999.0,0.05,0.5,0.3,-999.0,1.4);
107#my @shapeData2 = (-999.0,-999.0,-999.0,1.23,1.4,0.05,0.5,0.3,-999.0,1.4);
108my @shapeAli = ($shapeSeq1,$shapeSeq2);
109$fc= new RNA::fold_compound(\@shapeAli);
110
111my @assoc = (-1,1,2);
112my $ret = $fc->sc_add_SHAPE_deigan_ali(\@shapeAli, \@assoc,1.8,-0.6);
113($ss,$mfe) = $fc->mfe();
114printf("%s [%6.2f] \n",$ss,$mfe);
115is($ret,1);
116##################################
117##test_sc_add_SHAPE_zarringhalam
118print("test_sc_add_SHAPE_zarringhalam");
119$seq_ribo  =  getShapeSequenceFromFile($datadir . "TPP_riboswitch_E.coli.db");
120@reactivities_ribo = getShapeDataFromFile($datadir . "TPP_riboswitch_E.coli.shape_2rows");
121
122$fc= new RNA::fold_compound($seq_ribo);
123print "@reactivities_ribo \n";
124
125$fc->sc_add_SHAPE_zarringhalam(\@reactivities_ribo,0.5,0.5,"O");
126($ss,$mfe) = $fc->mfe();
127printf("%s [%6.2f] \n",$ss,$mfe);
128is(sprintf ("%6.2f",$mfe), sprintf("%6.2f",-5.28));
129
130undef $fc;
131