uchime_src/writechhit.cpp

   1 #include "myutils.h"\r
   2 #include "chime.h"\r
   3 \r
   4 void WriteChimeFileHdr(FILE *f)\r
   5         {\r
   6         if (f == 0)\r
   7                 return;\r
   8 \r
   9         fprintf(f,\r
  10                 "\tQuery"               // 1\r
  11                 "\tA"                   // 2\r
  12                 "\tB"                   // 3\r
  13                 "\tIdQM"                // 4\r
  14                 "\tIdQA"                // 5\r
  15                 "\tIdQB"                // 6\r
  16                 "\tIdAB"                // 7\r
  17                 "\tIdQT"                // 8\r
  18                 "\tLY"                  // 9\r
  19                 "\tLN"                  // 10\r
  20                 "\tLA"                  // 11\r
  21                 "\tRY"                  // 12\r
  22                 "\tRN"                  // 13\r
  23                 "\tRA"                  // 14\r
  24                 "\tDiv"                 // 15\r
  25                 "\tY"                   // 16\r
  26                 "\n"\r
  27                 );\r
  28         }\r
  29 \r
  30 void WriteChimeHit(FILE *f, const ChimeHit2 &Hit)\r
  31         {\r
  32         if (f == 0)\r
  33                 return;\r
  34 \r
  35         if (Hit.Div <= 0.0)\r
  36                 {\r
  37                 fprintf(f, "0.0000");           // 0\r
  38 \r
  39                 fprintf(f,\r
  40                   "\t%s", Hit.QLabel.c_str());  // 1\r
  41 \r
  42                 fprintf(f,\r
  43                   "\t*"                                         // 2\r
  44                   "\t*"                                         // 3\r
  45                   "\t*"                                         // 4\r
  46                   "\t*"                                         // 5\r
  47                   "\t*"                                         // 6\r
  48                   "\t*"                                         // 7\r
  49                   "\t*"                                         // 8\r
  50                   "\t*"                                         // 9\r
  51                   "\t*"                                         // 10\r
  52                   "\t*"                                         // 11\r
  53                   "\t*"                                         // 12\r
  54                   "\t*"                                         // 13\r
  55                   "\t*"                                         // 14\r
  56                   "\t*"                                         // 15\r
  57                   "\tN"                                         // 16\r
  58                   "\n"\r
  59                   );\r
  60                 return;\r
  61                 }\r
  62 \r
  63         fprintf(f, "%.4f", Hit.Score);          // 0\r
  64 \r
  65         fputc('\t', f);\r
  66         fputs(Hit.QLabel.c_str(), f);           // 1\r
  67 \r
  68         fputc('\t', f);\r
  69         fputs(Hit.ALabel.c_str(), f);           // 2\r
  70 \r
  71         fputc('\t', f);\r
  72         fputs(Hit.BLabel.c_str(), f);           // 3\r
  73 \r
  74         fprintf(f, "\t%.1f", Hit.PctIdQM);      // 4\r
  75         fprintf(f, "\t%.1f", Hit.PctIdQA);      // 5\r
  76         fprintf(f, "\t%.1f", Hit.PctIdQB);      // 6\r
  77         fprintf(f, "\t%.1f", Hit.PctIdAB);      // 7\r
  78         fprintf(f, "\t%.1f", Hit.PctIdQT);      // 8\r
  79 \r
  80         fprintf(f, "\t%u", Hit.CS_LY);          // 9\r
  81         fprintf(f, "\t%u", Hit.CS_LN);          // 10\r
  82         fprintf(f, "\t%u", Hit.CS_LA);          // 11\r
  83 \r
  84         fprintf(f, "\t%u", Hit.CS_RY);          // 12\r
  85         fprintf(f, "\t%u", Hit.CS_RN);          // 13\r
  86         fprintf(f, "\t%u", Hit.CS_RA);          // 14\r
  87 \r
  88         fprintf(f, "\t%.2f", Hit.Div);          // 15\r
  89 \r
  90         fprintf(f, "\t%c", yon(Hit.Accept())); // 16\r
  91         fputc('\n', f);\r
  92         }\r
  93 \r
  94 unsigned GetUngappedLength(const byte *Seq, unsigned L)\r
  95         {\r
  96         unsigned UL = 0;\r
  97         for (unsigned i = 0; i < L; ++i)\r
  98                 if (!isgap(Seq[i]))\r
  99                         ++UL;\r
 100         return UL;\r
 101         }\r
 102 \r
 103 void WriteChimeHitX(FILE *f, const ChimeHit2 &Hit)\r
 104         {\r
 105         if (f == 0)\r
 106                 return;\r
 107 \r
 108         if (Hit.Div <= 0.0)\r
 109                 return;\r
 110 \r
 111         const string &Q3 = Hit.Q3;\r
 112         const string &A3 = Hit.A3;\r
 113         const string &B3 = Hit.B3;\r
 114 \r
 115         const byte *Q3Seq = (const byte *) Q3.c_str();\r
 116         const byte *A3Seq = (const byte *) A3.c_str();\r
 117         const byte *B3Seq = (const byte *) B3.c_str();\r
 118 \r
 119 // Aligned\r
 120         unsigned ColCount = SIZE(Q3);\r
 121         asserta(SIZE(A3) == ColCount && SIZE(B3) == ColCount);\r
 122 \r
 123         unsigned LQ = GetUngappedLength(Q3Seq, ColCount);\r
 124         unsigned LA = GetUngappedLength(A3Seq, ColCount);\r
 125         unsigned LB = GetUngappedLength(B3Seq, ColCount);\r
 126 \r
 127         fprintf(f, "\n");\r
 128         fprintf(f, "------------------------------------------------------------------------\n");\r
 129         fprintf(f, "Query   (%5u nt) %s\n", LQ, Hit.QLabel.c_str());\r
 130         fprintf(f, "ParentA (%5u nt) %s\n", LA, Hit.ALabel.c_str());\r
 131         fprintf(f, "ParentB (%5u nt) %s\n", LB, Hit.BLabel.c_str());\r
 132 \r
 133 // Strip terminal gaps in query\r
 134         unsigned FromCol = UINT_MAX;\r
 135         unsigned ToCol = UINT_MAX;\r
 136         for (unsigned Col = 0; Col < ColCount; ++Col)\r
 137                 {\r
 138                 if (!isgap(Q3Seq[Col]))\r
 139                         {\r
 140                         if (FromCol == UINT_MAX)\r
 141                                 FromCol = Col;\r
 142                         ToCol = Col;\r
 143                         }\r
 144                 }\r
 145 \r
 146         unsigned QPos = 0;\r
 147         unsigned APos = 0;\r
 148         unsigned BPos = 0;\r
 149         for (unsigned Col = 0; Col < FromCol; ++Col)\r
 150                 {\r
 151                 if (!isgap(A3Seq[Col]))\r
 152                         ++APos;\r
 153                 if (!isgap(B3Seq[Col]))\r
 154                         ++BPos;\r
 155                 }\r
 156 \r
 157         unsigned Range = ToCol - FromCol + 1;\r
 158         unsigned RowCount = (Range + 79)/80;\r
 159         unsigned RowFromCol = FromCol;\r
 160         for (unsigned RowIndex = 0; RowIndex < RowCount; ++RowIndex)\r
 161                 {\r
 162                 fprintf(f, "\n");\r
 163                 unsigned RowToCol = RowFromCol + 79;\r
 164                 if (RowToCol > ToCol)\r
 165                         RowToCol = ToCol;\r
 166 \r
 167         // A row\r
 168                 fprintf(f, "A %5u ", APos + 1);\r
 169                 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)\r
 170                         {\r
 171                         char q = Q3Seq[Col];\r
 172                         char a = A3Seq[Col];\r
 173                         if (a != q)\r
 174                                 a = tolower(a);\r
 175                         fprintf(f, "%c", a);\r
 176                         if (!isgap(a))\r
 177                                 ++APos;\r
 178                         }\r
 179                 fprintf(f, " %u\n", APos);\r
 180 \r
 181         // Q row\r
 182                 fprintf(f, "Q %5u ", QPos + 1);\r
 183                 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)\r
 184                         {\r
 185                         char q = Q3Seq[Col];\r
 186                         fprintf(f, "%c", q);\r
 187                         if (!isgap(q))\r
 188                                 ++QPos;\r
 189                         }\r
 190                 fprintf(f, " %u\n", QPos);\r
 191 \r
 192         // B row\r
 193                 fprintf(f, "B %5u ", BPos + 1);\r
 194                 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)\r
 195                         {\r
 196                         char q = Q3Seq[Col];\r
 197                         char b = B3Seq[Col];\r
 198                         if (b != q)\r
 199                                 b = tolower(b);\r
 200                         fprintf(f, "%c", b);\r
 201                         if (!isgap(b))\r
 202                                 ++BPos;\r
 203                         }\r
 204                 fprintf(f, " %u\n", BPos);\r
 205 \r
 206         // Diffs\r
 207                 fprintf(f, "Diffs   ");\r
 208                 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)\r
 209                         {\r
 210                         char q = Q3Seq[Col];\r
 211                         char a = A3Seq[Col];\r
 212                         char b = B3Seq[Col];\r
 213 \r
 214                         char c = ' ';\r
 215                         if (isgap(q) || isgap(a) || isgap(b))\r
 216                                 c = ' ';\r
 217                         else if (Col < Hit.ColXLo)\r
 218                                 {\r
 219                                 if (q == a && q == b)\r
 220                                         c = ' ';\r
 221                                 else if (q == a && q != b)\r
 222                                         c = 'A';\r
 223                                 else if (q == b && q != a)\r
 224                                         c = 'b';\r
 225                                 else if (a == b && q != a)\r
 226                                         c = 'N';\r
 227                                 else\r
 228                                         c = '?';\r
 229                                 }\r
 230                         else if (Col > Hit.ColXHi)\r
 231                                 {\r
 232                                 if (q == a && q == b)\r
 233                                         c = ' ';\r
 234                                 else if (q == b && q != a)\r
 235                                         c = 'B';\r
 236                                 else if (q == a && q != b)\r
 237                                         c = 'a';\r
 238                                 else if (a == b && q != a)\r
 239                                         c = 'N';\r
 240                                 else\r
 241                                         c = '?';\r
 242                                 }\r
 243 \r
 244                         fprintf(f, "%c", c);\r
 245                         }\r
 246                 fprintf(f, "\n");\r
 247 \r
 248         // SNPs\r
 249                 fprintf(f, "Votes   ");\r
 250                 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)\r
 251                         {\r
 252                         char q = Q3Seq[Col];\r
 253                         char a = A3Seq[Col];\r
 254                         char b = B3Seq[Col];\r
 255 \r
 256                         bool PrevGap = Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1]));\r
 257                         bool NextGap = Col+1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1]));\r
 258 \r
 259                         char c = ' ';\r
 260                         if (isgap(q) || isgap(a) || isgap(b) || PrevGap || NextGap)\r
 261                                 c = ' ';\r
 262                         else if (Col < Hit.ColXLo)\r
 263                                 {\r
 264                                 if (q == a && q == b)\r
 265                                         c = ' ';\r
 266                                 else if (q == a && q != b)\r
 267                                         c = '+';\r
 268                                 else if (q == b && q != a)\r
 269                                         c = '!';\r
 270                                 else\r
 271                                         c = '0';\r
 272                                 }\r
 273                         else if (Col > Hit.ColXHi)\r
 274                                 {\r
 275                                 if (q == a && q == b)\r
 276                                         c = ' ';\r
 277                                 else if (q == b && q != a)\r
 278                                         c = '+';\r
 279                                 else if (q == a && q != b)\r
 280                                         c = '!';\r
 281                                 else\r
 282                                         c = '0';\r
 283                                 }\r
 284 \r
 285                         fprintf(f, "%c", c);\r
 286                         }\r
 287                 fprintf(f, "\n");\r
 288 \r
 289         // LR row\r
 290                 fprintf(f, "Model   ");\r
 291                 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)\r
 292                         {\r
 293                         if (Col < Hit.ColXLo)\r
 294                                 fprintf(f, "A");\r
 295                         else if (Col >= Hit.ColXLo && Col <= Hit.ColXHi)\r
 296                                 fprintf(f, "x");\r
 297                         else\r
 298                                 fprintf(f, "B");\r
 299                         }\r
 300 \r
 301                 fprintf(f, "\n");\r
 302 \r
 303                 RowFromCol += 80;\r
 304                 }\r
 305         fprintf(f, "\n");\r
 306 \r
 307         double PctIdBestP = max(Hit.PctIdQA, Hit.PctIdQB);\r
 308         double Div = (Hit.PctIdQM - PctIdBestP)*100.0/PctIdBestP;\r
 309 \r
 310         unsigned LTot = Hit.CS_LY + Hit.CS_LN + Hit.CS_LA;\r
 311         unsigned RTot = Hit.CS_RY + Hit.CS_RN + Hit.CS_RA;\r
 312 \r
 313         double PctL = Pct(Hit.CS_LY, LTot);\r
 314         double PctR = Pct(Hit.CS_RY, RTot);\r
 315 \r
 316         fprintf(f,\r
 317           "Ids.  QA %.1f%%, QB %.1f%%, AB %.1f%%, QModel %.1f%%, Div. %+.1f%%\n",\r
 318           Hit.PctIdQA,\r
 319           Hit.PctIdQB,\r
 320           Hit.PctIdAB,\r
 321           Hit.PctIdQM,\r
 322           Div);\r
 323 \r
 324         fprintf(f,\r
 325           "Diffs Left %u: N %u, A %u, Y %u (%.1f%%); Right %u: N %u, A %u, Y %u (%.1f%%), Score %.4f\n",\r
 326           LTot, Hit.CS_LN, Hit.CS_LA, Hit.CS_LY, PctL,\r
 327           RTot, Hit.CS_RN, Hit.CS_RA, Hit.CS_RY, PctR,\r
 328           Hit.Score);\r
 329         }\r