t/psm.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib 't/lib';
   8     use BioperlTest;
   9
  10     test_begin(-tests => 63);
  11
  12         use_ok ('Bio::Matrix::PSM::IO');
  13 }
  14
  15 my $mmt= "chr04q        170164  170208  strong  -       0       Motif 3 occurrance in chr04q
  16 chr04q  215755  215799  strong  +       0       Motif 4 occurrance in chr04q
  17 chr04q  532530  532574  strong  +       2       Motif 2 occurrance in chr04q
  18 chr04q  539492  539536  strong  -       1       Motif 1 occurrance in chr04q
  19 chr04q  586113  586157  strong  +       2       Motif 2 occurrance in chr04q
  20 chr04q  698245  698289  strong  -       0       Motif 4 occurrance in chr04q
  21 chr04q  804412  804456  strong  -       0       Motif 3 occurrance in chr04q
  22 chr04q  858870  858914  strong  -       2       Motif 3 occurrance in chr04q
  23 chr04q  861561  861605  strong  -       2       Motif 3 occurrance in chr04q
  24 chr04q  916898  916942  strong  -       1       Motif 1 occurrance in chr04q
  25 chr04q  1146916 1146960 strong  -       0       Motif 1 occurrance in chr04q
  26 chr04q  1315772 1315816 strong  +       1       Motif 1 occurrance in chr04q
  27 chr04q  1636119 1636163 strong  +       2       Motif 3 occurrance in chr04q
  28 chr04q  1636200 1636244 strong  +       2       Motif 1 occurrance in chr04q
  29 chr04q  1636437 1636481 strong  +       2       Motif 4 occurrance in chr04q
  30 chr04q  1637361 1637405 strong  +       2       Motif 2 occurrance in chr04q
  31 chr04q  1652447 1652491 strong  +       1       Motif 4 occurrance in chr04q";
  32 my @mmt=split(/\n/,$mmt);
  33
  34 #Let's try meme here
  35 my $psmIO =  Bio::Matrix::PSM::IO->new(-format=>'meme',
  36              -file=>test_input_file('meme.dat'));
  37 ok $psmIO;
  38
  39 my @inputfile=grep(/datafile/i,$psmIO->unstructured);
  40 ok @inputfile;
  41
  42 my $release=$psmIO->release;
  43 ok $release;
  44
  45 my @ids=$psmIO->hid;
  46 is @ids,4;
  47
  48 my %weights=$psmIO->weight;
  49 ok %weights;
  50
  51 my %seq = $psmIO->seq;
  52 is %seq,'0';#Meme doesn't have seq
  53
  54 is $psmIO->version,'3.0';
  55
  56 my $psm = $psmIO->next_psm;
  57 ok $psm;
  58
  59 #Lets try to compress and uncompress the log odds and the frequencies, see if there is no
  60 #considerable loss of data.
  61 my $fA=$psm->get_compressed_freq('A');
  62 my @check=Bio::Matrix::PSM::SiteMatrix::_uncompress_string($fA,1,1);
  63 my @A=$psm->get_array('A');
  64 my ($var,$max) = (0,0);
  65 for (my $i = 0; $i<@check;$i++) {
  66   my $diff=abs(abs($check[$i])-abs($A[$i]));
  67   $var += $diff;
  68   $max=$diff if ($diff>$max);
  69 }
  70 my $avg=$var/@check;
  71 cmp_ok $avg,'<',0.01; #Loss of data under 1 percent
  72 #print $avg,"\n";
  73 is $psm->sequence_match_weight('CAGAAAAATAAAATGGCCACCACCC'),2015;
  74
  75 my $lA=$psm->get_compressed_logs('A');
  76 @check=Bio::Matrix::PSM::SiteMatrix::_uncompress_string($lA,1000,2);
  77 @A=$psm->get_logs_array('A');
  78 ($var,$max) = (0,0);
  79 for (my $i = 0;$i<@check;$i++) {
  80   my $diff=abs(abs($check[$i])-abs($A[$i]));
  81   $var += $diff;
  82   $max=$diff if ($diff>$max);
  83 }
  84 $avg=$var/@check;
  85 cmp_ok $avg,'<',10; #Loss of data under 1 percent
  86
  87 my $matrix=$psm->matrix;
  88 ok $matrix;
  89 my $psm2=$psm;
  90 $psm2->matrix($matrix);
  91 is $psm,$psm2;
  92
  93 my %psm_header=$psm->header;
  94 is $psm_header{IC},38.1;
  95 is $psm_header{sites},4;
  96 is $psm_header{width},25;
  97 is $psm_header{e_val},'1.2e-002';
  98
  99
 100 #Quick check if returned object works
 101 my $IUPAC=$psm->IUPAC;
 102 is $IUPAC,'CMKWMAAAKWVAWTYCMCASCHCCM';
 103 is $IUPAC,$psm2->IUPAC;
 104 is $IUPAC,$matrix->IUPAC;
 105
 106 my $instances=$psm->instances;
 107 ok $instances;
 108
 109 foreach my $instance (@{$instances}) {
 110   my $id=$instance->primary_id;
 111   is $instance->strand,1;
 112   last if (ok $id);
 113 }
 114
 115 ok $psm->header('e_val');
 116 #Meme parser should be OK if tests passed
 117
 118
 119 #Now we are going to try transfac
 120
 121 $psmIO =  Bio::Matrix::PSM::IO->new(-format=>'transfac',
 122           -file=> test_input_file('transfac.dat'));
 123 ok $psmIO;
 124
 125 my $version=$psmIO->version;
 126 ok !$version;
 127
 128 is $psmIO->release, '6.4--2002-12-02';
 129
 130 $psm     = $psmIO->next_psm;
 131 ok $psm;
 132
 133 # Lets try to compress and uncompress the the frequencies, see if
 134 # there is no considerable loss of data.
 135 $fA=$psm->get_compressed_freq('A');
 136 @check=Bio::Matrix::PSM::SiteMatrix::_uncompress_string($fA,1,1);
 137 @A=$psm->get_array('A');
 138 ($var,$max) = (0,0);
 139 for (my $i = 0; $i<@check;$i++) {
 140   my $diff=abs(abs($check[$i])-abs($A[$i]));
 141   $var += $diff;
 142   $max=$diff if ($diff>$max);
 143 }
 144 $avg=$var/@check;
 145 cmp_ok $avg,'<',0.01; #Loss of data under 1 percent
 146
 147 %weights = $psmIO->weight;
 148 ok !$weights{''};
 149
 150 %seq     = $psmIO->seq;
 151 is scalar keys %seq, 0;
 152
 153 #Quick check if returned object works
 154 $IUPAC   = $psm->IUPAC;
 155 is $IUPAC,'VVDCAKSTGBYD';
 156
 157 #Now we are going to try mast
 158 $psmIO =  Bio::Matrix::PSM::IO->new(-format=>'mast',
 159           -file=>test_input_file('mast.dat'));
 160 ok $psmIO;
 161
 162 @inputfile = grep(/datafile/i,$psmIO->unstructured);
 163 ok !@inputfile;
 164
 165 is( $psmIO->release, '2002/04/02 0:11:59');
 166
 167 @ids     = $psmIO->hid;
 168 is @ids,4;
 169
 170 %weights = $psmIO->weight;
 171 ok !%weights; #Mast doesn't have weights
 172
 173 ok %seq    = $psmIO->seq;
 174
 175 foreach my $id ($psmIO->hid) {
 176     ok $seq{$id};
 177 }
 178 ok $psm=$psmIO->next_psm;
 179
 180 my %instances=$psmIO->instances;
 181 ok %instances;
 182
 183 is $psmIO->version, '3.0';
 184
 185 my $mmastIO=Bio::Matrix::PSM::IO->new(-format=>'mast',-file=>test_input_file('mixedmast.dat'));
 186
 187 $psm = $mmastIO->next_psm;
 188 my $lastinstances = $psm->instances();
 189 my $i=0;
 190 foreach my $hit (@$lastinstances) {
 191     $hit -> end ( $hit-> start () + length ($hit->seq) - 1 ) ; # fix an old bug in InstanceSite.pm
 192     my $d=join("\t",$hit->{accession_number},$hit -> start () , $hit-> end (),$hit -> score (),
 193     $hit -> strand == 1 ? '+' : '-' , $hit -> frame,  $hit -> desc ( ));
 194     is $d,$mmt[$i];
 195     $i++;
 196     last if ($hit -> start == 1652447);
 197 }