← Back to branch summary

~ubuntu-branches/ubuntu/vivid/idba/vivid

« back to all changes in this revision

Viewing changes to src/tools/validate_reads_blat.cpp

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2014-02-13 13:57:55 UTC
  • Revision ID: package-import@ubuntu.com-20140213135755-c97jpl9wl5yzfwfw
Tags: upstream-1.1.1
Import upstream version 1.1.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/validate_reads_blat.cpp
 
1
/**
 
2
 * @file validate_reads_blat.cpp
 
3
 * @brief Use blat alignment result to validate reads.
 
4
 * @author Yu Peng (ypeng@cs.hku.hk)
 
5
 * @version 1.0.3
 
6
 * @date 2011-09-06
 
7
 */
 
8
 
 
9
#include <algorithm>
 
10
#include <cstdio>
 
11
#include <cstring>
 
12
#include <iostream>
 
13
#include <map>
 
14
#include <set>
 
15
#include <string>
 
16
#include <vector>
 
17
 
 
18
#include "misc/blat_record.h"
 
19
#include "misc/options_description.h"
 
20
#include "misc/utils.h"
 
21
#include "sequence/sequence.h"
 
22
#include "sequence/sequence_io.h"
 
23
 
 
24
using namespace std;
 
25
 
 
26
const int MaxLine = (1 << 20);
 
27
 
 
28
char line[MaxLine];
 
29
char buf[MaxLine];
 
30
 
 
31
int main(int argc, char *argv[])
 
32
{
 
33
    int min_contig = 100;
 
34
    double similar = 0.95;
 
35
    double complete_rate = 0.8;
 
36
    bool is_local = false;
 
37
    
 
38
    OptionsDescription desc;
 
39
    desc.AddOption("min_contig", "", min_contig, "minimum contigs");
 
40
    desc.AddOption("similar", "", similar, "similarity");
 
41
    desc.AddOption("complete_rate", "", complete_rate, "completeness");
 
42
    desc.AddOption("is_local", "", is_local, "local align");
 
43
 
 
44
    try
 
45
    {
 
46
        desc.Parse(argc, argv);
 
47
    }
 
48
    catch (exception &e)
 
49
    {
 
50
        cerr << e.what() << endl;
 
51
        cerr << "validate_contigs_blat - validate contigs by blat." << endl;
 
52
        cerr << "Usage: validate_contigs_blat ref.fa contigs.fa." << endl;
 
53
        cerr << "Allowed Options: " << endl;
 
54
        cerr << desc << endl;
 
55
        exit(1);
 
56
    }
 
57
 
 
58
    deque<Sequence> refs;
 
59
    deque<string> ref_names;
 
60
    ReadSequence(argv[1], refs, ref_names);
 
61
 
 
62
    deque<Sequence> contigs;
 
63
    deque<string> contig_names;
 
64
    ReadSequence(argv[2], contigs, contig_names);
 
65
 
 
66
    vector<int> is_found(refs.size());
 
67
    vector<vector<double> > flags(refs.size());
 
68
    map<string, int> dict;
 
69
    for (unsigned i = 0; i < refs.size(); ++i)
 
70
    {
 
71
        flags[i].resize(refs[i].size(), false);
 
72
        size_t index = ref_names[i].find(' ');
 
73
        if (index != string::npos)
 
74
            ref_names[i].resize(index);
 
75
        dict[ref_names[i]] = i;
 
76
    }
 
77
 
 
78
    int num_gaps = 0;
 
79
    for (unsigned i = 0; i < contigs.size(); ++i)
 
80
    {
 
81
        size_t index = contig_names[i].find(' ');
 
82
        if (index != string::npos)
 
83
            contig_names[i].resize(index);
 
84
 
 
85
        bool is_new_gap = true;
 
86
        for (unsigned j = 0; j < contigs[i].size(); ++j)
 
87
        {
 
88
            if (contigs[i][j] == 4)
 
89
            {
 
90
                if (is_new_gap)
 
91
                {
 
92
                    is_new_gap = false;
 
93
                    ++num_gaps;
 
94
                }
 
95
            }
 
96
            else
 
97
                is_new_gap = true;
 
98
        }
 
99
    }
 
100
 
 
101
    string blat_file = string(argv[2]) + ".blat";
 
102
    FILE *fblat = OpenFile(blat_file, "rb");
 
103
 
 
104
    map<string, int> valid_contigs;
 
105
    deque<int> valid_lengths;
 
106
    int64_t num_mismatch = 0;
 
107
    while (fgets(line, MaxLine, fblat) != NULL)
 
108
    {
 
109
        BlatRecord record;
 
110
        record.Parse(line);
 
111
 
 
112
        deque<BlatRecord> records;
 
113
        records.push_back(record);
 
114
 
 
115
        while (fgets(line, MaxLine, fblat) != NULL)
 
116
        {
 
117
            record.Parse(line);
 
118
            if (record.query_name == records.back().query_name)
 
119
                records.push_back(record);
 
120
            else
 
121
            {
 
122
                fseek(fblat, -strlen(line), SEEK_CUR);
 
123
                break;
 
124
            }
 
125
        }
 
126
 
 
127
        int index = 0;
 
128
        for (unsigned i = 0; i < records.size(); ++i)
 
129
        {
 
130
            if (records[i].match_count > similar * records[i].query_length
 
131
                    && records[i].match_count > similar * abs(record.ref_to - record.ref_from))
 
132
                records[index++] = records[i];
 
133
        }
 
134
        records.resize(index);
 
135
 
 
136
        for (unsigned i = 0; i < records.size(); ++i)
 
137
        {
 
138
            record = records[i];
 
139
            int ref_id = dict[record.ref_name];
 
140
 
 
141
            //if (record.match_count > similar * record.query_length && record.query_length >= min_contig
 
142
            if ((record.match_count > similar * record.query_length 
 
143
                        || (is_local && record.match_count > similar * abs(record.query_to - record.query_from)))
 
144
            //if (record.match_count > similar * abs(record.query_to - record.query_from)
 
145
                    && abs(record.query_to - record.query_from) >= min_contig
 
146
                    && record.match_count > similar * abs(record.ref_to - record.ref_from)
 
147
               )
 
148
            {
 
149
                //if (record.match_count >= similar * record.ref_length)
 
150
                if (record.match_count >= complete_rate * record.ref_length)
 
151
                    is_found[ref_id] = true;
 
152
    //            else
 
153
    //                continue;
 
154
 
 
155
                int not_used = 0;
 
156
                for (unsigned i = 0; i < record.blocks.size(); ++i)
 
157
                {
 
158
                    BlatBlock block = record.blocks[i];
 
159
                    for (unsigned j = block.ref_from; j < block.ref_from + block.size; ++j)
 
160
                    {
 
161
                        if (flags[ref_id][j] == false)
 
162
                        {
 
163
                            //flags[ref_id][j] = true;
 
164
                            not_used++;
 
165
                        }
 
166
                        
 
167
                        flags[ref_id][j] += 1.0 / records.size();
 
168
                    }
 
169
                }
 
170
 
 
171
                if (valid_contigs.find(record.query_name) == valid_contigs.end())
 
172
                {
 
173
                    valid_contigs[record.query_name] = record.mismatch_count;
 
174
                    valid_lengths.push_back(record.query_to - record.query_from);
 
175
                }
 
176
                else
 
177
                {
 
178
                    valid_contigs[record.query_name] = min(record.mismatch_count, (int64_t)valid_contigs[record.query_name]);
 
179
 
 
180
                    if (not_used > similar * record.query_length)
 
181
                        valid_lengths.push_back(record.query_to - record.query_from);
 
182
                }
 
183
            }
 
184
        }
 
185
    }
 
186
 
 
187
    for (map<string, int>::iterator p = valid_contigs.begin(); p != valid_contigs.end(); ++p)
 
188
    {
 
189
        num_mismatch += p->second;
 
190
    }
 
191
 
 
192
    long long count = 0;
 
193
    long long total = 0;
 
194
    for (unsigned k = 0; k < flags.size(); ++k)
 
195
    {
 
196
        for (unsigned i = 0; i < flags[k].size(); ++i)
 
197
        {
 
198
            if (flags[k][i])
 
199
                ++count;
 
200
            ++total;
 
201
        }
 
202
    }
 
203
 
 
204
    //valid_lengths.push_back(60000);
 
205
    sort(valid_lengths.begin(), valid_lengths.end());
 
206
    reverse(valid_lengths.begin(), valid_lengths.end());
 
207
 
 
208
    long long n50 = 0;
 
209
    long long sum = 0;
 
210
    long long n80 = 0;
 
211
 
 
212
    for (unsigned i = 0; i < valid_lengths.size(); ++i)
 
213
    {
 
214
        sum += valid_lengths[i];
 
215
        if (sum >= 0.5 * total && n50 == 0)
 
216
            n50 = valid_lengths[i];
 
217
        if (sum >= 0.8 * total && n80 == 0)
 
218
            n80 = valid_lengths[i];
 
219
    }
 
220
    cout << "total " << total << " " << sum << endl;
 
221
 
 
222
    long long maximum = 0;
 
223
    long long mean = 0;
 
224
    if (valid_lengths.size() > 0)
 
225
    {
 
226
        maximum = valid_lengths[0];
 
227
        mean = sum / valid_lengths.size();
 
228
    }
 
229
 
 
230
    long long sum_wrong = 0;
 
231
    long long num_wrong = 0;
 
232
    long long corret_contigs = 0;
 
233
    long long sum_corret = 0;
 
234
    int last_id = 0;
 
235
    int last_error = 0;
 
236
    deque<int> contig_flags(contigs.size(), false);
 
237
    FastaWriter error_writer(FormatString("%s.error.fa", argv[2]));
 
238
    for (unsigned i = 0; i < contigs.size(); ++i)
 
239
    {
 
240
        if ((int)contigs[i].size() < min_contig)
 
241
            continue;
 
242
 
 
243
        if (valid_contigs.find(contig_names[i]) == valid_contigs.end())
 
244
        {
 
245
            ++num_wrong;
 
246
            sum_wrong += contigs[i].size();
 
247
            error_writer.Write(contigs[i], contig_names[i]);
 
248
        }
 
249
        else
 
250
        {
 
251
            last_id = i;
 
252
            last_error = sum_wrong;
 
253
 
 
254
            ++corret_contigs;
 
255
            sum_corret += contigs[i].size();
 
256
            contig_flags[i] = true;
 
257
            //correct_writer.Write(contigs[i], contig_names[i]);
 
258
        }
 
259
    }
 
260
 
 
261
    printf("last id %d %d total contigs %d gaps %d\n", last_id, last_error, (int)(num_wrong + corret_contigs), num_gaps);
 
262
    printf("contigs: %lld N50: %lld coverage: %.2f%% max: %lld mean: %lld total: %lld/%lld N80: %lld\n",
 
263
            (long long)valid_contigs.size(), n50, count * 100.0 / total, maximum, mean, count, total, n80);
 
264
    printf("substitution error: %.4f%% wrong contigs: %lld %lld correct: %lld %lld %s\n", 
 
265
            num_mismatch * 100.0 /sum, num_wrong, sum_wrong, corret_contigs, sum_corret, argv[2]);
 
266
 
 
267
    deque<int> lengths;
 
268
    for (unsigned i = 0; i < refs.size(); ++i)
 
269
    {
 
270
        int last = 0;
 
271
        for (unsigned j = 0; j < refs[i].size(); ++j)
 
272
        {
 
273
            if (flags[i][j] == 0)
 
274
            {
 
275
                if (flags[i][last])
 
276
                {
 
277
                    lengths.push_back(j - last);
 
278
                    last = j;
 
279
                }
 
280
            }
 
281
            else
 
282
            {
 
283
                if (flags[i][last] == 0)
 
284
                    last = j;
 
285
            }
 
286
        }
 
287
    }
 
288
    sort(lengths.begin(), lengths.end());
 
289
    reverse(lengths.begin(), lengths.end());
 
290
 
 
291
    deque<Sequence> gaps;
 
292
    deque<bool> is_no_long_gaps(refs.size());
 
293
    for (unsigned i = 0; i < refs.size(); ++i)
 
294
    {
 
295
        deque<int> tmp;
 
296
        Sequence gap;
 
297
 
 
298
        if (flags[i][0])
 
299
            tmp.push_back(0);
 
300
 
 
301
        for (unsigned j = 0; j < refs[i].size(); ++j)
 
302
        {
 
303
            if (flags[i][j] == false)
 
304
            {
 
305
                gap.Append(refs[i][j]);
 
306
            }
 
307
            else
 
308
            {
 
309
                if (gap.size() > 0)
 
310
                {
 
311
                    gaps.push_back(gap);
 
312
                    tmp.push_back(gap.size());
 
313
                }
 
314
                gap.resize(0);
 
315
            }
 
316
        }
 
317
 
 
318
        if (gap.size() > 0)
 
319
        {
 
320
            gaps.push_back(gap);
 
321
            tmp.push_back(gap.size());
 
322
        }
 
323
        else
 
324
            tmp.push_back(0);
 
325
 
 
326
        is_no_long_gaps[i] = true;
 
327
        for (unsigned j = 1; j+1 < tmp.size(); ++j)
 
328
        {
 
329
            if (tmp[j] > 50)
 
330
                is_no_long_gaps[i] = false;
 
331
        }
 
332
    }
 
333
    WriteSequence(FormatString("%s.gap.fa", argv[2]), gaps, "gap");
 
334
 
 
335
    FastaWriter ref_writer(argv[1] + string(".found.fa"));
 
336
    FILE *ffcound_list = OpenFile(argv[1] + string(".found.fa.list"), "wb");
 
337
    int found = 0;
 
338
    int covered = 0;
 
339
    int total_contigs = 0;
 
340
    for (unsigned i = 0; i < refs.size(); ++i)
 
341
    {
 
342
        int count = 0;
 
343
        double total_hit = 0;
 
344
        for (unsigned j = 0; j < flags[i].size(); ++j)
 
345
        {
 
346
            if (flags[i][j])
 
347
                ++count;
 
348
            total_hit += flags[i][j];
 
349
        }
 
350
 
 
351
        if (count > complete_rate * refs[i].size() 
 
352
                //&& is_no_long_gaps[i]
 
353
                //&& 1.0 * total_hit / count > 2
 
354
                )
 
355
        {
 
356
            ++covered;
 
357
            ref_writer.Write(refs[i], ref_names[i]);
 
358
            fprintf(ffcound_list, "%s %.4f\n", ref_names[i].c_str(), 1.0 * total_hit / count);
 
359
        }
 
360
        if (is_found[i])
 
361
            ++found;
 
362
    }
 
363
 
 
364
    int64_t total_bases = 0;
 
365
    for (unsigned i = 0; i < contigs.size(); ++i)
 
366
    {
 
367
        if ((int)contigs[i].size() >= min_contig)
 
368
        {
 
369
            ++total_contigs;
 
370
            total_bases += contigs[i].size();
 
371
        }
 
372
    }
 
373
 
 
374
    cout << corret_contigs << " " << total_contigs << " " << total_bases << endl;
 
375
    printf("cover ref: %d %d\n", covered, (int)refs.size());
 
376
    printf("found ref: %d %d\n", found, (int)refs.size());
 
377
    printf("precision: %.2f%% %d %d\n", 100.0 * corret_contigs / total_contigs, (int)corret_contigs, total_contigs);
 
378
 
 
379
    FastaWriter correct_writer(FormatString("%s.correct.fa", argv[2]));
 
380
    for (unsigned i = 0; i < contigs.size(); i += 2)
 
381
    {
 
382
        if (contig_flags[i] || contig_flags[i+1])
 
383
        {
 
384
            correct_writer.Write(contigs[i], contig_names[i]);
 
385
            correct_writer.Write(contigs[i+1], contig_names[i+1]);
 
386
        }
 
387
    }
 
388
    
 
389
    return 0;
 
390
}
 
391