6 #include "hspfinder.h"
\r
9 #include "mothurout.h"
\r
11 bool SearchChime(Ultra &U, const SeqData &QSD, float QAb,
\r
12 const AlnParams &AP, const AlnHeuristics &AH, HSPFinder &HF,
\r
13 float MinFractId, ChimeHit2 &Hit);
\r
16 FILE *g_fUChimeAlns;
\r
17 const vector<float> *g_SortVecFloat;
\r
18 bool g_UchimeDeNovo = false;
\r
23 //printf("UCHIME %s by Robert C. Edgar\n", MY_VERSION);
\r
24 //printf("http://www.drive5.com/uchime\n");
\r
26 //printf("This software is donated to the public domain\n");
\r
36 Die("SetBLOSUM62 not implemented");
\r
39 void ReadSubstMx(const string &/*FileName*/, Mx<float> &/*Mxf*/)
\r
41 Die("ReadSubstMx not implemented");
\r
49 static bool CmpDescVecFloat(unsigned i, unsigned j)
\r
51 return (*g_SortVecFloat)[i] > (*g_SortVecFloat)[j];
\r
54 void Range(vector<unsigned> &v, unsigned N)
\r
58 for (unsigned i = 0; i < N; ++i)
\r
62 void SortDescending(const vector<float> &Values, vector<unsigned> &Order)
\r
65 const unsigned N = SIZE(Values);
\r
67 g_SortVecFloat = &Values;
\r
68 sort(Order.begin(), Order.end(), CmpDescVecFloat);
\r
72 float GetAbFromLabel(const string &Label)
\r
74 vector<string> Fields;
\r
75 Split(Label, Fields, '/');
\r
76 const unsigned N = SIZE(Fields);
\r
77 for (unsigned i = 0; i < N; ++i)
\r
79 const string &Field = Fields[i];
\r
80 if (Field.substr(0, 3) == "ab=")
\r
82 string a = Field.substr(3, string::npos);
\r
83 return (float) atof(a.c_str());
\r
87 Die("Missing abundance /ab=xx/ in label >%s", Label.c_str());
\r
91 int uchime_main(int argc, char *argv[])
\r
94 m = MothurOut::getInstance();
\r
96 MyCmdLine(argc, argv);
\r
106 printf("uchime v" MY_VERSION ".%s\n", SVN_VERSION);
\r
110 //printf("uchime v" MY_VERSION ".%s\n", SVN_VERSION);
\r
111 //printf("by Robert C. Edgar\n");
\r
112 //printf("http://drive5.com/uchime\n");
\r
113 //printf("This code is donated to the public domain.\n");
\r
118 float MinFractId = 0.95f;
\r
120 MinFractId = (float) opt_id;
\r
122 Log("%8.2f minh\n", opt_minh);
\r
123 Log("%8.2f xn\n", opt_xn);
\r
124 Log("%8.2f dn\n", opt_dn);
\r
125 Log("%8.2f xa\n", opt_xa);
\r
126 Log("%8.2f mindiv\n", opt_mindiv);
\r
127 Log("%8u maxp\n", opt_maxp);
\r
129 if (opt_input == "" && opt_uchime != "")
\r
130 opt_input = opt_uchime;
\r
132 if (opt_input == "")
\r
133 Die("Missing --input");
\r
135 g_UchimeDeNovo = (opt_db == "");
\r
137 if (opt_uchimeout != "")
\r
138 g_fUChime = CreateStdioFile(opt_uchimeout);
\r
140 if (opt_uchimealns != "")
\r
141 g_fUChimeAlns = CreateStdioFile(opt_uchimealns);
\r
146 Input.FromFasta(opt_input);
\r
147 if (!Input.IsNucleo())
\r
148 Die("Input contains amino acid sequences");
\r
150 const unsigned QuerySeqCount = Input.GetSeqCount();
\r
151 vector<unsigned> Order;
\r
152 for (unsigned i = 0; i < QuerySeqCount; ++i)
\r
153 Order.push_back(i);
\r
155 if (g_UchimeDeNovo)
\r
158 for (unsigned i = 0; i < QuerySeqCount; ++i)
\r
160 const char *Label = Input.GetLabel(i);
\r
161 float Ab = GetAbFromLabel(Label);
\r
164 SortDescending(Abs, Order);
\r
165 DB.m_IsNucleoSet = true;
\r
166 DB.m_IsNucleo = true;
\r
170 DB.FromFasta(opt_db);
\r
171 if (!DB.IsNucleo())
\r
172 Die("Database contains amino acid sequences");
\r
175 vector<ChimeHit2> Hits;
\r
176 unsigned HitCount = 0;
\r
177 for (unsigned i = 0; i < QuerySeqCount; ++i)
\r
180 if (m->control_pressed) { break; }
\r
182 unsigned QuerySeqIndex = Order[i];
\r
185 Input.GetSeqData(QuerySeqIndex, QSD);
\r
188 if (g_UchimeDeNovo)
\r
189 QAb = GetAbFromLabel(QSD.Label);
\r
192 AlnParams &AP = *(AlnParams *) 0;
\r
193 AlnHeuristics &AH = *(AlnHeuristics *) 0;
\r
194 HSPFinder &HF = *(HSPFinder *) 0;
\r
195 bool Found = SearchChime(DB, QSD, QAb, AP, AH, HF, MinFractId, Hit);
\r
200 if (g_UchimeDeNovo)
\r
201 DB.AddSeq(QSD.Label, QSD.Seq, QSD.L);
\r
204 WriteChimeHit(g_fUChime, Hit);
\r
206 ProgressStep(i, QuerySeqCount, "%u/%u chimeras found (%.1f%%)", HitCount, i, Pct(HitCount, i+1));
\r
211 Log("%s: %u/%u chimeras found (%.1f%%)\n",
\r
212 opt_input.c_str(), HitCount, QuerySeqCount, Pct(HitCount, QuerySeqCount));
\r
214 CloseStdioFile(g_fUChime);
\r
215 CloseStdioFile(g_fUChimeAlns);
\r