uchime_src/seqdb.cpp

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