本文整理汇总了C++中VariantCallFile类的典型用法代码示例。如果您正苦于以下问题:C++ VariantCallFile类的具体用法?C++ VariantCallFile怎么用?C++ VariantCallFile使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了VariantCallFile类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
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;
}
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:56,代码来源:vcfhethomratio.cpp
示例2: main
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;
}
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:43,代码来源:vcfuniq.cpp
示例3: main
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;
}
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:43,代码来源:vcfafpath.cpp
示例4: main
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;
}
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:33,代码来源:vcfgenosamplenames.cpp
示例5: main
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;
}
开发者ID:ShujiaHuang,项目名称:vcflib,代码行数:52,代码来源:vcfkeepinfo.cpp
示例6: main
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;
}
开发者ID:indapa,项目名称:vcflib,代码行数:36,代码来源:vcfparsealts.cpp
示例7: main
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;
}
开发者ID:ShujiaHuang,项目名称:vcflib,代码行数:41,代码来源:vcflength.cpp
示例8: main
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;
}
开发者ID:ShujiaHuang,项目名称:vcflib,代码行数:47,代码来源:vcfremovesamples.cpp
示例9: main
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;
}
开发者ID:vcflib,项目名称:vcflib,代码行数:88,代码来源:vcfnormalizesvs.cpp
示例10: main
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;
//.........这里部分代码省略.........
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:101,代码来源:vcfleftalign.cpp
示例11: main
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;
}
开发者ID:alimanfoo,项目名称:vcflib,代码行数:82,代码来源:vcfdistance.cpp
示例12: main
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;
}
开发者ID:ShujiaHuang,项目名称:vcflib,代码行数:61,代码来源:vcfrandom.cpp
示例13: main
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
//.........这里部分代码省略.........
开发者ID:hannespetur,项目名称:vcflib,代码行数:101,代码来源:vcfcheck.cpp
示例14: main
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;
}
//.........这里部分代码省略.........
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:101,代码来源:vcfannotate.cpp
示例15: main
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;
//.........这里部分代码省略.........
开发者ID:IvantheDugtrio,项目名称:vcflib,代码行数:101,代码来源:bFst.cpp
示例16: main
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);
//.........这里部分代码省略.........
开发者ID:MMesbahU,项目名称:vcflib,代码行数:101,代码来源:vcfentropy.cpp
示例17: main
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;
}
开发者ID:biocyberman,项目名称:vcflib,代码行数:73,代码来源:vcffixup.cpp
示例18: main
//.........这里部分代码省略.........
// add an option for dumping
// for(st
|
请发表评论