← Back to branch summary

~ubuntu-branches/debian/sid/libvcflib/sid

« back to all changes in this revision

Viewing changes to src/vcfremap.cpp

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2016-09-16 15:52:29 UTC
  • Revision ID: package-import@ubuntu.com-20160916155229-24mxrntfylvsshsg
Tags: upstream-1.0.0~rc1+dfsg
Import upstream version 1.0.0~rc1+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/vcfremap.cpp
 
1
#include "Variant.h"
 
2
#include "BedReader.h"
 
3
#include "IntervalTree.h"
 
4
#include <getopt.h>
 
5
#include "Fasta.h"
 
6
#include <algorithm>
 
7
#include <list>
 
8
#include <set>
 
9
 
 
10
using namespace std;
 
11
using namespace vcflib;
 
12
 
 
13
 
 
14
void printSummary(char** argv) {
 
15
    cerr << "usage: " << argv[0] << " [options] [<vcf file>]" << endl
 
16
         << endl
 
17
         << "options:" << endl 
 
18
         << "    -w, --ref-window-size N      align using this many bases flanking each side of the reference allele" << endl
 
19
         << "    -s, --alt-window-size N      align using this many flanking bases from the reference around each alternate allele" << endl
 
20
         << "    -r, --reference FILE         FASTA reference file, required with -i and -u" << endl
 
21
         << "    -m, --match-score N          match score for SW algorithm" << endl
 
22
         << "    -x, --mismatch-score N       mismatch score for SW algorithm" << endl
 
23
         << "    -o, --gap-open-penalty N     gap open penalty for SW algorithm" << endl
 
24
         << "    -e, --gap-extend-penalty N   gap extension penalty for SW algorithm" << endl
 
25
         << "    -z, --entropy-gap-open       use entropy scaling for the gap open penalty" << endl
 
26
         << "    -R, --repeat-gap-extend N    penalize non-repeat-unit gaps in repeat sequence" << endl
 
27
         << "    -a, --adjust-vcf TAG         supply a new cigar as TAG in the output VCF" << endl
 
28
         << endl
 
29
         << "For each alternate allele, attempt to realign against the reference with lowered gap open penalty." << endl
 
30
         << "If realignment is possible, adjust the cigar and reference/alternate alleles." << endl;
 
31
    exit(0);
 
32
}
 
33
 
 
34
int main(int argc, char** argv) {
 
35
 
 
36
    string vcfFileName;
 
37
    string fastaFileName;
 
38
    int windowsize = 100;
 
39
    bool includePreviousBaseForIndels = false;
 
40
    bool useMNPs = true;
 
41
    int altwindowsize = 50;
 
42
 
 
43
    // constants for SmithWaterman algorithm
 
44
    float matchScore = 10.0f;
 
45
    float mismatchScore = -9.0f;
 
46
    float gapOpenPenalty = 15.0f;
 
47
    float gapExtendPenalty = 6.66f;
 
48
 
 
49
    bool useEntropy = false;
 
50
    bool useRepeatGapExtendPenalty = false;
 
51
    float repeatGapExtendPenalty = 1;
 
52
 
 
53
    bool adjustVcf = false;
 
54
    string adjustedTag = "remappedCIGAR";
 
55
 
 
56
    if (argc == 1)
 
57
        printSummary(argv);
 
58
 
 
59
    int c;
 
60
    while (true) {
 
61
        static struct option long_options[] =
 
62
            {
 
63
                /* These options set a flag. */
 
64
                //{"verbose", no_argument,       &verbose_flag, 1},
 
65
                {"help", no_argument, 0, 'h'},
 
66
                {"ref-window-size", required_argument, 0, 'w'},
 
67
                {"reference", required_argument, 0, 'r'},
 
68
                {"match-score", required_argument, 0, 'm'},
 
69
                {"mismatch-score", required_argument, 0, 'x'},
 
70
                {"gap-open-penalty", required_argument, 0, 'o'},
 
71
                {"gap-extend-penalty", required_argument, 0, 'e'},
 
72
                {"alt-window-size", required_argument, 0, 's'},
 
73
                {"entropy-gap-open", no_argument, 0, 'z'},
 
74
                {"repeat-gap-extend", no_argument, 0, 'R'},
 
75
                {"adjust-vcf", required_argument, 0, 'a'},
 
76
                {0, 0, 0, 0}
 
77
            };
 
78
        /* getopt_long stores the option index here. */
 
79
        int option_index = 0;
 
80
 
 
81
        c = getopt_long (argc, argv, "hza:w:r:m:x:o:e:s:R:",
 
82
                         long_options, &option_index);
 
83
 
 
84
        if (c == -1)
 
85
            break;
 
86
 
 
87
        switch (c) {
 
88
 
 
89
            case 'w':
 
90
            windowsize = atoi(optarg);
 
91
            break;
 
92
 
 
93
            case 'a':
 
94
                adjustVcf = true;
 
95
            adjustedTag = optarg;
 
96
            break;
 
97
 
 
98
            case 'r':
 
99
            fastaFileName = string(optarg);
 
100
            break;
 
101
 
 
102
        case 'h':
 
103
            printSummary(argv);
 
104
            break;
 
105
 
 
106
            case 'm':
 
107
            matchScore = atof(optarg);
 
108
                break;
 
109
 
 
110
            case 'x':
 
111
            mismatchScore = atof(optarg);
 
112
                break;
 
113
 
 
114
            case 'o':
 
115
            gapOpenPenalty = atof(optarg);
 
116
                break;
 
117
 
 
118
            case 'e':
 
119
            gapExtendPenalty = atof(optarg);
 
120
                break;
 
121
 
 
122
            case 's':
 
123
            altwindowsize = atoi(optarg);
 
124
            break;
 
125
 
 
126
            case 'z':
 
127
            useEntropy = true;
 
128
            break;
 
129
 
 
130
            case 'R':
 
131
            useRepeatGapExtendPenalty = true;
 
132
            repeatGapExtendPenalty = atof(optarg);
 
133
            break;
 
134
 
 
135
        case '?':
 
136
            printSummary(argv);
 
137
            exit(1);
 
138
            break;
 
139
 
 
140
        default:
 
141
            abort ();
 
142
        }
 
143
    }
 
144
 
 
145
    VariantCallFile variantFile;
 
146
    string inputFilename;
 
147
    if (optind == argc - 1) {
 
148
        inputFilename = argv[optind];
 
149
        variantFile.open(inputFilename);
 
150
    } else {
 
151
        variantFile.open(std::cin);
 
152
    }
 
153
 
 
154
    if (!variantFile.is_open()) {
 
155
        cerr << "could not open VCF file" << endl;
 
156
        exit(1);
 
157
    }
 
158
 
 
159
    FastaReference freference;
 
160
    if (fastaFileName.empty()) {
 
161
        cerr << "a reference is required" << endl;
 
162
        exit(1);
 
163
    } else {
 
164
        freference.open(fastaFileName);
 
165
    }
 
166
    
 
167
    if (adjustVcf) {
 
168
        vector<string> commandline;
 
169
        for (int i = 0; i < argc; ++i)
 
170
            commandline.push_back(argv[i]);
 
171
        variantFile.addHeaderLine("##INFO=<ID=" + adjustedTag + ",Number=A,Type=String,Description=\"CIGAR when remapped using"+ join(commandline, " ") +"\">");
 
172
    }
 
173
 
 
174
    cout << variantFile.header << endl;
 
175
 
 
176
    Variant var(variantFile);
 
177
    while (variantFile.getNextVariant(var)) {
 
178
        //if (!adjustVcf) {
 
179
            cout << endl;
 
180
            cout << var << endl;
 
181
            //}
 
182
        map<string, vector<VariantAllele> > variantAlleles;
 
183
        vector<vector<pair<int, char> > > cigars;
 
184
        vector<int> positionDiffs;
 
185
        for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
 
186
            //if (!adjustVcf) cout << endl;
 
187
            cout << endl;
 
188
 
 
189
            // try to remap locally
 
190
 
 
191
            string reference = freference.getSubSequence(var.sequenceName, var.position - 1 - windowsize, windowsize * 2 + var.ref.size());
 
192
            
 
193
            // passed to sw align
 
194
            unsigned int referencePos;
 
195
            string cigar;
 
196
 
 
197
            string& alternate = *a;
 
198
 
 
199
            vector<VariantAllele>& variants = variantAlleles[alternate];
 
200
 
 
201
            string alternateQuery = reference.substr(windowsize - altwindowsize, altwindowsize) + alternate + reference.substr(reference.size() - windowsize, altwindowsize);
 
202
 
 
203
            //cout << "REF:\t" << reference << endl;
 
204
            //cout << "ALT:\t" << string(windowsize - altwindowsize, ' ') << alternateQuery << endl;
 
205
            
 
206
            CSmithWatermanGotoh sw(matchScore, mismatchScore, gapOpenPenalty, gapExtendPenalty);
 
207
            if (useEntropy) sw.EnableEntropyGapPenalty(1);
 
208
            if (useRepeatGapExtendPenalty) sw.EnableRepeatGapExtensionPenalty(repeatGapExtendPenalty);
 
209
            sw.Align(referencePos, cigar, reference, alternateQuery);
 
210
 
 
211
            int altpos = 0;
 
212
            int refpos = 0;
 
213
            int len;
 
214
            string slen;
 
215
            vector<pair<int, char> > cigarData;
 
216
 
 
217
            string ref = reference.substr(referencePos);
 
218
            positionDiffs.push_back(referencePos); // TODO this... is borked
 
219
 
 
220
            stringstream refss;
 
221
            stringstream altss;
 
222
 
 
223
            if (!adjustVcf) cout << cigar << endl;
 
224
            cout << cigar << endl;
 
225
            for (string::iterator c = cigar.begin(); c != cigar.end(); ++c) {
 
226
                switch (*c) {
 
227
                case 'I':
 
228
                    len = atoi(slen.c_str());
 
229
                    slen.clear();
 
230
                    if (altpos < altwindowsize) {
 
231
                        cigarData.push_back(make_pair(len, 'M'));
 
232
                    } else {
 
233
                        cigarData.push_back(make_pair(len, *c));
 
234
                    }
 
235
                    altss << alternateQuery.substr(altpos, len);
 
236
                    refss << string(len, '-');
 
237
                    altpos += len;
 
238
                    break;
 
239
                case 'D':
 
240
                    len = atoi(slen.c_str());
 
241
                    slen.clear();
 
242
                    if (altpos < altwindowsize) {
 
243
                    } else {
 
244
                        cigarData.push_back(make_pair(len, *c));
 
245
                    }
 
246
                    refss << ref.substr(refpos, len);
 
247
                    altss << string(len, '-');
 
248
                    refpos += len;
 
249
                    break;
 
250
                case 'M':
 
251
                    len = atoi(slen.c_str());
 
252
                    slen.clear();
 
253
                    {
 
254
                        for (int i = 0; i < len; ++i) {
 
255
                            if (ref.at(refpos + i) == alternateQuery.at(altpos + i)) {
 
256
                                if (!cigarData.empty() && cigarData.back().second == 'M') {
 
257
                                    cigarData.back().first++;
 
258
                                } else {
 
259
                                    cigarData.push_back(make_pair(1, 'M'));
 
260
                                }
 
261
                            } else {
 
262
                                if (!cigarData.empty() && cigarData.back().second == 'X') {
 
263
                                    cigarData.back().first++;
 
264
                                } else {
 
265
                                    cigarData.push_back(make_pair(1, 'X'));
 
266
                                }
 
267
                            }
 
268
                        }
 
269
                    }
 
270
                    refss << ref.substr(refpos, len);
 
271
                    altss << alternateQuery.substr(altpos, len);
 
272
                    refpos += len;
 
273
                    altpos += len;
 
274
                    break;
 
275
                case 'S':
 
276
                    len = atoi(slen.c_str());
 
277
                    slen.clear();
 
278
                    cigarData.push_back(make_pair(len, *c));
 
279
                    refss << ref.substr(refpos, len);
 
280
                    //altss << alternateQuery.substr(altpos, len); // TODO deal with soft clipping, weird behavior
 
281
                    refpos += len;
 
282
                    altpos += len;
 
283
                    break;
 
284
                default:
 
285
                    len = 0;
 
286
                    slen += *c;
 
287
                    break;
 
288
                }
 
289
            }
 
290
 
 
291
            if (!adjustVcf) {
 
292
                cout << "ref:\t" << refss.str() << endl;
 
293
                cout << "alt:\t" << altss.str() << endl;
 
294
            } else {
 
295
                cout << "ref:\t" << refss.str() << endl;
 
296
                cout << "alt:\t" << altss.str() << endl;
 
297
                cigars.push_back(cigarData);
 
298
            }
 
299
 
 
300
        }
 
301
 
 
302
        if (adjustVcf) {
 
303
            int substart = cigars.front().front().first;
 
304
            int subend = cigars.front().back().first;
 
305
 
 
306
            // find the min and max match
 
307
            for (vector<vector<pair<int, char> > >::iterator c = cigars.begin(); c != cigars.end(); ++c) {
 
308
                if (c->front().second == 'M' && c->front().first <= substart) {
 
309
                    substart = c->front().first;
 
310
                    if (c->size() > 1 && c->at(1).second != 'X') {
 
311
                        --substart;
 
312
                    }
 
313
                }
 
314
                if (c->back().second == 'M' && c->back().first <= subend) {
 
315
                    subend = c->back().first;
 
316
                }
 
317
            }
 
318
            
 
319
            // adjust the cigars and get the new reference length
 
320
            int reflen = 0;
 
321
            for (vector<vector<pair<int, char> > >::iterator c = cigars.begin(); c != cigars.end(); ++c) {
 
322
                c->front().first -= substart;
 
323
                c->back().first -= subend;
 
324
                int crf = cigarRefLen(*c);
 
325
                if (crf > reflen)
 
326
                    reflen = crf;
 
327
                var.info[adjustedTag].push_back(joinCigar(*c));
 
328
            }
 
329
 
 
330
            // find the lowest positional difference
 
331
            int pdiff = 0;
 
332
            for (vector<int>::iterator d = positionDiffs.begin(); d != positionDiffs.end(); ++d) {
 
333
                if (*d + altwindowsize < pdiff)
 
334
                    pdiff = *d + altwindowsize;
 
335
            }
 
336
 
 
337
            // adjust the reference string
 
338
            var.position += pdiff;
 
339
 
 
340
            // adjust the variant position
 
341
            var.ref = freference.getSubSequence(var.sequenceName, var.position - 1, reflen);
 
342
 
 
343
            cout << var << endl;
 
344
        }
 
345
    }
 
346
 
 
347
    return 0;
 
348
 
 
349
}
 
350