getparents.cpp

   1 //uchime by Robert C. Edgar http://drive5.com/uchime This code is donated to the public domain.\r
   2 \r
   3 #include "myutils.h"\r
   4 #include "chime.h"\r
   5 #include "ultra.h"\r
   6 #include <set>\r
   7 \r
   8 void AddTargets(Ultra &U, const SeqData &Query, set<unsigned> &TargetIndexes);\r
   9 \r
  10 void GetChunkInfo(unsigned L, unsigned &Length, vector<unsigned> &Los)\r
  11         {\r
  12         Los.clear();\r
  13 \r
  14         if (L <= opt_minchunk)\r
  15                 {\r
  16                 Length = L;\r
  17                 Los.push_back(0);\r
  18                 return;\r
  19                 }\r
  20 \r
  21         Length = (L - 1)/opt_chunks + 1;\r
  22         if (Length < opt_minchunk)\r
  23                 Length = opt_minchunk;\r
  24 \r
  25         unsigned Lo = 0;\r
  26         for (;;)\r
  27                 {\r
  28                 if (Lo + Length >= L)\r
  29                         {\r
  30                         Lo = L - Length - 1;\r
  31                         Los.push_back(Lo);\r
  32                         return;\r
  33                         }\r
  34                 Los.push_back(Lo);\r
  35                 Lo += Length;\r
  36                 }\r
  37         }\r
  38 \r
  39 void GetCandidateParents(Ultra &U, const SeqData &QSD, float AbQ,\r
  40   vector<unsigned> &Parents)\r
  41         {\r
  42         Parents.clear();\r
  43 \r
  44         set<unsigned> TargetIndexes;\r
  45 \r
  46         unsigned QL = QSD.L;\r
  47 \r
  48         SeqData QuerySD = QSD;\r
  49 \r
  50         unsigned ChunkLength;\r
  51         vector<unsigned> ChunkLos;\r
  52         GetChunkInfo(QL, ChunkLength, ChunkLos);\r
  53         unsigned ChunkCount = SIZE(ChunkLos);\r
  54         for (unsigned ChunkIndex = 0; ChunkIndex < ChunkCount; ++ChunkIndex)\r
  55                 {\r
  56                 unsigned Lo = ChunkLos[ChunkIndex];\r
  57                 asserta(Lo + ChunkLength <= QL);\r
  58 \r
  59                 const byte *Chunk = QSD.Seq + Lo;\r
  60 \r
  61         // THIS MESSES UP --self!!\r
  62                 //char Prefix[32];\r
  63                 //sprintf(Prefix, "%u|", Lo);\r
  64                 //string ChunkLabel = string(Prefix) + string(QSD.Label);\r
  65 \r
  66                 //QuerySD.Label = ChunkLabel.c_str();\r
  67                 QuerySD.Seq = Chunk;\r
  68                 QuerySD.L = ChunkLength;\r
  69 \r
  70                 AddTargets(U, QuerySD, TargetIndexes);\r
  71 \r
  72                 Lo += ChunkLength;\r
  73                 }\r
  74 \r
  75         for (set<unsigned>::const_iterator p = TargetIndexes.begin();\r
  76           p != TargetIndexes.end(); ++p)\r
  77                 {\r
  78                 unsigned TargetIndex = *p;\r
  79                 bool Accept = true;\r
  80                 if (AbQ > 0.0f)\r
  81                         {\r
  82                         const char *TargetLabel = U.GetSeedLabel(TargetIndex);\r
  83                         float AbT = GetAbFromLabel(string(TargetLabel));\r
  84                         if (AbT > 0.0f && AbT < opt_abskew*AbQ)\r
  85                                 Accept = false;\r
  86                         }\r
  87 \r
  88                 if (Accept)\r
  89                         Parents.push_back(TargetIndex);\r
  90                 }\r
  91         }\r