Method calculateRatioUsingCG_Regression

src/GenomeCopyNumber.cpp:1205–1282

}

void GenomeCopyNumber::calculateRatioUsingCG_Regression( GenomeCopyNumber & controlCopyNumber) {

    int degree = 1;
    bool intercept=0;
    cout << "..polynomial degree set to 1\n";
    cout << "Initial guess for polynomial:\n";

    //first guess about parameters

    double a[MAXDEGREE+1];

    a[degree] = 0;
    a[0]=1;

    cout << "Y = " << a[0] << "*x+" << a[1] << "\n";


    int maximalNumberOfIterations = 300;
    int maximalNumberOfCopies = ploidy_*2;

    vector <float> y; //y ~ a0x+a1
    vector <float> x;
    vector <int> positions;

    //fill x and y:
    vector<ChrCopyNumber>::iterator it;
    for ( it=chrCopyNumber_.begin() ; it != chrCopyNumber_.end(); it++ ) {
        if (! (it->getChromosome().find("X")!=string::npos || it->getChromosome().find("Y")!=string::npos)) {
            vector <float> controlcounts = controlCopyNumber.getChrCopyNumber(it->getChromosome()).getRatio() ;
            //check that everything is all right:
            if (int(controlcounts.size())!=it->getLength()) {
                cerr << "Possible Error: calculateMedianAround ()\n";
            }
            for (int i = 0; i< it->getLength(); i++) {
                if (it->getRatioAtBin(i)>0 && controlcounts[i]>0) {
                    x.push_back(controlcounts[i]);
                    y.push_back(it->getRatioAtBin(i));
                    positions.push_back(it->getCoordinateAtBin(i));
                }
            }
            controlcounts.clear();
        }
    }
    if (degree==1) {
        const char * nametmp = "/home/vboeva/NAS_public/data/projects/FREEC/TEST_FREEC_CLBMA_WES/output_WES/xy.txt";
        std::ofstream file;
        file.open(nametmp);
        for ( unsigned int i=0 ;i<x.size(); i++ ) {
             file << positions[i]<<"\t"<<x[i] <<"\t" << y[i] <<"\n" ;
        }
        file.close();
    }
    float rmserror = runEM(x,y,a,degree,maximalNumberOfIterations,ploidy_,maximalNumberOfCopies, intercept, normalContamination_);
    if (rmserror == -1) {
        cerr << "Error in EM => unable to calculate normalized profile\n";
        return;
    }
    cout << "root mean square error = " << rmserror << "\n";