diff --git a/data/1000GP_Phase3.GBR.chr22.hap.gz b/data/1000GP_Phase3.GBR.chr22.hap.gz index 86435e8..bb457ae 100644 Binary files a/data/1000GP_Phase3.GBR.chr22.hap.gz and b/data/1000GP_Phase3.GBR.chr22.hap.gz differ diff --git a/data/chr22.map.gz b/data/chr22.map.gz index 5a5254b..f156bab 100644 Binary files a/data/chr22.map.gz and b/data/chr22.map.gz differ diff --git a/src/ihs.cpp b/src/ihs.cpp index eb14b31..1d5e62f 100644 --- a/src/ihs.cpp +++ b/src/ihs.cpp @@ -59,6 +59,7 @@ void calcIhsNoMpi( std::cout << "Calculations took " << std::chrono::duration(diff).count() << "ms" << std::endl; auto unStd = ihsfinder->unStdIHSByLine(); + auto freQs = ihsfinder->freqsByLine(); std::ofstream out2(outfile); out2 << "Index\tID\tFreq\tiHH_0\tiHH_1\tiHS\tStd iHS" << std::endl; @@ -66,7 +67,8 @@ void calcIhsNoMpi( for (const auto& it : res) { auto s = unStd[it.first]; - out2 << it.first << '\t' << hm.lineToId(it.first) << '\t' << s.freq << '\t' << s.iHH_0 << '\t' << s.iHH_1 << '\t' << s.iHS << "\t" << it.second << std::endl; + auto q = freQs[it.first]; + out2 << it.first << '\t' << hm.lineToId(it.first) << '\t' << q << '\t' << s.iHH_0 << '\t' << s.iHH_1 << '\t' << s.iHS << "\t" << it.second << std::endl; } std::cout << "# valid loci: " << res.size() << std::endl; std::cout << "# loci with MAF <= " << minMAF << ": " << ihsfinder->numOutsideMaf() << std::endl; diff --git a/src/ihs_mpi.cpp b/src/ihs_mpi.cpp index f7f0f69..f4e4ac2 100644 --- a/src/ihs_mpi.cpp +++ b/src/ihs_mpi.cpp @@ -152,6 +152,7 @@ void calcIhsMpi( IHSFinder::LineMap res = ihsfinder->normalize(); auto unStd = ihsfinder->unStdIHSByLine(); + auto freQs = ihsfinder->freqsByLine(); /*std::ofstream out(outfile); out << "Location\tiHH_0\tiHH_1\tiHS" << std::endl; @@ -165,7 +166,8 @@ void calcIhsMpi( for (const auto& it : res) { auto s = unStd[it.first]; - out2 << it.first << '\t' << hap.lineToId(it.first) << '\t' << s.freq << '\t' << s.iHH_0 << '\t' << s.iHH_1 << '\t' << s.iHS << "\t" << it.second << std::endl; + auto q = freQs[it.first]; + out2 << it.first << '\t' << hap.lineToId(it.first) << '\t' << q << '\t' << s.iHH_0 << '\t' << s.iHH_1 << '\t' << s.iHS << "\t" << it.second << std::endl; } std::cout << "# valid loci: " << res.size() << std::endl; std::cout << "# loci with MAF <= " << minMAF << ": " << ihsfinder->numOutsideMaf() << std::endl;