Commit 36dcc828 authored by Dylan Aïssi's avatar Dylan Aïssi

New upstream version 2.4.1

parent 2d8ac15b
+-----------------------------+
| |
| Eagle v2.4 |
| December 13, 2017 |
| Eagle v2.4.1 |
| November 18, 2018 |
| Po-Ru Loh |
| |
+-----------------------------+
Copyright (C) 2015-2017 Harvard University.
Copyright (C) 2015-2018 Harvard University.
Distributed under the GNU GPLv3+ open source license.
Command line options:
......@@ -56,25 +56,25 @@ Setting --histFactor=1.00
BEGINNING STEP 1
Time for step 1: 0.282519
Time for step 1 MN^2: 0.0201082
Time for step 1: 0.347481
Time for step 1 MN^2: 0.0270685
Making hard calls (time: 0.0164621)
Making hard calls (time: 0.0181408)
BEGINNING STEP 2
BATCH 1 OF 1
Building hash tables
.................................................................. (time: 0.0826769)
.................................................................. (time: 0.065696)
Phasing samples 1-379
Time for phasing batch: 0.335781
Time for phasing batch: 0.309231
Making hard calls (time: 0.0172811)
Making hard calls (time: 0.0192151)
Time for step 2: 0.435747
Time for step 2 MN^2: 0.0391194
Time for step 2: 0.394149
Time for step 2 MN^2: 0.0408464
BEGINNING STEP 3 (PBWT ITERS)
......@@ -87,108 +87,108 @@ BEGINNING PBWT ITER 1
BATCH 1 OF 10
Phasing samples 1-37
Time for phasing batch: 1.13401
Time for phasing batch: 0.977095
BATCH 2 OF 10
Phasing samples 38-75
Time for phasing batch: 1.04853
Time for phasing batch: 0.898498
BATCH 3 OF 10
Phasing samples 76-113
Time for phasing batch: 1.18956
Time for phasing batch: 0.950794
BATCH 4 OF 10
Phasing samples 114-151
Time for phasing batch: 1.1353
Time for phasing batch: 0.937661
BATCH 5 OF 10
Phasing samples 152-189
Time for phasing batch: 1.01419
Time for phasing batch: 0.909736
BATCH 6 OF 10
Phasing samples 190-227
Time for phasing batch: 1.05172
Time for phasing batch: 0.924443
BATCH 7 OF 10
Phasing samples 228-265
Time for phasing batch: 0.99322
Time for phasing batch: 0.9072
BATCH 8 OF 10
Phasing samples 266-303
Time for phasing batch: 0.99847
Time for phasing batch: 0.926957
BATCH 9 OF 10
Phasing samples 304-341
Time for phasing batch: 0.998471
Time for phasing batch: 0.888218
BATCH 10 OF 10
Phasing samples 342-379
Time for phasing batch: 1.05052
Time for phasing batch: 0.90299
Time for PBWT iter 1: 10.6142
Time for PBWT iter 1: 9.22363
BEGINNING PBWT ITER 2
BATCH 1 OF 10
Phasing samples 1-37
Time for phasing batch: 1.76314
Time for phasing batch: 1.49903
BATCH 2 OF 10
Phasing samples 38-75
Time for phasing batch: 1.79145
Time for phasing batch: 1.4698
BATCH 3 OF 10
Phasing samples 76-113
Time for phasing batch: 1.8253
Time for phasing batch: 1.51891
BATCH 4 OF 10
Phasing samples 114-151
Time for phasing batch: 1.54563
Time for phasing batch: 1.45062
BATCH 5 OF 10
Phasing samples 152-189
Time for phasing batch: 1.54299
Time for phasing batch: 1.42069
BATCH 6 OF 10
Phasing samples 190-227
Time for phasing batch: 1.59936
Time for phasing batch: 1.4851
BATCH 7 OF 10
Phasing samples 228-265
Time for phasing batch: 1.55847
Time for phasing batch: 1.44716
BATCH 8 OF 10
Phasing samples 266-303
Time for phasing batch: 1.51538
Time for phasing batch: 1.49382
BATCH 9 OF 10
Phasing samples 304-341
Time for phasing batch: 1.58919
Time for phasing batch: 1.43647
BATCH 10 OF 10
Phasing samples 342-379
Time for phasing batch: 1.64308
Time for phasing batch: 1.46241
Time for PBWT iter 2: 16.3741
Time for PBWT iter 2: 14.6841
Writing .haps.gz and .sample output
Time for writing output: 0.281509
Total elapsed time for analysis = 28.1508 sec
Time for writing output: 0.221324
Total elapsed time for analysis = 24.9776 sec
+-----------------------------+
| |
| Eagle v2.4 |
| December 13, 2017 |
| Eagle v2.4.1 |
| November 18, 2018 |
| Po-Ru Loh |
| |
+-----------------------------+
Copyright (C) 2015-2017 Harvard University.
Copyright (C) 2015-2018 Harvard University.
Distributed under the GNU GPLv3+ open source license.
Command line options:
......@@ -18,7 +18,7 @@ Command line options:
--outPrefix=target.phased
Setting number of threads to 1
Warning: The index file is older than the data file: ref.bcf.csi
[W::hts_idx_load2] The index file is older than the data file: ref.bcf.csi
Reference samples: Nref = 169
Target samples: Ntarget = 8
......@@ -43,7 +43,7 @@ Average # SNPs per cM: 42
Number of <=(64-SNP, 1cM) segments: 9
Average # SNPs per segment: 47
Time for reading input: 6.18623 sec
Time for reading input: 5.09263 sec
Fraction of heterozygous genotypes: 0.178114
Typical span of default 100-het history length: 13.44 cM
......@@ -59,15 +59,15 @@ PHASING ITER 1 OF 1
Phasing target samples
................................................................................
Time for phasing iter 1: 0.144591
Time for phasing iter 1: 0.132139
Writing vcf.gz output to target.phased.vcf.gz
[W::vcf_parse] INFO 'AC' is not defined in the header, assuming Type=String
[W::vcf_parse] INFO 'AN' is not defined in the header, assuming Type=String
[W::vcf_parse] INFO 'DP' is not defined in the header, assuming Type=String
[W::vcf_parse] INFO 'AFR_AF' is not defined in the header, assuming Type=String
[W::vcf_parse] INFO 'EX_TARGET' is not defined in the header, assuming Type=String
Time for writing output: 0.0239
Total elapsed time for analysis = 6.35535 sec
Time for writing output: 0.0356269
Total elapsed time for analysis = 5.26099 sec
Mean phase confidence of each target individual:
ID PHASE_CONFIDENCE
......
+-----------------------------+
| |
| Eagle v2.4 |
| December 13, 2017 |
| Eagle v2.4.1 |
| November 18, 2018 |
| Po-Ru Loh |
| |
+-----------------------------+
Copyright (C) 2015-2017 Harvard University.
Copyright (C) 2015-2018 Harvard University.
Distributed under the GNU GPLv3+ open source license.
Command line options:
......@@ -39,25 +39,25 @@ Setting --histFactor=1.00
BEGINNING STEP 1
Time for step 1: 0.267283
Time for step 1 MN^2: 0.0164427
Time for step 1: 0.263286
Time for step 1 MN^2: 0.0227702
Making hard calls (time: 0.0152462)
Making hard calls (time: 0.018074)
BEGINNING STEP 2
BATCH 1 OF 1
Building hash tables
.................................................................. (time: 0.0786691)
.................................................................. (time: 0.0552511)
Phasing samples 1-379
Time for phasing batch: 0.331454
Time for phasing batch: 0.292423
Making hard calls (time: 0.0164559)
Making hard calls (time: 0.0192568)
Time for step 2: 0.426586
Time for step 2 MN^2: 0.0375612
Time for step 2: 0.366938
Time for step 2 MN^2: 0.0381578
BEGINNING STEP 3 (PBWT ITERS)
......@@ -70,108 +70,108 @@ BEGINNING PBWT ITER 1
BATCH 1 OF 10
Phasing samples 1-37
Time for phasing batch: 1.0778
Time for phasing batch: 1.02615
BATCH 2 OF 10
Phasing samples 38-75
Time for phasing batch: 1.0477
Time for phasing batch: 0.964455
BATCH 3 OF 10
Phasing samples 76-113
Time for phasing batch: 1.11441
Time for phasing batch: 0.988999
BATCH 4 OF 10
Phasing samples 114-151
Time for phasing batch: 1.06495
Time for phasing batch: 0.978426
BATCH 5 OF 10
Phasing samples 152-189
Time for phasing batch: 1.07034
Time for phasing batch: 0.935981
BATCH 6 OF 10
Phasing samples 190-227
Time for phasing batch: 1.13122
Time for phasing batch: 0.941438
BATCH 7 OF 10
Phasing samples 228-265
Time for phasing batch: 1.07535
Time for phasing batch: 0.935864
BATCH 8 OF 10
Phasing samples 266-303
Time for phasing batch: 1.16462
Time for phasing batch: 0.926938
BATCH 9 OF 10
Phasing samples 304-341
Time for phasing batch: 1.15337
Time for phasing batch: 0.927558
BATCH 10 OF 10
Phasing samples 342-379
Time for phasing batch: 1.17788
Time for phasing batch: 0.932059
Time for PBWT iter 1: 11.0777
Time for PBWT iter 1: 9.55792
BEGINNING PBWT ITER 2
BATCH 1 OF 10
Phasing samples 1-37
Time for phasing batch: 1.87655
Time for phasing batch: 1.56033
BATCH 2 OF 10
Phasing samples 38-75
Time for phasing batch: 1.6772
Time for phasing batch: 1.49036
BATCH 3 OF 10
Phasing samples 76-113
Time for phasing batch: 1.67862
Time for phasing batch: 1.51976
BATCH 4 OF 10
Phasing samples 114-151
Time for phasing batch: 1.57155
Time for phasing batch: 1.48635
BATCH 5 OF 10
Phasing samples 152-189
Time for phasing batch: 1.57551
Time for phasing batch: 1.43754
BATCH 6 OF 10
Phasing samples 190-227
Time for phasing batch: 1.64418
Time for phasing batch: 1.50166
BATCH 7 OF 10
Phasing samples 228-265
Time for phasing batch: 1.60561
Time for phasing batch: 1.48078
BATCH 8 OF 10
Phasing samples 266-303
Time for phasing batch: 1.63349
Time for phasing batch: 1.45079
BATCH 9 OF 10
Phasing samples 304-341
Time for phasing batch: 1.78741
Time for phasing batch: 1.47995
BATCH 10 OF 10
Phasing samples 342-379
Time for phasing batch: 1.80977
Time for phasing batch: 1.42872
Time for PBWT iter 2: 16.86
Time for PBWT iter 2: 14.8363
Writing vcf.gz output to phased.vcf.gz
Time for writing output: 0.204556
Total elapsed time for analysis = 34.7048 sec
Time for writing output: 0.176335
Total elapsed time for analysis = 30.3964 sec
No preview for this file type
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
This diff is collapsed.
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......@@ -191,7 +191,7 @@ namespace EAGLE {
int argc, char**argv) const;
void writeVcfNonRef(const std::string &vcfFile, const std::string &outFile, int inputChrom,
int chromX, double bpStart, double bpEnd, const std::string &writeMode,
bool keepMissingPloidyX, int argc, char **argv) const;
bool noImpMissing, bool keepMissingPloidyX, int argc, char **argv) const;
void makeHardCalls(uint64 n0start, uint64 n0end, uint seed);
void initRefIter(int refIter);
void buildHashTables(int iter, int batch, int seed);
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......@@ -83,8 +83,6 @@ void phaseWithRef(EagleParams &params, Timer &timer, double t0, int argc, char *
uint64 Nref = vcfData.getNref(), Ntarget = vcfData.getNtarget();
int iters = params.pbwtIters;
if (params.noImpMissing)
iters = 1;
if (iters == 0) {
if (Ntarget < Nref/2)
iters = 1;
......@@ -95,6 +93,7 @@ void phaseWithRef(EagleParams &params, Timer &timer, double t0, int argc, char *
cout << endl << "Auto-selecting number of phasing iterations: setting --pbwtIters to "
<< iters << endl << endl;
}
bool internalImpMissing = !params.noImpMissing || iters > 1; // internally impute missing
cout << endl << "BEGINNING PHASING" << endl;
......@@ -110,7 +109,7 @@ void phaseWithRef(EagleParams &params, Timer &timer, double t0, int argc, char *
HapBitsT *hapBitsTptr = NULL;
HapHedge *hapHedgePtr = NULL;
#ifdef OLD_IMP_MISSING
if (!params.noImpMissing) {
if (internalImpMissing) {
cout << "Making HapHedge" << endl;
hapBitsTptr = new HapBitsT(eagle.getHaploBitsT(), 2*eagle.getNlib(2+iter),
eagle.getMseg64(), eagle.getMaskSnps64j());
......@@ -138,10 +137,10 @@ void phaseWithRef(EagleParams &params, Timer &timer, double t0, int argc, char *
}
confs[i-Nref] = eagle.runPBWT
(i, nF1, nF2, params.Kpbwt/(iter<iters?2:1), params.expectIBDcM, params.histFactor,
iter==iters, iter>1, !params.noImpMissing, params.usePS, conPSall[i-Nref],
iter==iters, iter>1, internalImpMissing, params.usePS, conPSall[i-Nref],
params.chrom==params.chromX);
#ifdef OLD_IMP_MISSING
if (!params.noImpMissing) {
if (internalImpMissing) {
Timer tim;
eagle.imputeMissing(*hapHedgePtr, i);
timeImpMissing += tim.update_time();
......@@ -155,7 +154,7 @@ void phaseWithRef(EagleParams &params, Timer &timer, double t0, int argc, char *
}
}
#ifdef OLD_IMP_MISSING
if (!params.noImpMissing) {
if (internalImpMissing) {
delete hapHedgePtr;
delete hapBitsTptr;
}
......@@ -177,7 +176,7 @@ void phaseWithRef(EagleParams &params, Timer &timer, double t0, int argc, char *
cout << "Time for phasing iter " << iter << " MN^2: " << timeMN2 / params.numThreads
<< endl;
#ifdef OLD_IMP_MISSING
else if (!params.noImpMissing)
else if (internalImpMissing)
cout << "Time for phasing iter " << iter << " impMissing: "
<< timeImpMissing / params.numThreads << endl;
#endif
......@@ -216,7 +215,7 @@ int main(int argc, char *argv[]) {
cout << " +-----------------------------+" << endl;
cout << endl;
cout << "Copyright (C) 2015-2017 Harvard University." << endl;
cout << "Copyright (C) 2015-2018 Harvard University." << endl;
cout << "Distributed under the GNU GPLv3+ open source license." << endl << endl;
//cout << "Boost version: " << BOOST_LIB_VERSION << endl;
......@@ -503,8 +502,8 @@ int main(int argc, char *argv[]) {
string outFile = params.outPrefix + "." + params.vcfOutSuffix;
cout << "Writing " << params.vcfOutSuffix << " output to " << outFile << endl;
eagle.writeVcfNonRef(params.vcfFile, outFile, params.chrom, params.chromX, params.bpStart,
params.bpEnd, params.vcfWriteMode, params.keepMissingPloidyX, argc,
argv);
params.bpEnd, params.vcfWriteMode, params.noImpMissing,
params.keepMissingPloidyX, argc, argv);
}
else {
cout << "Writing .haps.gz and .sample output" << endl; timer.update_time();
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......@@ -237,10 +237,6 @@ namespace EAGLE {
if (vm.count("vcfRef") || vm.count("vcfTarget") || vm.count("vcf")) { // VCF mode
if (vm.count("vcf")) { // non-ref mode
if (noImpMissing) {
cerr << "ERROR: --noImpMissing is only supported in ref-mode" << endl;
return false;
}
if (bpFlanking != 0) {
cerr << "ERROR: --bpFlanking is only supported in ref-mode" << endl;
return false;
......@@ -255,13 +251,14 @@ namespace EAGLE {
cerr << "ERROR: --vcfTarget must be specified in reference-based phasing mode" << endl;
return false;
}
if (pbwtIters > 1 && noImpMissing) {
cerr << "ERROR: --pbwtIters cannot be greater than 1 if --noImpMissing is set" << endl;
return false;
}
if (vcfRef.substr(vcfRef.length()-4) != string(".bcf")) {
cerr << "WARNING: --vcfRef does not end in '.bcf'; BCF input is fastest" << endl;
}
if (chromStr=="0" && (bpStart>0 || bpEnd<1e9)) {
cerr << "ERROR: --chrom is required when specifying --bpStart or --bpEnd in ref mode"
<< endl;
return false;
}
}
// vcf input checks for both ref and non-ref mode
......@@ -291,7 +288,7 @@ namespace EAGLE {
return false;
}
}
else { // non-ref mode
else { // non-ref mode; plink input
if (famFile.empty()) {
cerr << "ERROR: fam file must be specified either using --fam or --bfile"
<< endl;
......@@ -308,7 +305,7 @@ namespace EAGLE {
return false;
}
if (noImpMissing) {
cerr << "ERROR: --noImpMissing is only supported in ref-mode" << endl;
cerr << "ERROR: --noImpMissing is only supported with vcf input" << endl;
return false;
}
if (bpFlanking != 0) {
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......
/*
This file is part of the Eagle haplotype phasing software package
developed by Po-Ru Loh. Copyright (C) 2015-2016 Harvard University.
developed by Po-Ru Loh. Copyright (C) 2015-2018 Harvard University.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
......@@ -24,8 +24,11 @@
#include <cstring>
#include <cmath>
#include <htslib/thread_pool.h>
#include <htslib/vcf.h>
#include "omp.h"
#include "Types.hpp"
#include "FileUtils.hpp"
#include "MemoryUtils.hpp"
......@@ -156,7 +159,7 @@ namespace EAGLE {
}
if (!ret.empty() &&
(snp.chrom < ret.back().chrom ||
(snp.chrom == ret.back().chrom && (snp.physpos <= ret.back().physpos ||
(snp.chrom == ret.back().chrom && (snp.physpos < ret.back().physpos ||
snp.genpos < ret.back().genpos)))) {
if (numOutOfOrder < 5) {
cerr << "WARNING: Out-of-order snp in bim file: " << bimFile << endl;
......@@ -310,7 +313,7 @@ namespace EAGLE {
preQCsnpInds.push_back(m); cMvec.push_back(100 * snpsPreQC[m].genpos);
}
seg64preQCsnpInds.push_back(preQCsnpInds); seg64cMvecs.push_back(cMvec);
Mseg64 = seg64preQCsnpInds.size();
cout << "Number of <=(64-SNP, " << cMmax << "cM) segments: " << Mseg64 << endl;
cout << "Average # SNPs per segment: " << M / Mseg64 << endl;
......@@ -369,13 +372,13 @@ namespace EAGLE {
af.cond[2][0] = log10safe(0);
af.cond[1][0] = log10safe(p1half / (p0 + p1half));
af.cond[0][0] = log10safe(p0 / (p0 + p1half));
af.cond[0][0] = log10safe(p0 / (p0 + p1half));
// same orientation as het-het => p(hap=1) = 1-p
af.cond[2][4] = log10safe((1-p) * p2 / (p1half + p2));
af.cond[1][4] = log10safe((1-p) * p1half / (p1half + p2) + p * p1half / (p0 + p1half));
af.cond[0][4] = log10safe(p * p0 / (p0 + p1half));
// opp orientation to het-het => p(hap=1) = p
af.cond[2][5] = log10safe(p * p2 / (p1half + p2));
af.cond[1][5] = log10safe(p * p1half / (p1half + p2) + (1-p) * p1half / (p0 + p1half));
......@@ -493,7 +496,7 @@ namespace EAGLE {
SGEMM_MACRO(&TRANSA_, &TRANSB_, &M_, &N_, &K_, &ALPHA_, A_, &LDA_, B_, &LDB_,
&BETA_, C_, &LDC_);
}
for (uint64 mPlus = 0; mPlus < mBlockCrop; mPlus++) {
uint64 m = mchr0 + mPlus;
if (!chrMaskSnps[m]) continue;
......@@ -539,7 +542,7 @@ namespace EAGLE {
return M/cMrange;
}
/**
/**
* reads indiv info from fam file, snp info from bim file
* allocates memory, reads genotypes, and does QC
*/
......@@ -548,7 +551,7 @@ namespace EAGLE {
const vector <string> &excludeFiles, const vector <string> &removeFiles,
double maxMissingPerSnp, double maxMissingPerIndiv, bool noMapCheck,
double cMmax) {
// indivsPreQC (without --remove indivs)
vector <bool> bedIndivRemoved = processIndivs(famFile, removeFiles);
// snpsPreQC (restricted to chrom:bpStart-bpEnd and without --exclude snps)
......@@ -645,7 +648,15 @@ namespace EAGLE {
cout << endl;
cout << "Total post-QC indivs: N = " << N << endl;
cout << "Total post-QC SNPs: M = " << M << endl;
if (N < 2) {
cerr << "ERROR: N>=2 individuals are required" << endl;
exit(1);
}
if (M == 0) {
cerr << "ERROR: At least 1 SNP is required" << endl;
exit(1);
}
cout << "MAF spectrum: " << endl;
const double mafBounds6[7] = {0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.500001};
vector <int> mafBinCounts(6);
......@@ -660,7 +671,7 @@ namespace EAGLE {
if (cMmax == 0) {
cMmax = std::min(1.0, std::max(N / 1e5, 0.25));
cout << "Auto-selecting --maxBlockLen: " << cMmax << " cM" << endl;
cout << "Auto-selecting --maxBlockLen: " << cMmax << " cM" << endl;
}
vector <bool> nullVec;
......@@ -669,7 +680,7 @@ namespace EAGLE {
ALIGNED_FREE(genosPreQC);
}
/**
/**
* reads genotypes from VCF/BCF file
* does not save indiv info (will be reread from VCF during output)
* only saves chrom, physpos, genpos in snp info (rest will be reread from VCF during output)
......@@ -678,12 +689,17 @@ namespace EAGLE {
void GenoData::initVcf(const string &vcfFile, const int inputChrom, const int chromX,
double bpStart, double bpEnd, const string &geneticMapFile,
bool noMapCheck, double cMmax) {
htsFile *fin = hts_open(vcfFile.c_str(), "r");
if (fin == NULL) {
cerr << "ERROR: Could not open " << vcfFile << " for reading" << endl;
exit(1);
}
htsThreadPool p = {hts_tpool_init(omp_get_max_threads()),
0};
hts_set_thread_pool(fin, &p);
bcf_hdr_t *hdr = bcf_hdr_read(fin);
bcf1_t *rec = bcf_init1();
int mgt = 0, *gt = NULL;
......@@ -710,7 +726,7 @@ namespace EAGLE {