sfasta.cpp

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