binaries/src/probcons/CompareToRef.cc

   1 /////////////////////////////////////////////////////////////////
   2 // CompareToRef.cc
   3 //
   4 // Program for scoring alignments according to the SUM-OF-PAIRS
   5 // or COLUMN score.
   6 /////////////////////////////////////////////////////////////////
   7
   8 #include "SafeVector.h"
   9 #include "MultiSequence.h"
  10 #include <string>
  11 #include <sstream>
  12 #include <iomanip>
  13 #include <iostream>
  14 #include <list>
  15 #include <set>
  16 #include <limits>
  17 #include <cstdio>
  18 #include <cstdlib>
  19 #include <cerrno>
  20 #include <iomanip>
  21
  22 const char CORE_BLOCK = 'h';
  23 typedef pair<int,int> PII;
  24 bool useCoreBlocks = false;
  25 bool useColScore = false;
  26 bool useCaps = false;
  27 bool useBaliAnnot = false;
  28 bool makeAnnot = false;
  29
  30 /////////////////////////////////////////////////////////////////
  31 // Function prototypes
  32 /////////////////////////////////////////////////////////////////
  33
  34 set<PII> ComputePairs (MultiSequence *align, bool isRef);
  35 set<VI> ComputeColumns (MultiSequence *align, bool isRef);
  36 string GetName (string s);
  37 set<int> coreCols;
  38
  39 set<VI> refCols, testCols;
  40 set<PII> refPairs, testPairs;
  41 VI annotation;
  42
  43 /////////////////////////////////////////////////////////////////
  44 // main()
  45 //
  46 // Main program.
  47 /////////////////////////////////////////////////////////////////
  48
  49 int main (int argc, char **argv){
  50
  51   // check arguments
  52   if (argc < 3){
  53     cerr << "Usage: score TEST_ALIGNMENT REFERENCE_ALIGNMENT [BALIBASE_ANNOT_FILE] [-col] [-core] [-caps] [-annot FILENAME]" << endl;
  54     exit (1);
  55   }
  56
  57   // try opening file
  58   FileBuffer infile (argv[1]);
  59
  60   MultiSequence *testAlign;
  61   if (infile.fail()){
  62     cerr << "ERROR: Could not open file '" << argv[1] << "' for reading." << endl;
  63     testAlign = NULL;
  64   }
  65   else {
  66     testAlign = new MultiSequence(); assert (testAlign);
  67     testAlign->LoadMFA (infile);
  68   }
  69   infile.close();
  70
  71   MultiSequence *refAlign = new MultiSequence (string (argv[2])); assert (refAlign);
  72
  73   string outFilename = "";
  74
  75   for (int i = 3; i < argc; i++){
  76     if (strcmp (argv[i], "-core") == 0)
  77       useCoreBlocks = true;
  78     else if (strcmp (argv[i], "-col") == 0)
  79       useColScore = true;
  80     else if (strcmp (argv[i], "-caps") == 0)
  81       useCaps = true;
  82     else if (strcmp (argv[i], "-annot") == 0){
  83       makeAnnot = true;
  84       outFilename = string (argv[++i]);
  85     }
  86     else { // annotation file
  87       useBaliAnnot = true;
  88
  89       ifstream annotFile (argv[i]);
  90       if (annotFile.fail()){
  91         cerr << "ERROR: Could not read BAliBASE annotation file." << endl;
  92         exit (1);
  93       }
  94
  95       SafeVector<int> *indices = refAlign->GetSequence(0)->GetMapping();
  96
  97       char buffer[10000];
  98       while (annotFile.getline (buffer, 10000)){
  99         istringstream ss;
 100         ss.str (string (buffer));
 101
 102         string s;
 103
 104         if ((ss >> s) && s == string ("BPOS")){
 105           while (ss >> s){
 106             int begin=-1, end=-1;
 107             if (sscanf (s.c_str(), "%d=%d", &begin, &end) == 2){
 108               for (int i = (*indices)[begin]; i <= (*indices)[end]; i++)
 109                 coreCols.insert (i);
 110             }
 111           }
 112         }
 113       }
 114
 115       delete indices;
 116
 117       annotFile.close();
 118     }
 119   }
 120
 121   if (useColScore) makeAnnot = false;
 122
 123   if (testAlign){
 124     for (int i = 0; i < testAlign->GetNumSequences(); i++){
 125       bool found = false;
 126
 127       for (int j = 0; !found && j < refAlign->GetNumSequences(); j++){
 128         if (testAlign->GetSequence(i)->GetHeader() == refAlign->GetSequence(j)->GetHeader())
 129           found = true;
 130       }
 131
 132       if (!found){
 133         testAlign->RemoveSequence (i);
 134         i--;
 135       }
 136     }
 137
 138     for (int i = 0; i < refAlign->GetNumSequences(); i++){
 139       bool found = false;
 140
 141       for (int j = 0; !found && j < testAlign->GetNumSequences(); j++){
 142         if (refAlign->GetSequence(i)->GetHeader() == testAlign->GetSequence(j)->GetHeader())
 143           found = true;
 144       }
 145
 146       if (!found){
 147         refAlign->RemoveSequence (i);
 148         i--;
 149       }
 150     }
 151
 152     testAlign->SortByHeader();
 153     refAlign->SortByHeader();
 154   }
 155
 156   int TP = 0;
 157   int TPFN = 0;
 158   int TPFP = 0;
 159   double FD, FM;
 160   if (useColScore){
 161     refCols = ComputeColumns (refAlign, true);
 162     if (testAlign) testCols = ComputeColumns (testAlign, false);
 163     set<VI> colIntersect;
 164     insert_iterator<set<VI> > colIntersectIter (colIntersect, colIntersect.begin());
 165     set_intersection (testCols.begin(), testCols.end(), refCols.begin(), refCols.end(), colIntersectIter);
 166     TP = (int) colIntersect.size();
 167     TPFN = (int) refCols.size();
 168     if (testAlign) TPFP = (int) testCols.size();
 169   }
 170   else {
 171     refPairs = ComputePairs (refAlign, true);
 172     if (testAlign) testPairs = ComputePairs (testAlign, false);
 173     set<PII> pairIntersect;
 174
 175     insert_iterator<set<PII> > pairIntersectIter (pairIntersect, pairIntersect.begin());
 176     set_intersection (testPairs.begin(), testPairs.end(), refPairs.begin(), refPairs.end(), pairIntersectIter);
 177     TP = (int) pairIntersect.size();
 178     TPFN = (int) refPairs.size();
 179     if (testAlign) TPFP = (int) testPairs.size();
 180   }
 181
 182   FD = (double) TP / TPFN;
 183   FM = (double) TP / TPFP;
 184
 185   cout << GetName(string (argv[2])) << " " << TP << " " << TPFN << " " << TPFP << " " << FD << " " << FM << endl;
 186
 187   if (makeAnnot){
 188     ofstream outfile (outFilename.c_str());
 189     for (int i = 0; i < (int) annotation.size(); i++){
 190       outfile << annotation[i] << endl;
 191     }
 192     outfile.close();
 193   }
 194
 195   if (testAlign) delete testAlign;
 196   delete refAlign;
 197 }
 198
 199 int GetOffset (Sequence *testSeq, Sequence *refSeq){
 200   string test = testSeq->GetString();
 201   string ref = refSeq->GetString();
 202
 203   for (int i = 0; i < (int) test.length(); i++) test[i] = toupper(test[i]);
 204   for (int i = 0; i < (int) ref.length(); i++) ref[i] = toupper(ref[i]);
 205
 206   size_t offset = test.find (ref, 0);
 207   if (offset == string::npos){
 208     cerr << "ERROR: Reference string not found in original sequence!" << endl;
 209     cerr << "       test = " << test << endl;
 210     cerr << "       ref = " << ref << endl;
 211     exit (1);
 212   }
 213
 214   cerr << "Offset found: " << offset << endl;
 215
 216   return (int) offset;
 217 }
 218
 219 string GetName (string s){
 220
 221   size_t index1 = s.rfind ('/');
 222   size_t index2 = s.rfind ('.');
 223
 224   if (index1 == string::npos) index1 = 0; else index1++;
 225   if (index2 == string::npos) index2 = s.length();
 226
 227   if (index2 < index1) index2 = s.length();
 228
 229   return s.substr (index1, index2 - index1);
 230 }
 231
 232 bool isCore (char ch, int col){
 233   if (ch == '-') return false;
 234   if (useBaliAnnot){
 235     return coreCols.find (col) != coreCols.end();
 236   }
 237   if (useCaps){
 238     return ch >= 'A' && ch <= 'Z';
 239   }
 240   return ch == CORE_BLOCK;
 241 }
 242
 243 /////////////////////////////////////////////////////////////////
 244 // ComputePairs
 245 //
 246 // Returns the set of all matching pairs.
 247 /////////////////////////////////////////////////////////////////
 248
 249 set<PII> ComputePairs (MultiSequence *align, bool isRef){
 250   int N = align->GetNumSequences();
 251   int L = align->GetSequence(0)->GetLength();
 252
 253   // retrieve all sequence data pointers
 254   SafeVector<SafeVector<char>::iterator> seqs (N);
 255   for (int i = 0; i < N; i++){
 256     seqs[i] = align->GetSequence(i)->GetDataPtr();
 257     assert (align->GetSequence(i)->GetLength() == L);
 258   }
 259
 260   set<PII> ret;
 261   VI ctr(N);
 262
 263   // compute pairs
 264   for (int i = 1; i <= L; i++){
 265
 266     // ctr keeps track of the current position in each sequence
 267     for (int j = 0; j < N; j++){
 268       ctr[j] += (seqs[j][i] != '-');
 269     }
 270
 271     int good = 0;
 272     int ct = 0;
 273
 274     // check for all matching pairs
 275     for (int j = 0; j < N - 1; j++){
 276       for (int k = j + 1; k < N; k++){
 277
 278         // skip if one of the sequences is gapped
 279         if (seqs[j][i] == '-' || seqs[k][i] == '-') continue;
 280
 281         // check for core blocks in the reference sequence
 282         if (isRef && useCoreBlocks)
 283           if (!isCore (seqs[j][i], i) || !isCore (seqs[k][i], i)) continue;
 284
 285         // if all ok, then add pair to list of pairs
 286         pair<int,int> p (10000 * j + ctr[j], 10000 * k + ctr[k]);
 287
 288         // if we're making an annotation, compute annotation statistics
 289         if (makeAnnot && !isRef){
 290           ct++;
 291           if (refPairs.find (p) != refPairs.end()) good++;
 292         }
 293         ret.insert (p);
 294       }
 295     }
 296
 297     // build annotation
 298     if (makeAnnot && !isRef){
 299       annotation.push_back ((ct == 0) ? 0 : 100 * good / ct);
 300     }
 301
 302   }
 303
 304   return ret;
 305 }
 306
 307 /////////////////////////////////////////////////////////////////
 308 // ComputeColumns
 309 //
 310 // Returns the set of all columns.
 311 /////////////////////////////////////////////////////////////////
 312
 313 set<VI> ComputeColumns (MultiSequence *align,  bool isRef){
 314   int N = align->GetNumSequences();
 315   int L = align->GetSequence(0)->GetLength();
 316
 317   // retrieve all sequence data pointers
 318   SafeVector<SafeVector<char>::iterator> seqs (N);
 319   for (int i = 0; i < N; i++){
 320     seqs[i] = align->GetSequence(i)->GetDataPtr();
 321   }
 322
 323   set<VI> ret;
 324   VI ctr(N);
 325
 326   // compute pairs
 327   for (int i = 1; i <= L; i++){
 328
 329     // ctr keeps track of the current position in each sequence
 330     for (int j = 0; j < N; j++){
 331       ctr[j] += (seqs[j][i] != '-');
 332     }
 333
 334     // add column, pick only positions that are matched
 335     SafeVector<int> column (N);
 336     bool useThisColumn = !useCoreBlocks;
 337
 338     for (int j = 0; j < N; j++){
 339       if (isCore (seqs[j][i], i)) useThisColumn = true;
 340       column[j] = (seqs[j][i] == '-') ? -1 : ctr[j];
 341     }
 342
 343     if (useThisColumn || !isRef)
 344       ret.insert (column);
 345   }
 346
 347   return ret;
 348 }