int main(int argc, char** argv) {

    if (argc != 2) {
        cerr << "usage: " << argv[0] << " <vcf file>" << endl
             << "outputs the het/hom ratio for each individual in the file" << endl;
        return 1;
    }

    string filename = argv[1];

    VariantCallFile variantFile;
    if (filename == "-") {
        variantFile.open(std::cin);
    } else {
        variantFile.open(filename);
    }
    if (!variantFile.is_open()) {
        cerr << "could not open " << filename << endl;
        return 1;
    }

    map<string, unsigned int> hetCounts;
    map<string, unsigned int> homCounts;
    for (vector<string>::iterator s = variantFile.sampleNames.begin(); s != variantFile.sampleNames.end(); ++s) {
        hetCounts[*s] = 0;
        homCounts[*s] = 0;
    }

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        //cout << var << endl;
        for (map<string, map<string, vector<string> > >::iterator s = var.samples.begin(); s != var.samples.end(); ++s) {
            string name = s->first;
            map<string, vector<string> >& sample = s->second;
            string& gt = sample["GT"].front();
            map<int, int> genotype = decomposeGenotype(gt);
            if (isHet(genotype)) {
                ++hetCounts[name];
            } else if (isHomNonRef(genotype)) {
                ++homCounts[name];
            }
        }
    }

    for (vector<string>::iterator s = variantFile.sampleNames.begin(); s != variantFile.sampleNames.end(); ++s) {
        cout << (s == variantFile.sampleNames.begin() ? "" : "\t") << *s;
    }
    cout << endl;
    for (vector<string>::iterator s = variantFile.sampleNames.begin(); s != variantFile.sampleNames.end(); ++s) {
        cout << (s == variantFile.sampleNames.begin() ? "" : "\t") << (double) hetCounts[*s] / (double) homCounts[*s];
    }
    cout << endl;

    return 0;

}

int main(int argc, char** argv) {

    VariantCallFile variantFile;

    if (argc > 1) {
        string filename = argv[1];
        variantFile.open(filename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    cout << variantFile.header << endl;

    string lastsn;
    long int lastpos;
    string lastref;
    vector<string> lastalt;

    variantFile.parseSamples = false;
    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        if (!lastsn.empty()
            && (lastsn == var.sequenceName
                && lastpos == var.position
                && lastref == var.ref
                && lastalt == var.alt)) {
            continue;
        } else {
            lastsn = var.sequenceName;
            lastpos = var.position;
            lastref = var.ref;
            lastalt = var.alt;
            cout << var.originalLine << endl;
        }
    }

    return 0;

}

int main(int argc, char** argv) {

    VariantCallFile variantFile;

    if (argc > 1) {
        string filename = argv[1];
        variantFile.open(filename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    //cout << variantFile.header << endl;

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        //cout << var << endl;
        double afref = 1;
        map<double, vector<string> > allelesByAf;
        vector<double> afd;
        vector<string>& afstr = var.info["AF"];
        for (vector<string>::iterator af = afstr.begin(); af != afstr.end(); ++af) {
            double r; convert(*af, r);
            afd.push_back(r);
        }
        vector<double>::iterator af = afd.begin();
        for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a, ++af) {
            afref -= *af;
            allelesByAf[*af].push_back(*a);
        }
        cout << var.ref;
        for (map<double, vector<string> >::reverse_iterator a = allelesByAf.rbegin(); a != allelesByAf.rend(); ++a) {
            cout << " -> " << join(a->second, ", ");
        }
        cout << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    VariantCallFile variantFile;

    if (argc > 1) {
        string filename = argv[1];
        variantFile.open(filename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    variantFile.addHeaderLine("##FORMAT=<ID=SN,Number=1,Type=String,Description=\"The name of the sample.\">");

    cout << variantFile.header << endl;

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        var.format.push_back("SN");
        for (map<string, map<string, vector<string> > >::iterator s = var.samples.begin();
             s != var.samples.end(); ++s) {
            s->second["SN"].clear();
            s->second["SN"].push_back(s->first);
        }
        cout << var << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    if (argc < 3) {
        cerr << "usage: " << argv[0] << " <vcf file> [FIELD1] [FIELD2] ..." << endl
             << "outputs each record in the vcf file, removing INFO fields not listed on the command line" << endl;
        return 1;
    }

    string filename = argv[1];

    set<string> fieldsToKeep;
    for (int i = 2; i < argc; ++i) {
        fieldsToKeep.insert(argv[i]);
    }

    VariantCallFile variantFile;
    if (filename == "-") {
        variantFile.open(std::cin);
    } else {
        variantFile.open(filename);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    Variant var(variantFile);

    vector<string> fieldsToErase;
    vector<string> infoIds = variantFile.infoIds();
    for (vector<string>::iterator i = infoIds.begin(); i != infoIds.end(); ++i) {
        if (!fieldsToKeep.count(*i)) {
            fieldsToErase.push_back(*i);
            variantFile.removeInfoHeaderLine(*i);
        }
    }

    // write the header
    cout << variantFile.header << endl;
 
    // print the records, filtering is done via the setting of varA's output sample names
    while (variantFile.getNextVariant(var)) {
        for (vector<string>::iterator f = fieldsToErase.begin(); f != fieldsToErase.end(); ++f) {
            var.info.erase(*f);
            var.infoFlags.erase(*f);
        }
        cout << var << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    VariantCallFile variantFile;

    if (argc > 1) {
        string filename = argv[1];
        variantFile.open(filename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    cout << variantFile.header;

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        map<string, vector<VariantAllele> > variants = var.parsedAlternates();
        cout << var << endl;
        for (map<string, vector<VariantAllele> >::iterator va = variants.begin(); va != variants.end(); ++va) {
            cout << " ( " << va->first << " :: ";
            vector<VariantAllele>& vars = va->second;
            vector<VariantAllele>::iterator g = vars.begin();
            for (; g != vars.end(); ++g) {
                cout << *g << "; ";
            }
            cout << " ) ";
        }
        cout << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    VariantCallFile variantFile;

    if (argc > 1) {
        string filename = argv[1];
        variantFile.open(filename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    variantFile.addHeaderLine("##INFO=<ID=length,Number=A,Type=Integer,Description=\"length(ALT) - length(REF) for each ALT\">");
    variantFile.addHeaderLine("##INFO=<ID=length.ref,Number=1,Type=Integer,Description=\"length(REF)\">");
    variantFile.addHeaderLine("##INFO=<ID=length.alt,Number=A,Type=Integer,Description=\"length(ALT) for each ALT\">");
    cout << variantFile.header << endl;

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        vector<string>& lengths = var.info["length"];
        lengths.clear();
        for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
            lengths.push_back(convert((int) a->size() - (int) var.ref.size()));
        }
        vector<string>& lengthsRef = var.info["length.ref"];
        lengthsRef.clear();
        lengthsRef.push_back(convert(var.ref.size()));
        vector<string>& lengthsAlt = var.info["length.alt"];
        lengthsAlt.clear();
        for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
            lengthsAlt.push_back(convert((int) a->size()));
        }
        cout << var << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    if (argc < 3) {
        cerr << "usage: " << argv[0] << " <vcf file> [SAMPLE1] [SAMPLE2] ..." << endl
             << "outputs each record in the vcf file, removing samples listed on the command line" << endl;
        return 1;
    }

    string filename = argv[1];

    vector<string> samplesToRemove;
    for (int i = 2; i < argc; ++i) {
        samplesToRemove.push_back(argv[i]);
    }

    VariantCallFile variantFile;
    if (filename == "-") {
        variantFile.open(std::cin);
    } else {
        variantFile.open(filename);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    Variant var(variantFile);

    vector<string> samplesToKeep = removeElems(samplesToRemove, variantFile.sampleNames);

    // update sample list in header
    variantFile.updateSamples(samplesToKeep);

    // and restrict the output sample names in the variant to those we are keeping
    var.setOutputSampleNames(samplesToKeep);
    
    // write the new header
    cout << variantFile.header << endl;
 
    // print the records, filtering is done via the setting of varA's output sample names
    while (variantFile.getNextVariant(var)) {
        cout << var << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    
    string ref_file = "";
    vector<string> insertion_files;
    int max_interval = -1;
    bool replace_sequences = true;

    int c = 0;
    while (true) {
        static struct option long_options[] =
            {
                {"insertions", no_argument, 0, 'i'},
                {"help", no_argument, 0, 'h'},
                {"reference", required_argument, 0, 'r'},
                {"no-replace-sequences", no_argument, 0, 's'},
                {0, 0, 0, 0}
            };
        int option_index = 0;

        c = getopt_long (argc, argv, "sr:i:h",
                         long_options, &option_index);
        if (c == -1)
            break;
        /* Detect the end of the options. */
        switch(c){
        case 's':
            replace_sequences = false;
            break;
        case 'r':
            ref_file = optarg;
            break;
        case 'i':
            insertion_files.push_back(optarg);
            break;
        case 'h':
        case '?':
            print_help(argv);
            exit(1);
        default:
            print_help(argv);
            abort();
        }
    }

    if (argc < 2){
        print_help(argv);
        exit(1);
    }

 

    VariantCallFile variantFile;
    string filename = argv[argc - 1];
    variantFile.open(filename);
    if (!variantFile.is_open()) {
        return 1;
    }

    vector<FastaReference*> insertions;
    if (!insertion_files.empty()){
        for (auto x : insertion_files){
            FastaReference* ins = new FastaReference();
            insertions.push_back(ins);
            ins->open(x);
        }
    }

    FastaReference ref;
    if(!ref_file.empty()){
        ref.open(ref_file);
    }


    cout << variantFile.header << endl;

    Variant var;
    while (variantFile.getNextVariant(var)) {
        bool valid = var.canonicalize_sv(ref, insertions, replace_sequences, max_interval);
        if (!valid){
            cerr << "Variant could not be normalized" << var << endl;
        }
        cout << var << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    int window = 150;
    VariantCallFile variantFile;
    string fastaFileName;

    int c;
    while (true) {
        static struct option long_options[] =
            {
                /* These options set a flag. */
                //{"verbose", no_argument,       &verbose_flag, 1},
                {"help", no_argument, 0, 'h'},
                {"reference", required_argument, 0, 'r'},
                {"window", required_argument, 0, 'w'},
                {0, 0, 0, 0}
            };
        /* getopt_long stores the option index here. */
        int option_index = 0;

        c = getopt_long (argc, argv, "hw:r:",
                         long_options, &option_index);

        if (c == -1)
            break;

        switch (c) {

	    case 'r':
            fastaFileName = optarg;
            break;

	    case 'w':
            window = atoi(optarg);
            break;

        case '?':
            printSummary(argv);
            exit(1);
            break;

        case 'h':
            printSummary(argv);
            break;

        default:
            abort ();
        }
    }

    if (optind < argc) {
        string filename = argv[optind];
        variantFile.open(filename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        cerr << "could not open VCF file" << endl;
        exit(1);
    }

    FastaReference fastaReference;
    if (fastaFileName.empty()) {
        cerr << "a reference is required" << endl;
        exit(1);
    } else {
        fastaReference.open(fastaFileName);
    }

    /*
    variantFile.addHeaderLine("##INFO=<ID=TYPE,Number=A,Type=String,Description=\"The type of allele, either snp, mnp, ins, del, or complex.\">");
    variantFile.addHeaderLine("##INFO=<ID=LEN,Number=A,Type=Integer,Description=\"allele length\">");
    if (!parseFlag.empty()) {
        variantFile.addHeaderLine("##INFO=<ID="+parseFlag+",Number=0,Type=Flag,Description=\"The allele was parsed using vcfallelicprimitives.\">");
    }
    */
    cout << variantFile.header << endl;

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {

        // if there is no indel, there is nothing to realign
        bool hasIndel = false;
        for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
            if (a->size() != var.ref.size()) {
                hasIndel = true;
                break;
            }
        }
        if (!hasIndel) {
            cout << var << endl;
            continue;
        }

        vector<AltAlignment> alignments;
        string ref;

        // determine window size to prevent mismapping with SW algorithm
        int currentWindow = window;
//.........这里部分代码省略.........

int main(int argc, char** argv) {

    if (argc != 2) {
        cerr << "usage: " << argv[0] << " <annotation-tag> <vcf file> <vcf file>" << endl
             << "adds a tag (BasesToNextVariant) to each variant record which indicates" << endl
	     << "the distance to the nearest variant" << endl;
        return 1;
    }

    string filename = argv[1];

    VariantCallFile variantFile;
    if (filename == "-") {
        variantFile.open(std::cin);
    } else {
        variantFile.open(filename);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    Variant varA(variantFile);
    Variant varB(variantFile);
    Variant varC(variantFile);

    vector<Variant*> vars;
    vars.push_back(&varA);
    vars.push_back(&varB);
    vars.push_back(&varC);
    
    for (vector<Variant*>::iterator v = vars.begin(); v != vars.end(); ++v) {
	variantFile.getNextVariant(**v);
    }

    string tag = "BasesToClosestVariant";
    string line = "##INFO=<ID=" + tag + ",Number=1,Type=Integer,Description=\"" \
	+ "Number of bases to the closest variant in the file.\">";
    variantFile.addHeaderLine(line);

    cout << variantFile.header << endl;

    // get the first distances
    if (vars.at(0)->sequenceName == vars.at(1)->sequenceName) {
	vars.at(0)->info[tag].push_back(convert(vars.at(1)->position - vars.at(0)->position));
    }

    while (variantFile.getNextVariant(*vars.back())) {

	if (vars.at(1)->sequenceName == vars.at(0)->sequenceName &&
	    vars.at(1)->sequenceName == vars.at(2)->sequenceName) {
	    vars.at(1)->info[tag].push_back(convert(min(vars.at(1)->position - vars.at(0)->position,
							vars.at(2)->position - vars.at(1)->position)));
	} else if (vars.at(1)->sequenceName == vars.at(0)->sequenceName) {
	    vars.at(1)->info[tag].push_back(convert(vars.at(1)->position - vars.at(0)->position));
	} else if (vars.at(2)->sequenceName == vars.at(1)->sequenceName) {
	    vars.at(1)->info[tag].push_back(convert(vars.at(2)->position - vars.at(1)->position));
	} else {
	    // don't add the tag
	}
	cout << *vars.front() << endl;
	// rotate
	Variant* v = vars.at(0);
	vars.at(0) = vars.at(1);
	vars.at(1) = vars.at(2);
	vars.at(2) = v;

    }

    // assign the last distances
    
    if (vars.at(0)->sequenceName == vars.at(1)->sequenceName) {
	vars.at(0)->info[tag].push_back(convert(vars.at(1)->position - vars.at(0)->position));
	cout << *vars.at(0) << endl;
	
	vars.at(1)->info[tag].push_back(convert(vars.at(1)->position - vars.at(0)->position));
	cout << *vars.at(1) << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    VariantCallFile variantFile;

    stringstream headerss;
    headerss << "##fileformat=VCFv4.0" << endl
             << "##source=vcfrandom" << endl
             << "##reference=/d2/data/references/build_37/human_reference_v37.fa" << endl
             << "##phasing=none" << endl
             << "##INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of samples with data\">" << endl
             << "##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Total read depth at the locus\">" << endl
             << "##INFO=<ID=AC,Number=1,Type=Integer,Description=\"Total number of alternate alleles in called genotypes\">" << endl
             << "##INFO=<ID=AN,Number=1,Type=Integer,Description=\"Total number of alleles in called genotypes\">" << endl
             << "##INFO=<ID=AF,Number=1,Type=Float,Description=\"Estimated allele frequency in the range (0,1]\">" << endl
             << "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">" << endl
             << "##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality, the Phred-scaled marginal (or unconditional) probability of the called genotype\">" << endl
             << "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Read Depth\">" << endl
             << "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tbill";

    string header = headerss.str();
    variantFile.openForOutput(header);

    cout << variantFile.header << endl;

    srand(time(NULL));

    vector<string> atgc;
    atgc.push_back("A");
    atgc.push_back("T");
    atgc.push_back("G");
    atgc.push_back("C");

    for (int i = 1; i < 10; ++i) {
        Variant var(variantFile);
        var.sequenceName = "one";
        var.id = ".";
        var.filter = ".";
        var.ref = atgc.at(rand() % 4);
        var.quality = 100;
        stringstream s;
        s << rand() % 100;
        var.info["DP"].push_back(s.str());
        var.format.push_back("GT");
        var.format.push_back("DP");
        var.position = i;
        for (vector<string>::iterator s = var.sampleNames.begin(); s != var.sampleNames.end(); ++s) {
            string& name = *s;
            var.alt.clear();
            var.alt.push_back(atgc.at(rand() % 4));
            var.alt.push_back(atgc.at(rand() % 4));
            var.samples[name]["GT"].push_back("0/1");
            stringstream dp;
            dp << floor(rand() % 100);
            var.samples[name]["DP"].push_back(dp.str());
        }
        cout << var << endl;
    }

    return 0;

}

int main(int argc, char** argv) {

    int c;
    string fastaRef;
    bool keepFailures = false;
    bool excludeFailures = false;

    if (argc == 1)
        printSummary(argv);

    while (true) {
        static struct option long_options[] =
            {
                /* These options set a flag. */
                //{"verbose", no_argument,       &verbose_flag, 1},
                {"help", no_argument, 0, 'h'},
                {"fasta-reference",  required_argument, 0, 'f'},
                {"exclude-failures",  no_argument, 0, 'x'},
                {"keep-failures",  no_argument, 0, 'k'},
                //{"length",  no_argument, &printLength, true},
                {0, 0, 0, 0}
            };
        /* getopt_long stores the option index here. */
        int option_index = 0;

        c = getopt_long (argc, argv, "hxkf:",
                         long_options, &option_index);

        /* Detect the end of the options. */
        if (c == -1)
            break;
 
        switch (c)
        {
        case 0:
            /* If this option set a flag, do nothing else now. */
            if (long_options[option_index].flag != 0)
                break;
            printf ("option %s", long_options[option_index].name);
            if (optarg)
                printf (" with arg %s", optarg);
            printf ("\n");
            break;

        case 'f':
            fastaRef = optarg;
            break;

        case 'x':
            excludeFailures = true;
            break;

        case 'k':
            keepFailures = true;
            break;
 
        case 'h':
            printSummary(argv);
            exit(0);
            break;

        case '?':
            /* getopt_long already printed an error message. */
            printSummary(argv);
            exit(1);
            break;
 
        default:
            abort ();
        }
    }

    if (fastaRef.empty()) {
        cerr << "a FASTA reference sequence must be specified" << endl;
        exit(1);
    }

    FastaReference ref;
    ref.open(fastaRef);

    VariantCallFile variantFile;
    string inputFilename;
    if (optind == argc - 1) {
        inputFilename = argv[optind];
        variantFile.open(inputFilename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    if (keepFailures || excludeFailures) {
        cout << variantFile.header << endl;
    }

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        int refstart = var.position - 1; // convert to 0-based
//.........这里部分代码省略.........

int main(int argc, char** argv) {

    string bedFileName;
    string annotationInfoKey;
    string defaultAnnotationValue;

    if (argc == 1)
        printSummary(argv);

    int c;
    while (true) {
        static struct option long_options[] =
        {
            /* These options set a flag. */
            //{"verbose", no_argument,       &verbose_flag, 1},
            {"help", no_argument, 0, 'h'},
            {"bed",  required_argument, 0, 'b'},
            {"key",  required_argument, 0, 'k'},
            {"default",  required_argument, 0, 'd'},
            {0, 0, 0, 0}
        };
        /* getopt_long stores the option index here. */
        int option_index = 0;

        c = getopt_long (argc, argv, "hb:k:d:",
                         long_options, &option_index);

        if (c == -1)
            break;

        switch (c) {
            case 'b':
                bedFileName = string(optarg);
                break;

            case 'k':
                annotationInfoKey = string(optarg);
                break;

            case 'd':
                defaultAnnotationValue = string(optarg);
                break;

            case 'h':
                printSummary(argv);
                break;

            case '?':
                printSummary(argv);
                exit(1);
                break;

            default:
                abort ();
        }
    }

    if (bedFileName.empty()) {
        cerr << "a BED file is required when intersecting" << endl;
        exit(1);
    }

    BedReader bed(bedFileName);

    VariantCallFile variantFile;
    string inputFilename;
    if (optind == argc - 1) {
        inputFilename = argv[optind];
        variantFile.open(inputFilename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        cout << "could not open VCF file" << endl;
        return 1;
    }

    string line = "##INFO=<ID=" + annotationInfoKey + ",Number=1,Type=String,Description=\"Annotation from "
        + bedFileName + " delimited by ':'\">";
    variantFile.addHeaderLine(line);

    cout << variantFile.header << endl;

    Variant var(variantFile);
    while (variantFile.getNextVariant(var)) {
        BedTarget record(var.sequenceName, var.position, var.position + var.ref.size() - 1, "");
        vector<BedTarget*> overlaps = bed.targetsOverlapping(record);
        vector<string> annotations;
        if (!overlaps.empty()) {
            for (vector<BedTarget*>::iterator t = overlaps.begin(); t != overlaps.end(); ++t) {
                annotations.push_back((*t)->desc);
            }
            var.info[annotationInfoKey].push_back(join(annotations, ":"));
        } else if (!defaultAnnotationValue.empty()) {
            var.info[annotationInfoKey].push_back(defaultAnnotationValue);
        }
        cout << var << endl;
    }

//.........这里部分代码省略.........

int main(int argc, char** argv) {

  // set the random seed for MCMC

  srand((unsigned)time(NULL));

  // the filename

  string filename = "NA";

  // using vcflib; thanks to Erik Garrison 
  
  VariantCallFile variantFile ;

  // zero based index for the target and background indivudals 
  
  map<int, int> it, ib;
  
  // deltaaf is the difference of allele frequency we bother to look at 

  string deltaaf ;
  double daf  = -1;

    const struct option longopts[] = 
      {
	{"version"   , 0, 0, 'v'},
	{"help"      , 0, 0, 'h'},
        {"file"      , 1, 0, 'f'},
	{"target"    , 1, 0, 't'},
	{"background", 1, 0, 'b'},
	{"deltaaf"   , 1, 0, 'd'},
	{0,0,0,0}
      };

    int index;
    int iarg = 0;

    while(iarg != -1)
      {
	iarg = getopt_long(argc, argv, "d:t:b:f:hv", longopts, &index);
	
	switch (iarg)
	  {
	  case 0:
	    break;
	  case 'h':
	    cerr << endl;

	    cerr << "INFO: help: " << endl << endl;

	    cerr << "     bFst is a Bayesian approach to Fst.  Importantly bFst account for genotype uncertainty in the model using genotype likelihoods."       << endl;
	    cerr << "     For a more detailed description see: Holsinger et al. Molecular Ecology Vol 11, issue 7 2002.  The likelihood function has been "	     << endl;
	    cerr << "     modified to use genotype likelihoods provided by variant callers. There are five free parameters estimated in the model: each "	     << endl;
	    cerr << "     subpopulation's allele frequency and Fis (fixation index, within each subpopulation), a free parameter for the total population\'s "  << endl;
	    cerr << "     allele frequency, and Fst. "                                                                                      << endl             << endl;
	
	      cerr << "Output : 11 columns :                          " << endl; 
	      cerr << "     1.  Seqid                                     " << endl;
	      cerr << "     2.  Position				     " << endl;
	      cerr << "     3.  Observed allele frequency in target.	     " << endl;
	      cerr << "     4.  Estimated allele frequency in target.     " << endl;
	      cerr << "     5.  Observed allele frequency in background.  " << endl;
	      cerr << "     6.  Estimated allele frequency in background. " << endl;
	      cerr << "     7.  Observed allele frequency combined. 	     " << endl;
	      cerr << "     8.  Estimated allele frequency in combined.   " << endl;
	      cerr << "     9.  ML estimate of Fst (mean)		     " << endl;
	      cerr << "     10. Lower bound of the 95% credible interval  " << endl;
	      cerr << "     11. Upper bound of the 95% credible interval  " << endl << endl;
											 

	    cerr << "INFO: usage:  bFst --target 0,1,2,3,4,5,6,7 --background 11,12,13,16,17,19,22 --file my.vcf --deltaaf 0.1" << endl;
	    cerr << endl;
	    cerr << "INFO: required: t,target     -- a zero bases comma separated list of target individuals corrisponding to VCF columns" << endl;
	    cerr << "INFO: required: b,background -- a zero bases comma separated list of background individuals corrisponding to VCF columns" << endl;
	    cerr << "INFO: required: f,file a     -- a proper formatted VCF file.  the FORMAT field MUST contain \"PL\"" << endl; 
	    cerr << "INFO: required: d,deltaaf    -- skip sites were the difference in allele frequency is less than deltaaf" << endl;
	    cerr << endl; 
	    printVersion();
	    cerr << endl << endl;
	    return 0;

	  case 'v':
	    printVersion();
	    return 0;

	  case 't':
	    loadIndices(ib, optarg);
	    cerr << "INFO: There are " << ib.size() << " individuals in the target" << endl;
	    break;

	  case 'b':
	    loadIndices(it, optarg);
	    cerr << "INFO: There are " << it.size() << " individuals in the background" << endl;
	    break;

	  case 'f':
	    cerr << "INFO: File: " << optarg  <<  endl;
	    filename = optarg;
	    break;

//.........这里部分代码省略.........

int main(int argc, char** argv) {

    int c;
    string fastaRef;
    int windowSize = 0;

    if (argc == 1)
        printSummary(argv);

    while (true) {
        static struct option long_options[] =
        {
            /* These options set a flag. */
            //{"verbose", no_argument,       &verbose_flag, 1},
            {"help", no_argument, 0, 'h'},
            {"fasta-reference",  required_argument, 0, 'f'},
            {"window-size", required_argument, 0, 'w'},
            //{"length",  no_argument, &printLength, true},
            {0, 0, 0, 0}
        };
        /* getopt_long stores the option index here. */
        int option_index = 0;

        c = getopt_long (argc, argv, "hf:w:",
                         long_options, &option_index);

      /* Detect the end of the options. */
          if (c == -1)
            break;
 
          switch (c)
            {
            case 0:
            /* If this option set a flag, do nothing else now. */
            if (long_options[option_index].flag != 0)
              break;
            printf ("option %s", long_options[option_index].name);
            if (optarg)
              printf (" with arg %s", optarg);
            printf ("\n");
            break;

          case 'f':
            fastaRef = optarg;
            break;

          case 'w':
            windowSize = atoi(optarg);
            break;
 
          case 'h':
            printSummary(argv);
            exit(0);
            break;

          case '?':
            /* getopt_long already printed an error message. */
            printSummary(argv);
            exit(1);
            break;
 
          default:
            abort ();
          }
      }

    if (windowSize == 0) {
        cerr << "a window size must be specified" << endl;
        exit(1);
    }
    if (fastaRef.empty()) {
        cerr << "a FASTA reference sequence must be specified" << endl;
        exit(1);
    }

    FastaReference ref;
    ref.open(fastaRef);

    VariantCallFile variantFile;
    string inputFilename;
    if (optind == argc - 1) {
        inputFilename = argv[optind];
        variantFile.open(inputFilename);
    } else {
        variantFile.open(std::cin);
    }

    if (!variantFile.is_open()) {
        return 1;
    }

    variantFile.addHeaderLine("##INFO=<ID=EntropyLeft,Number=1,Type=Float,Description=\"Entropy of left-flanking sequence of "+ convert(windowSize) +"bp\">");
    variantFile.addHeaderLine("##INFO=<ID=EntropyCenter,Number=1,Type=Float,Description=\"Entropy of centered sequence of "+ convert(windowSize) +"bp\">");
    variantFile.addHeaderLine("##INFO=<ID=EntropyRight,Number=1,Type=Float,Description=\"Entropy of right-flanking sequence of "+ convert(windowSize) +"bp\">");
    variantFile.addHeaderLine("##INFO=<ID=EntropyRef,Number=1,Type=Float,Description=\"Entropy of REF allele\">");
    variantFile.addHeaderLine("##INFO=<ID=EntropyAlt,Number=A,Type=Float,Description=\"Entropy of each ALT allele\">");

    cout << variantFile.header << endl;

    Variant var(variantFile);
//.........这里部分代码省略.........

int main(int argc, char** argv) {

    if (argc > 1 && (argv[1] == "-h" || argv[1] == "--help")) {
        cerr << "usage: " << argv[0] << " <vcf file>" << endl
             << "outputs a VCF stream where AC and NS have been generated for each record using sample genotypes" << endl;
        return 1;
    }

    VariantCallFile variantFile;
    if (argc == 1 || (argc == 2 && argv[1] == "-")) {
        variantFile.open(std::cin);
        if (!variantFile.is_open()) {
            cerr << "vcffixup: could not open stdin" << endl;
            return 1;
        }
    } else {
        string filename = argv[1];
        variantFile.open(filename);
        if (!variantFile.is_open()) {
            cerr << "vcffixup: could not open " << filename << endl;
            return 1;
        }
    }

    Variant var(variantFile);

    // remove header lines we're going to add
    variantFile.removeInfoHeaderLine("AC");
    variantFile.removeInfoHeaderLine("AF");
    variantFile.removeInfoHeaderLine("NS");
    variantFile.removeInfoHeaderLine("AN");

    // and add them back, so as not to duplicate them if they are already there
    variantFile.addHeaderLine("##INFO=<ID=AC,Number=A,Type=Integer,Description=\"Total number of alternate alleles in called genotypes\">");
    variantFile.addHeaderLine("##INFO=<ID=AF,Number=A,Type=Float,Description=\"Estimated allele frequency in the range (0,1]\">");
    variantFile.addHeaderLine("##INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of samples with data\">");
    variantFile.addHeaderLine("##INFO=<ID=AN,Number=1,Type=Integer,Description=\"Total number of alleles in called genotypes\">");

    // write the new header
    cout << variantFile.header << endl;
 
    // print the records, filtering is done via the setting of varA's output sample names
    while (variantFile.getNextVariant(var)) {
        stringstream ns;
        ns << var.samples.size();
        var.info["NS"].clear();
        var.info["NS"].push_back(ns.str());

        var.info["AC"].clear();
        var.info["AF"].clear();
        var.info["AN"].clear();

        int allelecount = countAlleles(var);
        stringstream an;
        an << allelecount;
        var.info["AN"].push_back(an.str());

        for (vector<string>::iterator a = var.alt.begin(); a != var.alt.end(); ++a) {
            string& allele = *a;
            int altcount = countAlts(var, var.getAltAlleleIndex(allele) + 1);
            stringstream ac;
            ac << altcount;
            var.info["AC"].push_back(ac.str());
            stringstream af;
            af << (double) altcount / (double) allelecount;
            var.info["AF"].push_back(af.str());
        }
        cout << var << endl;
    }

    return 0;

}

//.........这里部分代码省略.........
    

    // add an option for dumping

//    for(st