seqdb.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 "seqdb.h"\r
   5 #include "alpha.h"\r
   6 #include "timing.h"\r
   7 #include "sfasta.h"\r
   8 #include "seq.h"\r
   9 \r
  10 void SeqToFasta(FILE *f, const char *Label, const byte *Seq, unsigned L)\r
  11         {\r
  12         const unsigned ROWLEN = 80;\r
  13         if (Label != 0)\r
  14                 fprintf(f, ">%s\n", Label);\r
  15         unsigned BlockCount = (L + ROWLEN - 1)/ROWLEN;\r
  16         for (unsigned BlockIndex = 0; BlockIndex < BlockCount; ++BlockIndex)\r
  17                 {\r
  18                 unsigned From = BlockIndex*ROWLEN;\r
  19                 unsigned To = From + ROWLEN;\r
  20                 if (To >= L)\r
  21                         To = L;\r
  22                 for (unsigned Pos = From; Pos < To; ++Pos)\r
  23                         fputc(Seq[Pos], f);\r
  24                 fputc('\n', f);\r
  25                 }\r
  26         }\r
  27 \r
  28 SeqDB::~SeqDB()\r
  29         {\r
  30         Clear();\r
  31         }\r
  32 \r
  33 SeqDB::SeqDB()\r
  34         {\r
  35         Clear(true);\r
  36         }\r
  37 \r
  38 void SeqDB::Clear(bool ctor)\r
  39         {\r
  40         if (!ctor)\r
  41                 {\r
  42                 for (unsigned i = 0; i < m_SeqCount; ++i)\r
  43                         {\r
  44                         unsigned n = strlen(m_Labels[i]);\r
  45                         MYFREE(m_Labels[i], n, SeqDB);\r
  46                         MYFREE(m_Seqs[i], m_SeqLengths[i], SeqDB);\r
  47                         }\r
  48                 MYFREE(m_Labels, m_Size, SeqDB);\r
  49                 MYFREE(m_Seqs, m_Size, SeqDB);\r
  50                 MYFREE(m_SeqLengths, m_Size, SeqDB);\r
  51                 }\r
  52 \r
  53         m_FileName.clear();\r
  54         m_SeqCount = 0;\r
  55         m_Size = 0;\r
  56 \r
  57         m_Labels = 0;\r
  58         m_Seqs = 0;\r
  59         m_SeqLengths = 0;\r
  60 \r
  61         m_Aligned = false;\r
  62         m_IsNucleo = false;\r
  63         m_IsNucleoSet = false;\r
  64         }\r
  65 \r
  66 void SeqDB::InitEmpty(bool Nucleo)\r
  67         {\r
  68         Clear();\r
  69         m_IsNucleo = Nucleo;\r
  70         m_IsNucleoSet = true;\r
  71         }\r
  72 \r
  73 void SeqDB::FromFasta(const string &FileName, bool AllowGaps)\r
  74         {\r
  75         Clear();\r
  76         m_FileName = FileName;\r
  77         SFasta SF;\r
  78 \r
  79         SF.Open(FileName);\r
  80         SF.m_AllowGaps = AllowGaps;\r
  81 \r
  82         ProgressStep(0, 1000, "Reading %s", FileName.c_str());\r
  83         for (;;)\r
  84                 {\r
  85                 unsigned QueryPctDoneX10 = SF.GetPctDoneX10();\r
  86                 ProgressStep(QueryPctDoneX10, 1000, "Reading %s", FileName.c_str());\r
  87                 const byte *Seq = SF.GetNextSeq();\r
  88                 if (Seq == 0)\r
  89                         break;\r
  90 \r
  91                 const char *Label = SF.GetLabel();\r
  92                 unsigned L = SF.GetSeqLength();\r
  93                 AddSeq(Label, Seq, L);\r
  94                 }\r
  95         ProgressStep(999, 1000, "Reading %s", FileName.c_str());\r
  96 \r
  97         SetIsNucleo();\r
  98 \r
  99         Progress("%s sequences\n", IntToStr(GetSeqCount()));\r
 100         }\r
 101 \r
 102 void SeqDB::ToFasta(const string &FileName) const\r
 103         {\r
 104         FILE *f = CreateStdioFile(FileName);\r
 105         for (unsigned SeqIndex = 0; SeqIndex < GetSeqCount(); ++SeqIndex)\r
 106                 ToFasta(f, SeqIndex);\r
 107         CloseStdioFile(f);\r
 108         }\r
 109 \r
 110 void SeqDB::SeqToFasta(FILE *f, unsigned SeqIndex, bool WithLabel) const\r
 111         {\r
 112         if (WithLabel)\r
 113                 fprintf(f, ">%s\n", GetLabel(SeqIndex));\r
 114 \r
 115         const unsigned ROWLEN = 80;\r
 116 \r
 117         unsigned L = GetSeqLength(SeqIndex);\r
 118         const byte *Seq = GetSeq(SeqIndex);\r
 119         unsigned BlockCount = (L + ROWLEN - 1)/ROWLEN;\r
 120         for (unsigned BlockIndex = 0; BlockIndex < BlockCount; ++BlockIndex)\r
 121                 {\r
 122                 unsigned From = BlockIndex*ROWLEN;\r
 123                 unsigned To = From + ROWLEN;\r
 124                 if (To >= L)\r
 125                         To = L;\r
 126                 for (unsigned Pos = From; Pos < To; ++Pos)\r
 127                         fputc(Seq[Pos], f);\r
 128                 fputc('\n', f);\r
 129                 }\r
 130         }\r
 131 \r
 132 void SeqDB::ToFasta(FILE *f, unsigned SeqIndex) const\r
 133         {\r
 134         asserta(SeqIndex < m_SeqCount);\r
 135         fprintf(f, ">%s\n", GetLabel(SeqIndex));\r
 136         SeqToFasta(f, SeqIndex);\r
 137         }\r
 138 \r
 139 unsigned SeqDB::GetMaxLabelLength() const\r
 140         {\r
 141         const unsigned SeqCount = GetSeqCount();\r
 142         unsigned MaxL = 0;\r
 143         for (unsigned Index = 0; Index < SeqCount; ++Index)\r
 144                 {\r
 145                 unsigned L = (unsigned) strlen(m_Labels[Index]);\r
 146                 if (L > MaxL)\r
 147                         MaxL = L;\r
 148                 }\r
 149         return MaxL;\r
 150         }\r
 151 \r
 152 unsigned SeqDB::GetMaxSeqLength() const\r
 153         {\r
 154         const unsigned SeqCount = GetSeqCount();\r
 155         unsigned MaxL = 0;\r
 156         for (unsigned Index = 0; Index < SeqCount; ++Index)\r
 157                 {\r
 158                 unsigned L = m_SeqLengths[Index];\r
 159                 if (L > MaxL)\r
 160                         MaxL = L;\r
 161                 }\r
 162         return MaxL;\r
 163         }\r
 164 \r
 165 void SeqDB::LogMe() const\r
 166         {\r
 167         Log("\n");\r
 168         const unsigned SeqCount = GetSeqCount();\r
 169         Log("SeqDB %u seqs, aligned=%c\n", SeqCount, tof(m_Aligned));\r
 170         if (SeqCount == 0)\r
 171                 return;\r
 172 \r
 173         Log("Index             Label  Length  Seq\n");\r
 174         Log("-----  ----------------  ------  ---\n");\r
 175         for (unsigned Index = 0; Index < SeqCount; ++Index)\r
 176                 {\r
 177                 Log("%5u", Index);\r
 178                 Log("  %16.16s", m_Labels[Index]);\r
 179                 unsigned L = m_SeqLengths[Index];\r
 180                 Log("  %6u", L);\r
 181                 Log("  %*.*s", L, L, m_Seqs[Index]);\r
 182                 Log("\n");\r
 183                 }\r
 184         }\r
 185 \r
 186 void SeqDB::GetSeqData(unsigned Id, SeqData &Buffer) const\r
 187         {\r
 188         asserta(Id < m_SeqCount);\r
 189         Buffer.Seq = m_Seqs[Id];\r
 190         Buffer.Label = m_Labels[Id];\r
 191         Buffer.L = m_SeqLengths[Id];\r
 192         Buffer.Index = Id;\r
 193         Buffer.ORFParent = 0;\r
 194         Buffer.RevComp = false;\r
 195         Buffer.Nucleo = IsNucleo();\r
 196         }\r
 197 \r
 198 void SeqDB::SetIsNucleo()\r
 199         {\r
 200         const unsigned SeqCount = GetSeqCount();\r
 201         unsigned N = 0;\r
 202         for (unsigned i = 0; i < 100; ++i)\r
 203                 {\r
 204                 unsigned SeqIndex = unsigned(rand()%SeqCount);\r
 205                 const byte *Seq = GetSeq(SeqIndex);\r
 206                 unsigned L = GetSeqLength(SeqIndex);\r
 207                 const unsigned Pos = unsigned(rand()%L);\r
 208                 byte c = Seq[Pos];\r
 209 \r
 210                 if (g_IsNucleoChar[c])\r
 211                         ++N;\r
 212                 }\r
 213         m_IsNucleo = (N > 80);\r
 214         m_IsNucleoSet = true;\r
 215         }\r
 216 \r
 217 unsigned SeqDB::GetTotalLength() const\r
 218         {\r
 219         const unsigned SeqCount = GetSeqCount();\r
 220         unsigned TotalLength = 0;\r
 221         for (unsigned Id = 0; Id < SeqCount; ++Id)\r
 222                 TotalLength += GetSeqLength(Id);\r
 223         return TotalLength;\r
 224         }\r
 225 \r
 226 unsigned SeqDB::AddSeq(const char *Label, const byte *Seq, unsigned L)\r
 227         {\r
 228         StartTimer(AddSeq);\r
 229         if (m_SeqCount >= m_Size)\r
 230                 {\r
 231                 unsigned NewSize = unsigned(m_Size*1.5) + 1024;\r
 232                 char **NewLabels = MYALLOC(char *, NewSize, SeqDB);\r
 233                 byte **NewSeqs = MYALLOC(byte *, NewSize, SeqDB);\r
 234                 unsigned *NewSeqLengths = MYALLOC(unsigned, NewSize, SeqDB);\r
 235 \r
 236                 for (unsigned i = 0; i < m_SeqCount; ++i)\r
 237                         {\r
 238                         NewLabels[i] = m_Labels[i];\r
 239                         NewSeqs[i] = m_Seqs[i];\r
 240                         NewSeqLengths[i] = m_SeqLengths[i];\r
 241                         }\r
 242 \r
 243                 MYFREE(m_Labels, m_SeqCount, SeqDB);\r
 244                 MYFREE(m_Seqs, m_SeqCount, SeqDB);\r
 245                 MYFREE(m_SeqLengths, m_SeqCount, SeqDB);\r
 246 \r
 247                 m_Labels = NewLabels;\r
 248                 m_Seqs = NewSeqs;\r
 249                 m_SeqLengths = NewSeqLengths;\r
 250                 m_Size = NewSize;\r
 251                 }\r
 252 \r
 253         unsigned Index = m_SeqCount++;\r
 254         m_Seqs[Index] = MYALLOC(byte, L, SeqDB);\r
 255         memcpy(m_Seqs[Index], Seq, L);\r
 256 \r
 257         unsigned n = strlen(Label) + 1;\r
 258         m_Labels[Index] = MYALLOC(char, n, SeqDB);\r
 259         memcpy(m_Labels[Index], Label, n);\r
 260 \r
 261         if (Index == 0)\r
 262                 m_Aligned = true;\r
 263         else\r
 264                 m_Aligned = (m_Aligned && L == m_SeqLengths[0]);\r
 265 \r
 266         m_SeqLengths[Index] = L;\r
 267 \r
 268         EndTimer(AddSeq);\r
 269         return Index;\r
 270         }\r
 271 \r
 272 unsigned SeqDB::GetIndex(const char *Label) const\r
 273         {\r
 274         for (unsigned i = 0; i < m_SeqCount; ++i)\r
 275                 if (strcmp(Label, m_Labels[i]) == 0)\r
 276                         return i;\r
 277         Die("SeqDB::GetIndex(%s), not found", Label);\r
 278         return UINT_MAX;\r
 279         }\r
 280 \r
 281 void SeqDB::MakeLabelToIndex(map<string, unsigned> &LabelToIndex)\r
 282         {\r
 283         LabelToIndex.clear();\r
 284         for (unsigned i = 0; i < m_SeqCount; ++i)\r
 285                 {\r
 286                 const string &Label = string(GetLabel(i));\r
 287                 if (LabelToIndex.find(Label) != LabelToIndex.end())\r
 288                         Die("Duplicate label: %s", Label.c_str());\r
 289                 LabelToIndex[Label] = i;\r
 290                 }\r
 291         }\r