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