← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/vcfstreamsort.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/vcfstreamsort.cpp
 
1
#include "Variant.h"
 
2
#include <algorithm>
 
3
#include <getopt.h>
 
4
#include "convert.h"
 
5
 
 
6
using namespace std;
 
7
using namespace vcflib;
 
8
 
 
9
bool listContains(list<string>& l, string& v) {
 
10
    for (list<string>::iterator i = l.begin(); i != l.end(); ++i) {
 
11
        if (*i == v) return true;
 
12
    }
 
13
    return false;
 
14
}
 
15
 
 
16
void printSummary(char** argv) {
 
17
    cerr << "usage: " << argv[0] << " [options] [vcf file]" << endl
 
18
         << endl
 
19
         << "Sorts the input (either stdin or file) using a streaming sort algorithm."
 
20
         << endl
 
21
         << "options:" << endl
 
22
         << endl
 
23
         << "    -h, --help             this dialog" << endl
 
24
         << "    -w, --window N         number of sites to sort (default 10000)" << endl
 
25
         << "    -a, --all              load all sites and then sort in memory" << endl;
 
26
}
 
27
 
 
28
int main(int argc, char** argv) {
 
29
 
 
30
    VariantCallFile variantFile;
 
31
    int sortSitesWindow = 10000;
 
32
    bool sortAll = false;
 
33
 
 
34
    int c;
 
35
 
 
36
    while (true) {
 
37
        static struct option long_options[] =
 
38
        {  
 
39
            /* These options set a flag. */
 
40
            //{"verbose", no_argument,       &verbose_flag, 1},
 
41
            {"help", no_argument, 0, 'h'},
 
42
            {"window", required_argument, 0, 'w'},
 
43
            {"all", required_argument, 0, 'a'},
 
44
            {0, 0, 0, 0}
 
45
        };
 
46
        /* getopt_long stores the option index here. */
 
47
        int option_index = 0;
 
48
 
 
49
        c = getopt_long (argc, argv, "haw:",
 
50
                         long_options, &option_index);
 
51
 
 
52
        if (c == -1)
 
53
            break;
 
54
 
 
55
        string field;
 
56
 
 
57
        switch (c)
 
58
        {
 
59
 
 
60
        case 'w':
 
61
            if (!convert(optarg, sortSitesWindow)) {
 
62
                cerr << "could not parse --window, -w" << endl;
 
63
                exit(1);
 
64
            }
 
65
            break;
 
66
                
 
67
        case 'a':
 
68
            sortAll = true;
 
69
            break;
 
70
 
 
71
        case 'h':
 
72
            printSummary(argv);
 
73
            exit(0);
 
74
            break;
 
75
            
 
76
        default:
 
77
            break;
 
78
        }
 
79
    }
 
80
 
 
81
    if (optind == argc - 1) {
 
82
        string inputFilename = argv[optind];
 
83
        variantFile.open(inputFilename);
 
84
    } else {
 
85
        variantFile.open(std::cin);
 
86
    }
 
87
 
 
88
    if (!variantFile.is_open()) {
 
89
        return 1;
 
90
    }
 
91
 
 
92
    cout << variantFile.header << endl;
 
93
 
 
94
    map<string, map<long int, map<string, vector<Variant> > > > records;
 
95
    long int back = 0;
 
96
    int numrecords = 0;
 
97
    list<string> sequenceNames;
 
98
 
 
99
    variantFile.parseSamples = false;
 
100
    Variant var(variantFile);
 
101
    while (variantFile.getNextVariant(var)) {
 
102
        //cerr << "at position " << var.sequenceName << ":" << var.position << endl;
 
103
        if (!listContains(sequenceNames, var.sequenceName)) {
 
104
            //cerr << "adding new sequence name " << var.sequenceName << endl;
 
105
            sequenceNames.push_back(var.sequenceName);
 
106
        }
 
107
        records[var.sequenceName][var.position][var.vrepr()].push_back(var);
 
108
        if (records[var.sequenceName][var.position].size() == 1) ++numrecords;
 
109
        if (!sortAll && numrecords > sortSitesWindow) {
 
110
            //cerr << "outputting a position" << endl;
 
111
            if (records[sequenceNames.front()].empty()) {
 
112
                //cerr << "end of reference sequence " << sequenceNames.front() << endl;
 
113
                sequenceNames.pop_front();
 
114
            }
 
115
            map<long int, map<string, vector<Variant> > >& frecords = records[sequenceNames.front()];
 
116
            map<string, vector<Variant> >& vars = frecords.begin()->second;
 
117
            for (map<string, vector<Variant> >::iterator v = vars.begin(); v != vars.end(); ++v) {
 
118
                for (vector<Variant>::iterator s = v->second.begin(); s != v->second.end(); ++s) {
 
119
                    cout << s->originalLine << endl;
 
120
                }
 
121
            }
 
122
            frecords.erase(frecords.begin());
 
123
            --numrecords;
 
124
        }
 
125
    }
 
126
    //cerr << "done processing input, cleaning up" << endl;
 
127
    for (list<string>::iterator s = sequenceNames.begin(); s != sequenceNames.end(); ++s) {
 
128
        map<long int, map<string, vector<Variant> > >& q = records[*s];
 
129
        for (map<long int, map<string, vector<Variant> > >::iterator r = q.begin(); r != q.end(); ++r) {
 
130
            for (map<string, vector<Variant> >::iterator v = r->second.begin(); v != r->second.end(); ++v) {
 
131
                for (vector<Variant>::iterator s = v->second.begin(); s != v->second.end(); ++s) {
 
132
                    cout << s->originalLine << endl;
 
133
                }
 
134
            }
 
135
            --numrecords;
 
136
        }
 
137
    }
 
138
    //cerr << numrecords << " remain" << endl;
 
139
 
 
140
    return 0;
 
141
 
 
142
}
 
143