062add94-a5b5-4fd4-a562-e213ecc8da7c

Committed 20 Jan 2025 05:00PM UTC coverage: 86.049% (-0.1%) from 86.156%

Build # 062add94-a5b5-4fd4-a562-e213ecc8da7c

Build Type

push

circleci

Committed by

web-flow

Commit Message

Merge pull request #377 from DEploid-dev/shajoezhu-patch-1

Update config.yml

Run Details

5545 of 6444 relevant lines covered (86.05%)

5298085.9 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

16.46

/src/export/dEploidIOExport.cpp

/*
 * dEploid is used for deconvoluting Plasmodium falciparum genome from
 * mix-infected patient sample.
 *
 * Copyright (C) 2016-2017 University of Oxford
 *
 * Author: Sha (Joe) Zhu
 *
 * This file is part of dEploid.
 *
 * dEploid is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include "dEploidIO.hpp"
#include "mcmc.hpp"


void DEploidIO::wrapUp() {
    this->writeRecombProb( panel );

    // Get End time before writing the log
    this->getTime(false);

    if (this->useBestPractice()) {
        this->writeChooseKProportion();
    }

    this->writeLog (&std::cout);

    ofstreamExportTmp.open(strExportLog.c_str(),
        ios::out | ios::app | ios::binary);
    this->writeLog (&ofstreamExportTmp);
    ofstreamExportTmp.close();
}


void DEploidIO::writeRecombProb(Panel * panel) {
    if (!doExportRecombProb()) return;

    if ( panel != NULL ) {
        ofstreamExportTmp.open(strExportRecombProb.c_str(),
            ios::out | ios::app | ios::binary );
        ofstreamExportTmp << "p.recomb"       << "\t"
                                 << "p.each"         << "\t"
                                 << "p.no.recomb"    << "\t"
                                 << "p.rec.rec"      << "\t"
                                 << "p.rec.norec"    << "\t"
                                 << "p.norec.norec"  << "\n";
        for ( size_t i = 0; i < panel->pRec_.size(); i++ ) {
            ofstreamExportTmp << panel->pRec_[i]           << "\t"
                                     << panel->pRecEachHap_[i]    << "\t"
                                     << panel->pNoRec_[i]         << "\t"
                                     << panel->pRecRec_[i]        << "\t"
                                     << panel->pRecNoRec_[i]      << "\t"
                                     << panel->pNoRecNoRec_[i]    << "\n";
        }
        ofstreamExportTmp.close();
    }
}



void DEploidIO::writeLog ( ostream * writeTo ) {
    size_t nHash = 30 + string(VERSION).size();
    for ( size_t i = 0; i < nHash; i++) {
        (*writeTo) << "#";
    }
    (*writeTo) << "\n";
    (*writeTo) << "#        dEploid "<< setw(10) << VERSION << " log        #\n";
    for ( size_t i = 0; i < nHash; i++) {
        (*writeTo) << "#";
    }
    (*writeTo) << "\n";
    (*writeTo) << "Program was compiled on: " << compileTime_ << endl;
    (*writeTo) << "dEploid version: " << dEploidGitVersion_ << endl;
    (*writeTo) << "lasso version: " << lassoGitVersion_ << endl;
    (*writeTo) << "\n";
    (*writeTo) << "Input data: \n";
    if (panelFileName_.size() > 0) {
        (*writeTo) << setw(12) << "Panel: "     << panelFileName_  << "\n";
    }
    (*writeTo) << setw(12) << "PLAF: "      << plafFileName_   << "\n";
    if ( useVcf() ) (*writeTo) << setw(12) << "VCF: " << vcfFileName_    << "\n";
    if ( refFileName_.size()>0) (*writeTo) << setw(12) << "REF count: " << refFileName_    << "\n";
    if ( altFileName_.size()>0) (*writeTo) << setw(12) << "ALT count: " << altFileName_    << "\n";
    if ( excludeSites() ) { (*writeTo) << setw(12) << "Exclude: " << excludeFileName_    << "\n"; }
    (*writeTo) << "\n";
    if ( (this->doLsPainting() == false) & (this->doIbdPainting() == false) ) {
        (*writeTo) << "MCMC parameters: "<< "\n";
        (*writeTo) << setw(19) << " MCMC burn: " << mcmcBurn_.getValue() << "\n";
        (*writeTo) << setw(19) << " MCMC sample: " << nMcmcSample_.getValue() << "\n";
        (*writeTo) << setw(19) << " MCMC sample rate: " << mcmcMachineryRate_.getValue() <<"\n";
        (*writeTo) << setw(19) << " Random seed: " << this->randomSeed_.getValue() << "\n";
        if (this->useIBD()) {
            (*writeTo) << setw(19) << "  IBD Method used: YES" << "\n";
        }
        (*writeTo) << setw(19) << " Update Prop: "
                   << (this->doUpdateProp()  ? "YES":"NO") << "\n";
        (*writeTo) << setw(19) << " Update Single: "
                   << (this->doUpdateSingle()? "YES":"NO") << "\n";
        (*writeTo) << setw(19) << " Update Pair: "
                   << (this->doUpdatePair()  ? "YES":"NO") << "\n";
        (*writeTo) << "\n";
    }
    (*writeTo) << "Other parameters:"<< "\n";
    if ( forbidCopyFromSame_ ) { (*writeTo) << " Update pair haplotypes move forbid copying from the same strain!!! \n"; }
    (*writeTo) << setw(20) << " Miss copy prob: "   << this->missCopyProb_.getValue() << "\n";
    (*writeTo) << setw(20) << " Avrg Cent Morgan: " << this->averageCentimorganDistance_ << "\n";
    (*writeTo) << setw(20) << " G: "                << this->parameterG() << "\n";
    if (this->useIBD()) {
    (*writeTo) << setw(20) << " IBD sigma: "        << this->ibdSigma() << "\n";
    } else {
    (*writeTo) << setw(20) << " sigma: "            << this->parameterSigma_.getValue() << "\n";
    }
    (*writeTo) << setw(20) << " ScalingFactor: "    << this->scalingFactor() << "\n";
    (*writeTo) << setw(20) << " VQSLOD:        "    << this->vqslod() << "\n";
    if ( this->initialPropWasGiven() ) {
        (*writeTo) << setw(20) << " Initial prob: " ;
        for ( size_t i = 0; i < this->initialProp.size(); i++ ) {
            (*writeTo) << this->initialProp[i]
                       << ( ( i != (this->kStrain_.getValue()-1) ) ? " " : "\n" );
        }
    }
    (*writeTo) << "\n";
    if ((this->useLasso() == false) & (this->doLsPainting() == false) & (this->doIbdPainting() == false) & (this->doComputeLLK() == false) ) {
        (*writeTo) << "MCMC diagnostic:"<< "\n";
        (*writeTo) << setw(19) << " Accept_ratio: " << acceptRatio_ << "\n";
        (*writeTo) << setw(19) << " Max_llks: " << maxLLKs_ << "\n";
        (*writeTo) << setw(19) << " Final_theta_llks: " << meanThetallks_ << "\n";
        (*writeTo) << setw(19) << " Mean_llks: " << meanllks_ << "\n";
        (*writeTo) << setw(19) << " Stdv_llks: " << stdvllks_ << "\n";
        (*writeTo) << setw(19) << " DIC_by_Dtheta: " << dicByTheta_ << "\n";
        (*writeTo) << setw(19) << " DIC_by_varD: " << dicByVar_ << "\n";
        (*writeTo) << "\n";
    }
    (*writeTo) << "Run time:\n";
    (*writeTo) << setw(14) << "Start at: "  << startingTime_  ;
    (*writeTo) << setw(14) << "End at: "    << endTime_  ;
    (*writeTo) << "\n";
    if (this->useBestPractice()) {
        (*writeTo) << "Confidence: " << endl;
        (*writeTo) << setw(14) << "k: "
                   << this->chooseK.confidence() * 100 << "%" << endl;
        (*writeTo) << "\n";
    }

    if ( this->doComputeLLK() ) {
        (*writeTo) << "Input likelihood: " << llkFromInitialHap_;
        (*writeTo) << "\n";
    } else {
        (*writeTo) << "Output saved to:\n";
        if ( this->doLsPainting() ) {
            for ( size_t i = 0; i < kStrain_.getValue(); i++ ) {
                (*writeTo) << "Posterior probability of strain " << i << ": "<< strExportSingleFwdProbPrefix << i <<endl;
            }
        } else if (this->doIbdPainting()) {
            if (this->ibdProbsIntegrated.size()>1) {
                (*writeTo) << setw(14) << "IBD probs: "  << strIbdExportProbs  << "\n\n";
                (*writeTo) << " IBD probabilities:\n";
                for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ) {
                    (*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";
                }
            }
        } else {
            // if (this->useLasso() == false) {
                // (*writeTo) << setw(14) << "Likelihood: "  << strExportLLK  << "\n";
                // (*writeTo) << setw(14) << "Proportions: " << strExportProp << "\n";
            // }
            // (*writeTo) << setw(14) << "Haplotypes: "  << strExportHap  << "\n";
            if ( doExportVcf() ) { (*writeTo) << setw(14) << "Vcf: "  << strExportVcf  << "\n"; }
            if (this->useIBD()) {
                (*writeTo) << " IBD method output saved to:\n";
                (*writeTo) << setw(14) << "Likelihood: "  << this->prefix_ << ".*.llk\n";
                (*writeTo) << setw(14) << "Proportions: " << this->prefix_ << ".*.prop\n";
                (*writeTo) << setw(14) << "Haplotypes: "  << this->prefix_ << ".*.hap\n";
            } else if(this->useBestPractice()) {
              (*writeTo) << setw(14) << "ChooseK Haps: "  << this->prefix_ << ".chooseK.hap\n";
              (*writeTo) << setw(14) << "IBD Haps: "  << this->prefix_ << ".ibd.hap\n";
              (*writeTo) << setw(14) << "Final Haps: "  << this->prefix_ << ".final.hap\n";
              (*writeTo) << setw(14) << "IBD llks: "  << this->prefix_ << ".ibd.llk\n";
              (*writeTo) << setw(14) << "Proportions: " << this->prefix_ << ".chooseK.prop\n";
              (*writeTo) << setw(14) << "IBD props: " << this->prefix_ << ".ibd.prop\n";
            } else {
              (*writeTo) << setw(14) << "Likelihood: "  << this->prefix_ << ".llk\n";
              (*writeTo) << setw(14) << "Proportions: " << this->prefix_ << ".prop\n";
              (*writeTo) << setw(14) << "Haplotypes: "  << this->prefix_ << ".hap\n";
            }
            if (this->ibdProbsIntegrated.size()>1) {
                (*writeTo) << setw(14) << "IBD probs: "  << strIbdExportProbs  << "\n\n";
                (*writeTo) << " IBD probabilities:\n";
                for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ) {
                    (*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";
                }
            }
        }
        (*writeTo) << "\n";
        if (this->useLasso() == false) {
            (*writeTo) << " IBD best path llk: " << ibdLLK_ << "\n\n";
        }
        this->computeEffectiveKstrain(this->finalProp);
        (*writeTo) << "         Effective_K: " << this->effectiveKstrain_ <<"\n";
        this->computeInferredKstrain(this->finalProp);
        (*writeTo) << "          Inferred_K: " << this->inferredKstrain_ <<"\n";
        this->computeAdjustedEffectiveKstrain();
        (*writeTo) << "Adjusted_effective_K: " << this->adjustedEffectiveKstrain_ <<"\n";
    }
    (*writeTo) << "\n";
    (*writeTo) << "Proportions:\n";
    for ( size_t ii = 0; ii < this->finalProp.size(); ii++) {
        (*writeTo) << setw(10) << this->finalProp[ii];
        (*writeTo) << ((ii < (this->finalProp.size()-1)) ? "\t" : "\n") ;
    }
}


void DEploidIO::writeEventCount() {
    ofstreamExportTmp.open(strExportExtra.c_str(),
        ios::out | ios::app | ios::binary);

    // HEADER
    ofstreamExportTmp << "CHROM" << "\t"
                      << "POS" << "\t"
                      << "IBDpathChangeAt" << "\t"
                      << "finalIBDpathChangeAt" << "\t"

                      << "siteOfTwoSwitchOne" << "\t"
                      << "finalSiteOfTwoSwitchOne" << "\t"

                      << "siteOfTwoMissCopyOne" << "\t"
                      << "finalSiteOfTwoMissCopyOne" << "\t"

                      << "siteOfTwoSwitchTwo" << "\t"
                      << "finalSiteOfTwoSwitchTwo" << "\t"

                      << "siteOfTwoMissCopyTwo" << "\t"
                      << "finalSiteOfTwoMissCopyTwo" << "\t"

                      << "siteOfOneSwitchOne" << "\t"
                      << "finalSiteOfOneSwitchOne" << "\t"

                      << "siteOfOneMissCopyOne" << "\t"
                      << "finalSiteOfOneMissCopyOne" << endl;

    size_t siteIndex = 0;
    for ( size_t chromI = 0; chromI < chrom_.size(); chromI++ ) {
        for ( size_t posI = 0; posI < position_[chromI].size(); posI++) {
            ofstreamExportTmp << chrom_[chromI] << "\t"
                          << (int)position_[chromI][posI] << "\t"

                          << this->IBDpathChangeAt[siteIndex] << "\t"
                          << this->finalIBDpathChangeAt[siteIndex] << "\t"

                          << this->siteOfTwoSwitchOne[siteIndex] << "\t"
                          << this->finalSiteOfTwoSwitchOne[siteIndex] << "\t"

                          << this->siteOfTwoMissCopyOne[siteIndex] << "\t"
                          << this->finalSiteOfTwoMissCopyOne[siteIndex] << "\t"

                          << this->siteOfTwoSwitchTwo[siteIndex] << "\t"
                          << this->finalSiteOfTwoSwitchTwo[siteIndex] << "\t"

                          << this->siteOfTwoMissCopyTwo[siteIndex] << "\t"
                          << this->finalSiteOfTwoMissCopyTwo[siteIndex] << "\t"

                          << this->siteOfOneSwitchOne[siteIndex] << "\t"
                          << this->finalSiteOfOneSwitchOne[siteIndex] << "\t"

                          << this->siteOfOneMissCopyOne[siteIndex] << "\t"
                          << this->finalSiteOfOneMissCopyOne[siteIndex] << endl;
            siteIndex++;
        }
    }

    assert(siteIndex == this->IBDpathChangeAt.size());
    ofstreamExportTmp.close();
}


void DEploidIO::writeIBDpostProb(vector < vector <double> > & reshapedProbs,
    vector <string> header) {
    ostream * writeTo;
    #ifdef UNITTEST
    writeTo = &std::cout;
    #endif

    #ifndef UNITTEST
    ofstreamExportTmp.open(strIbdExportProbs.c_str(),
        ios::out | ios::app | ios::binary);
    writeTo = &ofstreamExportTmp;
    #endif

    (*writeTo) << "CHROM" << "\t" << "POS" << "\t";
    for (string tmp : header) {
        (*writeTo) << tmp << ((tmp!=header[header.size()-1])?"\t":"\n");
    }

    size_t siteIndex = 0;
    for ( size_t chromIndex = 0; chromIndex < position_.size(); chromIndex++) {
        for ( size_t posI = 0; posI < position_[chromIndex].size(); posI++) {
            (*writeTo) << chrom_[chromIndex] << "\t"
                       << (int)position_[chromIndex][posI] << "\t";
            for (size_t ij = 0; ij < reshapedProbs[siteIndex].size(); ij++) {
                (*writeTo) << reshapedProbs[siteIndex][ij] << "\t";
            }
            (*writeTo) << endl;
            siteIndex++;
        }
    }
    assert(siteIndex == nLoci());

    #ifndef UNITTEST
        ofstreamExportTmp.close();
    #endif
}


void DEploidIO::writeIBDviterbi(vector <size_t> & viterbiState) {
    ostream * writeTo;
    #ifdef UNITTEST
    writeTo = &std::cout;
    #endif

    #ifndef UNITTEST
    ofstreamExportTmp.open( strIbdExportViterbi.c_str(),
    ios::out | ios::app | ios::binary );
    writeTo = &ofstreamExportTmp;
    #endif

    (*writeTo) << "CHROM" << "\t" << "POS" << "\t" << "viterbi" << "\n" ;

    size_t siteIndex = 0;
    for ( size_t chromIndex = 0; chromIndex < position_.size(); chromIndex++) {
        for ( size_t posI = 0; posI < position_[chromIndex].size(); posI++) {
            (*writeTo) << chrom_[chromIndex] << "\t"
                       << (int)position_[chromIndex][posI] << "\t"
                       << (int)viterbiState[siteIndex] << "\n";
            siteIndex++;
        }
    }
    assert(siteIndex == nLoci());

    #ifndef UNITTEST
    ofstreamExportTmp.close();
    #endif
}


void DEploidIO::operation_paintIBD() {
    vector <double> goodProp;
    vector <size_t> goodStrainIdx;

    if ( this->doIbdPainting() ) {
        this->finalProp = this->initialProp;
    }

    for ( size_t i = 0; i < this->finalProp.size(); i++) {
        if (this->finalProp[i] > 0.01) {
            goodProp.push_back(this->finalProp[i]);
            goodStrainIdx.push_back(i);
        }
    }

    if (goodProp.size() == 1) {
        return;
    }

    DEploidIO tmpDEploidIO; // (*this);
    tmpDEploidIO.kStrain_.init(goodProp.size());
    tmpDEploidIO.setInitialPropWasGiven(true);
    tmpDEploidIO.initialProp = goodProp;
    tmpDEploidIO.finalProp = goodProp;
    tmpDEploidIO.refCount_ = this->refCount_;
    tmpDEploidIO.altCount_ = this->altCount_;
    tmpDEploidIO.plaf_ = this->plaf_;
    tmpDEploidIO.nLoci_ = this->nLoci();
    tmpDEploidIO.position_ = this->position_;
    tmpDEploidIO.chrom_ = this->chrom_;
    tmpDEploidIO.setParameterG(this->parameterG());
    // tmpDEploidIO.useConstRecomb_ = true;
    // tmpDEploidIO.constRecombProb_ = 0.000001;

    // tmpDEploidIO.writeLog (&std::cout);

    MersenneTwister tmpRg(this->randomSeed_.getValue());
    IBDpath tmpIBDpath;
    tmpIBDpath.init(tmpDEploidIO, &tmpRg);
    tmpIBDpath.buildPathProbabilityForPainting(goodProp);
    this->ibdLLK_ = tmpIBDpath.bestPath(goodProp);
    this->ibdProbsHeader = tmpIBDpath.getIBDprobsHeader();
    this->getIBDprobsIntegrated(tmpIBDpath.fwdbwd);


    for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ) {
        cout << setw(14) << this->ibdProbsHeader[stateI] << ": "
             << this->ibdProbsIntegrated[stateI] << "\n";
    }

    this->writeIBDpostProb(tmpIBDpath.fwdbwd, this->ibdProbsHeader);
}


void DEploidIO::operation_paintIBDviterbi() {
    vector <double> goodProp;
    vector <size_t> goodStrainIdx;

    if (this->doIbdViterbiPainting()) {
        this->finalProp = this->initialProp;
    }

    for (size_t i = 0; i < this->finalProp.size(); i++) {
        if (this->finalProp[i] > 0.01) {
            goodProp.push_back(this->finalProp[i]);
            goodStrainIdx.push_back(i);
        }
    }

    if (goodProp.size() == 1) {
        return;
    }

    DEploidIO tmpDEploidIO;  // (*this);
    tmpDEploidIO.kStrain_.init(goodProp.size());
    tmpDEploidIO.setInitialPropWasGiven(true);
    tmpDEploidIO.initialProp = goodProp;
    tmpDEploidIO.finalProp = goodProp;
    tmpDEploidIO.refCount_ = this->refCount_;
    tmpDEploidIO.altCount_ = this->altCount_;
    tmpDEploidIO.plaf_ = this->plaf_;
    tmpDEploidIO.nLoci_ = this->nLoci();
    tmpDEploidIO.position_ = this->position_;
    tmpDEploidIO.chrom_ = this->chrom_;
    tmpDEploidIO.setParameterG(this->parameterG());
    // tmpDEploidIO.useConstRecomb_ = true;
    // tmpDEploidIO.constRecombProb_ = 0.000001;

    // tmpDEploidIO.writeLog (&std::cout);

    MersenneTwister tmpRg(this->randomSeed_.getValue());
    IBDpath tmpIBDpath;
    tmpIBDpath.init(tmpDEploidIO, &tmpRg);

    // tmpIBDpath.buildPathProbabilityForPainting(goodProp);
    // tmpIBDpath.findViterbiPath(goodProp);
    this->ibdLLK_ = tmpIBDpath.findViterbiPath(goodProp);
    // ''tmpIBDpath.bestPath(goodProp);
    // this->ibdProbsHeader = tmpIBDpath.getIBDprobsHeader();
    // this->getIBDprobsIntegrated(tmpIBDpath.fm);

    this->writeIBDviterbi(tmpIBDpath.viterbiPath);
}


void DEploidIO::writeChooseKProportion() {
    string fileName = this->prefix_ + ".chooseK.prop";
    ofstreamExportTmp.open(fileName.c_str(), ios::out | ios::binary);
    ostream * writeTo = &ofstreamExportTmp;
    for (size_t i = 0; i < this->chooseK.proportions_.size(); i++) {
        for (size_t ii = 0; ii < this->chooseK.proportions_[i].size(); ii++) {
            (*writeTo)  << this->chooseK.proportions_[i][ii];
            (*writeTo)  << (ii == (this->chooseK.proportions_[i].size()-1) ?
                        "\n":"\t");
        }
    }
    ofstreamExportTmp.close();
}

1	/*
2	* dEploid is used for deconvoluting Plasmodium falciparum genome from
3	* mix-infected patient sample.
4	*
5	* Copyright (C) 2016-2017 University of Oxford
6	*
7	* Author: Sha (Joe) Zhu
8	*
9	* This file is part of dEploid.
10	*
11	* dEploid is free software: you can redistribute it and/or modify
12	* it under the terms of the GNU General Public License as published by
13	* the Free Software Foundation, either version 3 of the License, or
14	* (at your option) any later version.
15	*
16	* This program is distributed in the hope that it will be useful,
17	* but WITHOUT ANY WARRANTY; without even the implied warranty of
18	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19	* GNU General Public License for more details.
20	*
21	* You should have received a copy of the GNU General Public License
22	* along with this program. If not, see <http://www.gnu.org/licenses/>.
23	*
24	*/
25
26	#include "dEploidIO.hpp"
27	#include "mcmc.hpp"
28
29
30	void DEploidIO::wrapUp() {	×
31	this->writeRecombProb( panel );	×
32
33	// Get End time before writing the log
34	this->getTime(false);	×
35
36	if (this->useBestPractice()) {	×
37	this->writeChooseKProportion();	×
38	}
39
40	this->writeLog (&std::cout);	×
41
42	ofstreamExportTmp.open(strExportLog.c_str(),	×
43	ios::out \| ios::app \| ios::binary);
44	this->writeLog (&ofstreamExportTmp);	×
45	ofstreamExportTmp.close();	×
46	}	×
47
48
49	void DEploidIO::writeRecombProb(Panel * panel) {	×
50	if (!doExportRecombProb()) return;	×
51
52	if ( panel != NULL ) {	×
53	ofstreamExportTmp.open(strExportRecombProb.c_str(),	×
54	ios::out \| ios::app \| ios::binary );
55	ofstreamExportTmp << "p.recomb" << "\t"
56	<< "p.each" << "\t"
57	<< "p.no.recomb" << "\t"
58	<< "p.rec.rec" << "\t"
59	<< "p.rec.norec" << "\t"
60	<< "p.norec.norec" << "\n";	×
61	for ( size_t i = 0; i < panel->pRec_.size(); i++ ) {	×
62	ofstreamExportTmp << panel->pRec_[i] << "\t"	×
63	<< panel->pRecEachHap_[i] << "\t"	×
64	<< panel->pNoRec_[i] << "\t"	×
65	<< panel->pRecRec_[i] << "\t"	×
66	<< panel->pRecNoRec_[i] << "\t"	×
67	<< panel->pNoRecNoRec_[i] << "\n";	×
68	}
69	ofstreamExportTmp.close();	×
70	}
71	}
72
73
74
75	void DEploidIO::writeLog ( ostream * writeTo ) {	×
76	size_t nHash = 30 + string(VERSION).size();	×
77	for ( size_t i = 0; i < nHash; i++) {	×
78	(*writeTo) << "#";	×
79	}
80	(*writeTo) << "\n";	×
81	(*writeTo) << "# dEploid "<< setw(10) << VERSION << " log #\n";	×
82	for ( size_t i = 0; i < nHash; i++) {	×
83	(*writeTo) << "#";	×
84	}
85	(*writeTo) << "\n";	×
86	(*writeTo) << "Program was compiled on: " << compileTime_ << endl;
87	(*writeTo) << "dEploid version: " << dEploidGitVersion_ << endl;
88	(*writeTo) << "lasso version: " << lassoGitVersion_ << endl;
89	(*writeTo) << "\n";	×
90	(*writeTo) << "Input data: \n";	×
91	if (panelFileName_.size() > 0) {	×
92	(*writeTo) << setw(12) << "Panel: " << panelFileName_ << "\n";	×
93	}
94	(*writeTo) << setw(12) << "PLAF: " << plafFileName_ << "\n";	×
95	if ( useVcf() ) (*writeTo) << setw(12) << "VCF: " << vcfFileName_ << "\n";	×
96	if ( refFileName_.size()>0) (*writeTo) << setw(12) << "REF count: " << refFileName_ << "\n";	×
97	if ( altFileName_.size()>0) (*writeTo) << setw(12) << "ALT count: " << altFileName_ << "\n";	×
98	if ( excludeSites() ) { (*writeTo) << setw(12) << "Exclude: " << excludeFileName_ << "\n"; }	×
99	(*writeTo) << "\n";	×
100	if ( (this->doLsPainting() == false) & (this->doIbdPainting() == false) ) {	×
101	(*writeTo) << "MCMC parameters: "<< "\n";	×
102	(*writeTo) << setw(19) << " MCMC burn: " << mcmcBurn_.getValue() << "\n";	×
103	(*writeTo) << setw(19) << " MCMC sample: " << nMcmcSample_.getValue() << "\n";	×
104	(*writeTo) << setw(19) << " MCMC sample rate: " << mcmcMachineryRate_.getValue() <<"\n";	×
105	(*writeTo) << setw(19) << " Random seed: " << this->randomSeed_.getValue() << "\n";	×
106	if (this->useIBD()) {	×
107	(*writeTo) << setw(19) << " IBD Method used: YES" << "\n";	×
108	}
109	(*writeTo) << setw(19) << " Update Prop: "
110	<< (this->doUpdateProp() ? "YES":"NO") << "\n";	×
111	(*writeTo) << setw(19) << " Update Single: "
112	<< (this->doUpdateSingle()? "YES":"NO") << "\n";	×
113	(*writeTo) << setw(19) << " Update Pair: "
114	<< (this->doUpdatePair() ? "YES":"NO") << "\n";	×
115	(*writeTo) << "\n";	×
116	}
117	(*writeTo) << "Other parameters:"<< "\n";	×
118	if ( forbidCopyFromSame_ ) { (*writeTo) << " Update pair haplotypes move forbid copying from the same strain!!! \n"; }	×
119	(*writeTo) << setw(20) << " Miss copy prob: " << this->missCopyProb_.getValue() << "\n";	×
120	(*writeTo) << setw(20) << " Avrg Cent Morgan: " << this->averageCentimorganDistance_ << "\n";	×
121	(*writeTo) << setw(20) << " G: " << this->parameterG() << "\n";	×
122	if (this->useIBD()) {	×
123	(*writeTo) << setw(20) << " IBD sigma: " << this->ibdSigma() << "\n";	×
124	} else {
125	(*writeTo) << setw(20) << " sigma: " << this->parameterSigma_.getValue() << "\n";	×
126	}
127	(*writeTo) << setw(20) << " ScalingFactor: " << this->scalingFactor() << "\n";	×
128	(*writeTo) << setw(20) << " VQSLOD: " << this->vqslod() << "\n";	×
129	if ( this->initialPropWasGiven() ) {	×
130	(*writeTo) << setw(20) << " Initial prob: " ;	×
131	for ( size_t i = 0; i < this->initialProp.size(); i++ ) {	×
132	(*writeTo) << this->initialProp[i]	×
133	<< ( ( i != (this->kStrain_.getValue()-1) ) ? " " : "\n" );	×
134	}
135	}
136	(*writeTo) << "\n";	×
137	if ((this->useLasso() == false) & (this->doLsPainting() == false) & (this->doIbdPainting() == false) & (this->doComputeLLK() == false) ) {	×
138	(*writeTo) << "MCMC diagnostic:"<< "\n";	×
139	(*writeTo) << setw(19) << " Accept_ratio: " << acceptRatio_ << "\n";	×
140	(*writeTo) << setw(19) << " Max_llks: " << maxLLKs_ << "\n";	×
141	(*writeTo) << setw(19) << " Final_theta_llks: " << meanThetallks_ << "\n";	×
142	(*writeTo) << setw(19) << " Mean_llks: " << meanllks_ << "\n";	×
143	(*writeTo) << setw(19) << " Stdv_llks: " << stdvllks_ << "\n";	×
144	(*writeTo) << setw(19) << " DIC_by_Dtheta: " << dicByTheta_ << "\n";	×
145	(*writeTo) << setw(19) << " DIC_by_varD: " << dicByVar_ << "\n";	×
146	(*writeTo) << "\n";	×
147	}
148	(*writeTo) << "Run time:\n";	×
149	(*writeTo) << setw(14) << "Start at: " << startingTime_ ;
150	(*writeTo) << setw(14) << "End at: " << endTime_ ;
151	(*writeTo) << "\n";	×
152	if (this->useBestPractice()) {	×
153	(*writeTo) << "Confidence: " << endl;
154	(*writeTo) << setw(14) << "k: "
155	<< this->chooseK.confidence() * 100 << "%" << endl;	×
156	(*writeTo) << "\n";	×
157	}
158
159	if ( this->doComputeLLK() ) {	×
160	(*writeTo) << "Input likelihood: " << llkFromInitialHap_;	×
161	(*writeTo) << "\n";	×
162	} else {
163	(*writeTo) << "Output saved to:\n";	×
164	if ( this->doLsPainting() ) {	×
165	for ( size_t i = 0; i < kStrain_.getValue(); i++ ) {	×
166	(*writeTo) << "Posterior probability of strain " << i << ": "<< strExportSingleFwdProbPrefix << i <<endl;
167	}
168	} else if (this->doIbdPainting()) {	×
169	if (this->ibdProbsIntegrated.size()>1) {	×
170	(*writeTo) << setw(14) << "IBD probs: " << strIbdExportProbs << "\n\n";	×
171	(*writeTo) << " IBD probabilities:\n";	×
172	for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ) {	×
173	(*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";	×
174	}
175	}
176	} else {
177	// if (this->useLasso() == false) {
178	// (*writeTo) << setw(14) << "Likelihood: " << strExportLLK << "\n";
179	// (*writeTo) << setw(14) << "Proportions: " << strExportProp << "\n";
180	// }
181	// (*writeTo) << setw(14) << "Haplotypes: " << strExportHap << "\n";
182	if ( doExportVcf() ) { (*writeTo) << setw(14) << "Vcf: " << strExportVcf << "\n"; }	×
183	if (this->useIBD()) {	×
184	(*writeTo) << " IBD method output saved to:\n";	×
185	(writeTo) << setw(14) << "Likelihood: " << this->prefix_ << "..llk\n";	×
186	(writeTo) << setw(14) << "Proportions: " << this->prefix_ << "..prop\n";	×
187	(writeTo) << setw(14) << "Haplotypes: " << this->prefix_ << "..hap\n";	×
188	} else if(this->useBestPractice()) {	×
189	(*writeTo) << setw(14) << "ChooseK Haps: " << this->prefix_ << ".chooseK.hap\n";	×
190	(*writeTo) << setw(14) << "IBD Haps: " << this->prefix_ << ".ibd.hap\n";	×
191	(*writeTo) << setw(14) << "Final Haps: " << this->prefix_ << ".final.hap\n";	×
192	(*writeTo) << setw(14) << "IBD llks: " << this->prefix_ << ".ibd.llk\n";	×
193	(*writeTo) << setw(14) << "Proportions: " << this->prefix_ << ".chooseK.prop\n";	×
194	(*writeTo) << setw(14) << "IBD props: " << this->prefix_ << ".ibd.prop\n";	×
195	} else {
196	(*writeTo) << setw(14) << "Likelihood: " << this->prefix_ << ".llk\n";	×
197	(*writeTo) << setw(14) << "Proportions: " << this->prefix_ << ".prop\n";	×
198	(*writeTo) << setw(14) << "Haplotypes: " << this->prefix_ << ".hap\n";	×
199	}
200	if (this->ibdProbsIntegrated.size()>1) {	×
201	(*writeTo) << setw(14) << "IBD probs: " << strIbdExportProbs << "\n\n";	×
202	(*writeTo) << " IBD probabilities:\n";	×
203	for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ) {	×
204	(*writeTo) << setw(14) << this->ibdProbsHeader[stateI] << ": " << this->ibdProbsIntegrated[stateI] << "\n";	×
205	}
206	}
207	}
208	(*writeTo) << "\n";	×
209	if (this->useLasso() == false) {	×
210	(*writeTo) << " IBD best path llk: " << ibdLLK_ << "\n\n";	×
211	}
212	this->computeEffectiveKstrain(this->finalProp);	×
213	(*writeTo) << " Effective_K: " << this->effectiveKstrain_ <<"\n";	×
214	this->computeInferredKstrain(this->finalProp);	×
215	(*writeTo) << " Inferred_K: " << this->inferredKstrain_ <<"\n";	×
216	this->computeAdjustedEffectiveKstrain();	×
217	(*writeTo) << "Adjusted_effective_K: " << this->adjustedEffectiveKstrain_ <<"\n";	×
218	}
219	(*writeTo) << "\n";	×
220	(*writeTo) << "Proportions:\n";	×
221	for ( size_t ii = 0; ii < this->finalProp.size(); ii++) {	×
222	(*writeTo) << setw(10) << this->finalProp[ii];	×
223	(*writeTo) << ((ii < (this->finalProp.size()-1)) ? "\t" : "\n") ;	×
224	}
225	}	×
226
227
228	void DEploidIO::writeEventCount() {	×
229	ofstreamExportTmp.open(strExportExtra.c_str(),	×
230	ios::out \| ios::app \| ios::binary);
231
232	// HEADER
233	ofstreamExportTmp << "CHROM" << "\t"
234	<< "POS" << "\t"
235	<< "IBDpathChangeAt" << "\t"
236	<< "finalIBDpathChangeAt" << "\t"
237
238	<< "siteOfTwoSwitchOne" << "\t"
239	<< "finalSiteOfTwoSwitchOne" << "\t"
240
241	<< "siteOfTwoMissCopyOne" << "\t"
242	<< "finalSiteOfTwoMissCopyOne" << "\t"
243
244	<< "siteOfTwoSwitchTwo" << "\t"
245	<< "finalSiteOfTwoSwitchTwo" << "\t"
246
247	<< "siteOfTwoMissCopyTwo" << "\t"
248	<< "finalSiteOfTwoMissCopyTwo" << "\t"
249
250	<< "siteOfOneSwitchOne" << "\t"
251	<< "finalSiteOfOneSwitchOne" << "\t"
252
253	<< "siteOfOneMissCopyOne" << "\t"
254	<< "finalSiteOfOneMissCopyOne" << endl;	×
255
256	size_t siteIndex = 0;
257	for ( size_t chromI = 0; chromI < chrom_.size(); chromI++ ) {	×
258	for ( size_t posI = 0; posI < position_[chromI].size(); posI++) {	×
259	ofstreamExportTmp << chrom_[chromI] << "\t"
260	<< (int)position_[chromI][posI] << "\t"	×
261
262	<< this->IBDpathChangeAt[siteIndex] << "\t"	×
263	<< this->finalIBDpathChangeAt[siteIndex] << "\t"	×
264
265	<< this->siteOfTwoSwitchOne[siteIndex] << "\t"	×
266	<< this->finalSiteOfTwoSwitchOne[siteIndex] << "\t"	×
267
268	<< this->siteOfTwoMissCopyOne[siteIndex] << "\t"	×
269	<< this->finalSiteOfTwoMissCopyOne[siteIndex] << "\t"	×
270
271	<< this->siteOfTwoSwitchTwo[siteIndex] << "\t"	×
272	<< this->finalSiteOfTwoSwitchTwo[siteIndex] << "\t"	×
273
274	<< this->siteOfTwoMissCopyTwo[siteIndex] << "\t"	×
275	<< this->finalSiteOfTwoMissCopyTwo[siteIndex] << "\t"	×
276
277	<< this->siteOfOneSwitchOne[siteIndex] << "\t"	×
278	<< this->finalSiteOfOneSwitchOne[siteIndex] << "\t"	×
279
280	<< this->siteOfOneMissCopyOne[siteIndex] << "\t"	×
281	<< this->finalSiteOfOneMissCopyOne[siteIndex] << endl;	×
282	siteIndex++;	×
283	}
284	}
285
286	assert(siteIndex == this->IBDpathChangeAt.size());
287	ofstreamExportTmp.close();	×
288	}	×
289
290
291	void DEploidIO::writeIBDpostProb(vector < vector <double> > & reshapedProbs,	5✔
292	vector <string> header) {
293	ostream * writeTo;
294	#ifdef UNITTEST
295	writeTo = &std::cout;
296	#endif
297
298	#ifndef UNITTEST
299	ofstreamExportTmp.open(strIbdExportProbs.c_str(),
300	ios::out \| ios::app \| ios::binary);
301	writeTo = &ofstreamExportTmp;
302	#endif
303
304	(*writeTo) << "CHROM" << "\t" << "POS" << "\t";	5✔
305	for (string tmp : header) {	28✔
306	(*writeTo) << tmp << ((tmp!=header[header.size()-1])?"\t":"\n");	28✔
307	}
308
309	size_t siteIndex = 0;
310	for ( size_t chromIndex = 0; chromIndex < position_.size(); chromIndex++) {	75✔
311	for ( size_t posI = 0; posI < position_[chromIndex].size(); posI++) {	2,685✔
312	(*writeTo) << chrom_[chromIndex] << "\t"
313	<< (int)position_[chromIndex][posI] << "\t";	2,615✔
314	for (size_t ij = 0; ij < reshapedProbs[siteIndex].size(); ij++) {	15,567✔
315	(*writeTo) << reshapedProbs[siteIndex][ij] << "\t";	12,952✔
316	}
317	(*writeTo) << endl;
318	siteIndex++;	2,615✔
319	}
320	}
321	assert(siteIndex == nLoci());
322
323	#ifndef UNITTEST
324	ofstreamExportTmp.close();
325	#endif
326	}	5✔
327
328
329	void DEploidIO::writeIBDviterbi(vector <size_t> & viterbiState) {	×
330	ostream * writeTo;
331	#ifdef UNITTEST
332	writeTo = &std::cout;
333	#endif
334
335	#ifndef UNITTEST
336	ofstreamExportTmp.open( strIbdExportViterbi.c_str(),
337	ios::out \| ios::app \| ios::binary );
338	writeTo = &ofstreamExportTmp;
339	#endif
340
341	(*writeTo) << "CHROM" << "\t" << "POS" << "\t" << "viterbi" << "\n" ;	×
342
343	size_t siteIndex = 0;
344	for ( size_t chromIndex = 0; chromIndex < position_.size(); chromIndex++) {	×
345	for ( size_t posI = 0; posI < position_[chromIndex].size(); posI++) {	×
346	(*writeTo) << chrom_[chromIndex] << "\t"
347	<< (int)position_[chromIndex][posI] << "\t"	×
348	<< (int)viterbiState[siteIndex] << "\n";	×
349	siteIndex++;	×
350	}
351	}
352	assert(siteIndex == nLoci());
353
354	#ifndef UNITTEST
355	ofstreamExportTmp.close();
356	#endif
357	}	×
358
359
360	void DEploidIO::operation_paintIBD() {	5✔
361	vector <double> goodProp;
362	vector <size_t> goodStrainIdx;
363
364	if ( this->doIbdPainting() ) {	5✔
365	this->finalProp = this->initialProp;	1✔
366	}
367
368	for ( size_t i = 0; i < this->finalProp.size(); i++) {	20✔
369	if (this->finalProp[i] > 0.01) {	15✔
370	goodProp.push_back(this->finalProp[i]);	12✔
371	goodStrainIdx.push_back(i);	12✔
372	}
373	}
374
375	if (goodProp.size() == 1) {	5✔
376	return;
377	}
378
379	DEploidIO tmpDEploidIO; // (*this);	5✔
380	tmpDEploidIO.kStrain_.init(goodProp.size());
381	tmpDEploidIO.setInitialPropWasGiven(true);
382	tmpDEploidIO.initialProp = goodProp;	5✔
383	tmpDEploidIO.finalProp = goodProp;	5✔
384	tmpDEploidIO.refCount_ = this->refCount_;	5✔
385	tmpDEploidIO.altCount_ = this->altCount_;	5✔
386	tmpDEploidIO.plaf_ = this->plaf_;	5✔
387	tmpDEploidIO.nLoci_ = this->nLoci();	5✔
388	tmpDEploidIO.position_ = this->position_;	5✔
389	tmpDEploidIO.chrom_ = this->chrom_;	5✔
390	tmpDEploidIO.setParameterG(this->parameterG());
391	// tmpDEploidIO.useConstRecomb_ = true;
392	// tmpDEploidIO.constRecombProb_ = 0.000001;
393
394	// tmpDEploidIO.writeLog (&std::cout);
395
396	MersenneTwister tmpRg(this->randomSeed_.getValue());	5✔
397	IBDpath tmpIBDpath;	5✔
398	tmpIBDpath.init(tmpDEploidIO, &tmpRg);	5✔
399	tmpIBDpath.buildPathProbabilityForPainting(goodProp);	5✔
400	this->ibdLLK_ = tmpIBDpath.bestPath(goodProp);	5✔
401	this->ibdProbsHeader = tmpIBDpath.getIBDprobsHeader();	5✔
402	this->getIBDprobsIntegrated(tmpIBDpath.fwdbwd);	5✔
403
404
405	for ( size_t stateI = 0; stateI < this->ibdProbsHeader.size(); stateI++ ) {	28✔
406	cout << setw(14) << this->ibdProbsHeader[stateI] << ": "
407	<< this->ibdProbsIntegrated[stateI] << "\n";	23✔
408	}
409
410	this->writeIBDpostProb(tmpIBDpath.fwdbwd, this->ibdProbsHeader);	5✔
411	}	10✔
412
413
414	void DEploidIO::operation_paintIBDviterbi() {	×
415	vector <double> goodProp;
416	vector <size_t> goodStrainIdx;
417
418	if (this->doIbdViterbiPainting()) {	×
419	this->finalProp = this->initialProp;	×
420	}
421
422	for (size_t i = 0; i < this->finalProp.size(); i++) {	×
423	if (this->finalProp[i] > 0.01) {	×
424	goodProp.push_back(this->finalProp[i]);	×
425	goodStrainIdx.push_back(i);	×
426	}
427	}
428
429	if (goodProp.size() == 1) {	×
430	return;
431	}
432
433	DEploidIO tmpDEploidIO; // (*this);	×
434	tmpDEploidIO.kStrain_.init(goodProp.size());
435	tmpDEploidIO.setInitialPropWasGiven(true);
436	tmpDEploidIO.initialProp = goodProp;	×
437	tmpDEploidIO.finalProp = goodProp;	×
438	tmpDEploidIO.refCount_ = this->refCount_;	×
439	tmpDEploidIO.altCount_ = this->altCount_;	×
440	tmpDEploidIO.plaf_ = this->plaf_;	×
441	tmpDEploidIO.nLoci_ = this->nLoci();	×
442	tmpDEploidIO.position_ = this->position_;	×
443	tmpDEploidIO.chrom_ = this->chrom_;	×
444	tmpDEploidIO.setParameterG(this->parameterG());
445	// tmpDEploidIO.useConstRecomb_ = true;
446	// tmpDEploidIO.constRecombProb_ = 0.000001;
447
448	// tmpDEploidIO.writeLog (&std::cout);
449
450	MersenneTwister tmpRg(this->randomSeed_.getValue());	×
451	IBDpath tmpIBDpath;	×
452	tmpIBDpath.init(tmpDEploidIO, &tmpRg);	×
453
454	// tmpIBDpath.buildPathProbabilityForPainting(goodProp);
455	// tmpIBDpath.findViterbiPath(goodProp);
456	this->ibdLLK_ = tmpIBDpath.findViterbiPath(goodProp);	×
457	// ''tmpIBDpath.bestPath(goodProp);
458	// this->ibdProbsHeader = tmpIBDpath.getIBDprobsHeader();
459	// this->getIBDprobsIntegrated(tmpIBDpath.fm);
460
461	this->writeIBDviterbi(tmpIBDpath.viterbiPath);	×
462	}	×
463
464
465	void DEploidIO::writeChooseKProportion() {	×
466	string fileName = this->prefix_ + ".chooseK.prop";	×
467	ofstreamExportTmp.open(fileName.c_str(), ios::out \| ios::binary);	×
468	ostream * writeTo = &ofstreamExportTmp;	×
469	for (size_t i = 0; i < this->chooseK.proportions_.size(); i++) {	×
470	for (size_t ii = 0; ii < this->chooseK.proportions_[i].size(); ii++) {	×
471	(*writeTo) << this->chooseK.proportions_[i][ii];	×
472	(*writeTo) << (ii == (this->chooseK.proportions_[i].size()-1) ?	×
473	"\n":"\t");	×
474	}
475	}
476	ofstreamExportTmp.close();	×
477	}	×

DEploid-dev / DEploid / 062add94-a5b5-4fd4-a562-e213ecc8da7c

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous