Method evaluateContaminationwithLR

src/GenomeCopyNumber.cpp:3452–3525

}

float GenomeCopyNumber::evaluateContaminationwithLR () {
	string::size_type pos = 0;
	vector <float> observed_values;
    vector <float> expected_values;
	map<string,int>::iterator it;

	for ( it=chromosomesInd_.begin() ; it != chromosomesInd_.end(); it++ ) {
		string chrNumber = (*it).first;
		if ( ( pos = chrNumber.find("chr")) != string::npos )
			chrNumber.replace( pos, 3, "" );
        if ( ( pos = chrNumber.find("X")) != string::npos ) 		//exclude X and Y from the analysis
            continue;
        if ( ( pos = chrNumber.find("Y")) != string::npos )
            continue;
		int index = findIndex(chrNumber);
		if (index == NA) {
		    cerr << "An error occurred in GenomeCopyNumber::evaluateContamination: could not find an index for "<<chrNumber<<"\n";
		    return 0;
        }
		int length = chrCopyNumber_[index].getLength();
		for (int i = 0; i< length; i++) {
			float observed = chrCopyNumber_[index].getRatioAtBin(i);
			if (observed!=NA) {
                if (isRatioLogged_)
                    observed=pow(2,observed);
                float expected = observed;
                if (chrCopyNumber_[index].isMedianCalculated()) {
                    expected = round_by_ploidy(chrCopyNumber_[index].getMedianProfileAtI(i),ploidy_) ;
                }
                if (expected!=NA && expected<1+2.0/ploidy_){
                    observed_values.push_back(observed-1);
                    expected_values.push_back(expected-1);
                }
			}
		}
	}

    vector <float> x = expected_values;
    vector <float> y = observed_values;

    int n = x.size();
    float sum_x = 0;
    float sum_y = 0;
    float sum_x2 = 0;
    float sum_y2 = 0;
    float sum_xy = 0;
    for (int i = 0; i < n; i++)     {
        sum_x += x[i];
        sum_y += y[i];
        sum_x2 += x[i] * x[i];
        sum_y2 += y[i] * y[i];
        sum_xy += x[i] * y[i];
    }

    float a = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x * sum_x); //CARINO
    a =sum_xy/sum_x2; //VALENTINA
      //  float b = (sum_x2 * sum_y - sum_x * sum_xy) / (n * sum_x2 - sum_x * sum_x);
    float contam_a =  (1-a)/(1-a+2*a/(float)ploidy_);