uchime_src/sfasta.cpp

   1 #include "sfasta.h"\r
   2 #include "orf.h"\r
   3 #include "alpha.h"\r
   4 #include "timing.h"\r
   5 \r
   6 static inline bool isgap(byte c)\r
   7         {\r
   8         return c == '-' || c == '.';\r
   9         }\r
  10 \r
  11 const unsigned BufferSize = 16*1024*1024;\r
  12 \r
  13 static unsigned GetMaxPoly(const byte *Seq, unsigned L)\r
  14         {\r
  15         byte CurrChar = Seq[0];\r
  16         unsigned Start = 0;\r
  17         unsigned MaxLen = 1;\r
  18         for (unsigned i = 1; i < L; ++i)\r
  19                 {\r
  20                 char c = Seq[i];\r
  21                 if (c != CurrChar || i+1 == L)\r
  22                         {\r
  23                         unsigned Len = i - Start;\r
  24                         if (Len > MaxLen)\r
  25                                 MaxLen = Len;\r
  26                         CurrChar = c;\r
  27                         Start = i;\r
  28                         }\r
  29                 }\r
  30         return MaxLen;\r
  31         }\r
  32 \r
  33 SFasta::SFasta()\r
  34         {\r
  35         m_FileName = "";\r
  36         m_File = 0;\r
  37         m_Buffer = 0;\r
  38         m_BufferSize = 0;\r
  39         m_BufferOffset = 0;\r
  40         m_BufferBytes = 0;\r
  41         m_FilePos = 0;\r
  42         m_FileSize = 0;\r
  43         m_Label = 0;\r
  44         m_SeqLength = 0;\r
  45         m_TooShortCount = 0;\r
  46         m_TooLongCount = 0;\r
  47         m_ShortestLength = 0;\r
  48         m_LongestLength = 0;\r
  49         m_IsNucleo = false;\r
  50         m_IsNucleoSet = false;\r
  51         }\r
  52 \r
  53 SFasta::~SFasta()\r
  54         {\r
  55         Clear();\r
  56         }\r
  57 \r
  58 void SFasta::Clear()\r
  59         {\r
  60         MYFREE(m_Buffer, m_BufferSize, SFasta);\r
  61         if (m_File != 0)\r
  62                 CloseStdioFile(m_File);\r
  63 \r
  64         m_FileName = "";\r
  65         m_File = 0;\r
  66         m_Buffer = 0;\r
  67         m_BufferSize = 0;\r
  68         m_BufferOffset = 0;\r
  69         m_BufferBytes = 0;\r
  70         m_FilePos = 0;\r
  71         m_FileSize = 0;\r
  72         m_Label = 0;\r
  73         m_SeqLength = 0;\r
  74         m_SeqIndex = UINT_MAX;\r
  75         m_AllowGaps = false;\r
  76         m_IsNucleo = false;\r
  77         m_IsNucleoSet = false;\r
  78         m_TooShortCount = 0;\r
  79         m_TooLongCount = 0;\r
  80         m_ShortestLength = 0;\r
  81         m_LongestLength = 0;\r
  82         m_TooPolyCount = 0;\r
  83         }\r
  84 \r
  85 void SFasta::LogMe() const\r
  86         {\r
  87         Log("\n");\r
  88         Log("SFasta::LogMe()\n");\r
  89         Log("FileName=%s\n", m_FileName.c_str());\r
  90         Log("FileSize=%u\n", (unsigned) m_FileSize);\r
  91         Log("FilePos=%u\n", (unsigned) m_FilePos);\r
  92         Log("BufferSize=%u\n", m_BufferSize);\r
  93         Log("BufferPos=%u\n", m_BufferOffset);\r
  94         Log("BufferBytes=%u\n", m_BufferBytes);\r
  95         if (m_Label == 0)\r
  96                 Log("Label=NULL\n");\r
  97         else\r
  98                 Log("Label=%s\n", m_Label);\r
  99         Log("SeqLength=%u\n", m_SeqLength);\r
 100         }\r
 101 \r
 102 const byte *SFasta::GetNextSeq()\r
 103         {\r
 104         for (;;)\r
 105                 {\r
 106                 const byte *Seq = GetNextSeqLo();\r
 107                 if (Seq == 0)\r
 108                         {\r
 109                         if (m_TooShortCount > 0)\r
 110                                 Warning("%u short sequences (--minlen %u, shortest %u) discarded from %s",\r
 111                                   m_TooShortCount, opt_minlen, m_ShortestLength, m_FileName.c_str());\r
 112                         if (m_TooLongCount > 0)\r
 113                                 Warning("%u long sequences (--maxlen %u, longest %u) discarded from %s",\r
 114                                   m_TooLongCount, opt_maxlen, m_LongestLength, m_FileName.c_str());\r
 115                         if (m_TooPolyCount > 0)\r
 116                                 Warning("%u sequences with long homopolymers discarded (--maxpoly %u)",\r
 117                                   m_TooPolyCount, opt_maxpoly);\r
 118                         return 0;\r
 119                         }\r
 120                 if (m_SeqLength < opt_minlen)\r
 121                         {\r
 122                         ++m_TooShortCount;\r
 123                         if (m_ShortestLength == 0 || m_SeqLength < m_ShortestLength)\r
 124                                 m_ShortestLength = m_SeqLength;\r
 125                         continue;\r
 126                         }\r
 127                 if (m_SeqLength > opt_maxlen && opt_maxlen != 0)\r
 128                         {\r
 129                         if (m_LongestLength == 0 || m_SeqLength > m_LongestLength)\r
 130                                 m_LongestLength = m_SeqLength;\r
 131                         ++m_TooLongCount;\r
 132                         continue;\r
 133                         }\r
 134                 return Seq;\r
 135                 }\r
 136         }\r
 137 \r
 138 const byte *SFasta::GetNextSeqLo()\r
 139         {\r
 140 // End of cache?\r
 141         if (m_BufferOffset == m_BufferBytes)\r
 142                 {\r
 143         // End of file?\r
 144                 if (m_FilePos == m_FileSize)\r
 145                         return 0;\r
 146                 FillCache();\r
 147                 }\r
 148 \r
 149         StartTimer(SF_GetNextSeq);\r
 150         asserta(m_Buffer[m_BufferOffset] == '>');\r
 151         m_Label = (char *) (m_Buffer + m_BufferOffset + 1);\r
 152         \r
 153 //// Scan to end-of-line.\r
 154 //// Use dubious library function strchr() in the hope\r
 155 //// that it uses fast machine code.\r
 156 //      byte *ptr = (byte *) strchr(m_Label, '\n');\r
 157 //      asserta(ptr != 0);\r
 158 //      *ptr = 0;\r
 159 \r
 160         byte *ptr = 0;\r
 161         for (unsigned i = m_BufferOffset; i < m_BufferSize; ++i)\r
 162                 {\r
 163                 char c = m_Buffer[i];\r
 164                 if (c == '\n' || c == '\r')\r
 165                         {\r
 166                         ptr = m_Buffer + i;\r
 167                         break;\r
 168                         }\r
 169                 }\r
 170         asserta(ptr != 0);\r
 171 \r
 172         if (opt_trunclabels)\r
 173                 {\r
 174                 for (char *p = m_Label; *p; ++p)\r
 175                         if (isspace(*p))\r
 176                                 {\r
 177                                 *p = 0;\r
 178                                 break;\r
 179                                 }\r
 180                 }\r
 181         else\r
 182                 {\r
 183                 for (char *p = m_Label; *p; ++p)\r
 184                         {\r
 185                         if (*p == '\t')\r
 186                                 *p = ' ';\r
 187                         else if (*p == '\r' || *p == '\n')\r
 188                                 {\r
 189                                 *p = 0;\r
 190                                 char NextChar = *(p+1);\r
 191                                 if (NextChar == '\r' || NextChar == '\n')\r
 192                                         ++p;\r
 193                                 break;\r
 194                                 }\r
 195                         }\r
 196                 }\r
 197 \r
 198 // ptr points to end-of-line.\r
 199 // Move to start of sequence data.\r
 200         byte *Seq = ++ptr;\r
 201 \r
 202 // Delete white space in-place\r
 203         byte *To = ptr;\r
 204         m_BufferOffset = (unsigned) (ptr - m_Buffer);\r
 205         while (m_BufferOffset < m_BufferBytes)\r
 206                 {\r
 207                 byte c = m_Buffer[m_BufferOffset];\r
 208                 if (c == '>')\r
 209                         {\r
 210                         char prevc = '\n';\r
 211                         if (m_BufferOffset > 0)\r
 212                                 prevc = m_Buffer[m_BufferOffset-1];\r
 213                         if (prevc == '\n' || prevc == '\r')\r
 214                                 break;\r
 215                         }\r
 216                 ++m_BufferOffset;\r
 217                 if (isalpha(c) || (isgap(c) && m_AllowGaps))\r
 218                         *To++ = c;\r
 219                 else if (c == '\n' || c == '\r')\r
 220                         continue;\r
 221                 else\r
 222                         {\r
 223                         const char *Label = (m_Label == 0 ? "" : m_Label);\r
 224                         static bool WarningDone = false;\r
 225                         if (!WarningDone)\r
 226                                 {\r
 227                                 if (isgap(c))\r
 228                                         Warning("Ignoring gaps in FASTA file '%s'",\r
 229                                           m_FileName.c_str());\r
 230                                 else if (isprint(c))\r
 231                                         Warning("Invalid FASTA file '%s', non-letter '%c' in sequence >%s",\r
 232                                           m_FileName.c_str(), c, Label);\r
 233                                 else\r
 234                                         Warning("Invalid FASTA file '%s', non-printing byte (hex %02x) in sequence >%s",\r
 235                                           m_FileName.c_str(), c, Label);\r
 236                                 WarningDone = true;\r
 237                                 }\r
 238                         continue;\r
 239                         }\r
 240                 }\r
 241         m_SeqLength = unsigned(To - Seq);\r
 242 \r
 243         if (m_SeqIndex == UINT_MAX)\r
 244                 m_SeqIndex = 0;\r
 245         else\r
 246                 ++m_SeqIndex;\r
 247 \r
 248         EndTimer(SF_GetNextSeq);\r
 249         return Seq;\r
 250         }\r
 251 \r
 252 void SFasta::Open(const string &FileName)\r
 253         {\r
 254         Clear();\r
 255         m_FileName = FileName;\r
 256         m_File = OpenStdioFile(FileName);\r
 257         m_BufferSize = BufferSize;\r
 258         //m_Buffer = myalloc<byte>(m_BufferSize);\r
 259         m_Buffer = MYALLOC(byte, m_BufferSize, SFasta);\r
 260         m_FileSize = GetStdioFileSize(m_File);\r
 261         }\r
 262 \r
 263 void SFasta::Rewind()\r
 264         {\r
 265         m_BufferOffset = 0;\r
 266         m_BufferBytes = 0;\r
 267         m_FilePos = 0;\r
 268         }\r
 269 \r
 270 bool SFasta::SetIsNucleo()\r
 271         {\r
 272         if (m_FilePos != 0)\r
 273                 Die("SFasta::IsNucleo, not at BOF");\r
 274 \r
 275         unsigned LetterCount = 0;\r
 276         unsigned NucleoLetterCount = 0;\r
 277         for (;;)\r
 278                 {\r
 279                 const byte *Seq = GetNextSeq();\r
 280                 if (Seq == 0)\r
 281                         break;\r
 282                 unsigned L = GetSeqLength();\r
 283                 for (unsigned i = 0; i < L; ++i)\r
 284                         if (g_IsNucleoChar[Seq[i]])\r
 285                                 ++NucleoLetterCount;\r
 286                 LetterCount += L;\r
 287                 if (LetterCount > 256)\r
 288                         break;\r
 289                 }\r
 290         Rewind();\r
 291         if (LetterCount == 0)\r
 292                 {\r
 293                 m_IsNucleoSet = true;\r
 294                 m_IsNucleo = true;\r
 295                 return true;\r
 296                 }\r
 297 \r
 298 // Nucleo if more than 90% nucleo letters AGCTUN\r
 299         m_IsNucleo = double(NucleoLetterCount)/LetterCount > 0.9;\r
 300         m_IsNucleoSet = true;\r
 301         return m_IsNucleo;\r
 302         }\r
 303 \r
 304 void SFasta::FillCache()\r
 305         {\r
 306         StartTimer(SF_FillCache);\r
 307         asserta(m_FilePos < m_FileSize);\r
 308 \r
 309 // off_t may be larger type than unsigned, e.g. 64- vs. 32-bit.\r
 310         off_t otBytesToRead = m_FileSize - m_FilePos;\r
 311 \r
 312         bool FinalBuffer = true;\r
 313         if (otBytesToRead > (off_t) m_BufferSize)\r
 314                 {\r
 315                 FinalBuffer = false;\r
 316                 otBytesToRead = m_BufferSize;\r
 317                 }\r
 318 \r
 319         unsigned BytesToRead = unsigned(otBytesToRead);\r
 320         asserta(BytesToRead > 0);\r
 321         asserta(BytesToRead <= m_BufferSize);\r
 322 \r
 323         SetStdioFilePos(m_File, m_FilePos);\r
 324         ReadStdioFile(m_File, m_Buffer, BytesToRead);\r
 325         if (m_Buffer[0] != '>')\r
 326                 {\r
 327                 if (m_FilePos == 0)\r
 328                         Die("Input is not FASTA file");\r
 329                 else\r
 330                         Die("SFasta::FillCache() failed, expected '>'");\r
 331                 }\r
 332 \r
 333         m_BufferOffset = 0;\r
 334 \r
 335 // If last buffer in file, done\r
 336         if (FinalBuffer)\r
 337                 {\r
 338                 m_BufferBytes = BytesToRead;\r
 339                 m_FilePos += BytesToRead;\r
 340                 EndTimer(SF_FillCache);\r
 341                 return;\r
 342                 }\r
 343 \r
 344 // If not last buffer, truncate any partial sequence\r
 345 // at end of buffer. Search backwards to find last '>'.\r
 346         byte *ptr = m_Buffer + BytesToRead - 1;\r
 347         while (ptr > m_Buffer)\r
 348                 {\r
 349                 if (ptr[0] == '>' && (ptr[-1] == '\n' || ptr[-1] == '\r'))\r
 350                         break;\r
 351                 --ptr;\r
 352                 }\r
 353 \r
 354         if (ptr == m_Buffer)\r
 355                 {\r
 356                 LogMe();\r
 357                 if (*ptr != '>')\r
 358                         {\r
 359         // No '>' found.\r
 360         // This might techincally be legal FASTA if the entire\r
 361         // buffer is white space, but strange if not the last buffer\r
 362         // in the file, so quit anyway.\r
 363                         Die("Failed to find '>' (pos=%u, bytes=%u)",\r
 364                           (unsigned) m_FilePos, BytesToRead);\r
 365                         }\r
 366                 else\r
 367                         {\r
 368         // Entire buffer is one sequence which may be truncated.\r
 369                         Die("Sequence too long (pos=%u, bytes=%u)",\r
 370                           (unsigned) m_FilePos, BytesToRead);\r
 371                         }\r
 372                 }\r
 373 \r
 374         asserta(*ptr == '>');\r
 375 \r
 376         m_BufferBytes = unsigned(ptr - m_Buffer);\r
 377         m_FilePos += m_BufferBytes;\r
 378 \r
 379         EndTimer(SF_FillCache);\r
 380         }\r
 381 \r
 382 unsigned SFasta::GetPctDoneX10() const\r
 383         {\r
 384         if (m_FilePos == 0 || m_FileSize == 0)\r
 385                 return 0;\r
 386 \r
 387         assert(m_FilePos >= (off_t) m_BufferBytes);\r
 388         off_t BufferStart = m_FilePos - m_BufferBytes;\r
 389         off_t BufferPos = BufferStart + m_BufferOffset;\r
 390 \r
 391         unsigned iPctX10 = unsigned(10.0*double(BufferPos)*100.0/double(m_FileSize));\r
 392         if (iPctX10 == 0)\r
 393                 return 1;\r
 394         if (iPctX10 >= 999)\r
 395                 return 998;\r
 396         return iPctX10;\r
 397         }\r
 398 \r
 399 double SFasta::GetPctDone() const\r
 400         {\r
 401         if (m_FilePos == 0 || m_FileSize == 0)\r
 402                 return 0;\r
 403 \r
 404         assert(m_FilePos >= (off_t) m_BufferBytes);\r
 405         off_t BufferStart = m_FilePos - m_BufferBytes;\r
 406         off_t BufferPos = BufferStart + m_BufferOffset;\r
 407 \r
 408         return double(BufferPos)*100.0/double(m_FileSize);\r
 409         }\r
 410 \r
 411 bool SFasta::GetNextSD(SeqData &SD)\r
 412         {\r
 413         SD.Seq = GetNextSeq();\r
 414         if (SD.Seq == 0)\r
 415                 return false;\r
 416 \r
 417         SD.Label = GetLabel();\r
 418         SD.L = GetSeqLength();\r
 419         SD.Index = GetSeqIndex();\r
 420         SD.ORFParent = 0;\r
 421         SD.Nucleo = GetIsNucleo();\r
 422         SD.RevComp = false;\r
 423 \r
 424         return true;\r
 425         }\r
 426 \r
 427 #if     TEST\r
 428 void TestSFasta()\r
 429         {\r
 430         SFasta SF;\r
 431         SF.Open(opt_input);\r
 432 \r
 433         if (opt_verbose)\r
 434                 {\r
 435                 Log("  Index   Length  Label\n");\r
 436                 Log("-------  -------  -----\n");\r
 437                 }\r
 438 \r
 439         unsigned Index = 0;\r
 440         unsigned SeqCount = 0;\r
 441         double LetterCount = 0.0;\r
 442         ProgressStep(0, 1000, "Reading");\r
 443         for (;;)\r
 444                 {\r
 445                 const byte *Seq = SF.GetNextSeq();\r
 446                 if (Seq == 0)\r
 447                         break;\r
 448                 ProgressStep(SF.GetPctDoneX10(), 1000, "Reading");\r
 449                 const char *Label = SF.GetLabel();\r
 450                 unsigned L = SF.GetSeqLength();\r
 451                 ++SeqCount;\r
 452                 LetterCount += L;\r
 453 \r
 454                 if (opt_verbose)\r
 455                         {\r
 456                         Log(">%7u  %7u  '%s'\n", Index, L, Label);\r
 457                         Log("+%7.7s  %7.7s  \"%*.*s\"\n", "", "", L, L, Seq);\r
 458                         }\r
 459 \r
 460                 ++Index;\r
 461                 }\r
 462         ProgressStep(999, 1000, "Reading");\r
 463 \r
 464         Progress("%u seqs, %s letters\n", SeqCount, FloatToStr(LetterCount));\r
 465         Log("%u seqs, %s letters\n", SeqCount, FloatToStr(LetterCount));\r
 466         }\r
 467 #endif // TEST\r