| 3393 | } |
| 3394 | |
| 3395 | float GenomeCopyNumber::evaluateContamination () { |
| 3396 | float contam = 0; |
| 3397 | string::size_type pos = 0; |
| 3398 | vector <float> values; |
| 3399 | vector <float> weights; |
| 3400 | map<string,int>::iterator it; |
| 3401 | for ( it=chromosomesInd_.begin() ; it != chromosomesInd_.end(); it++ ) { |
| 3402 | string chrNumber = (*it).first; |
| 3403 | if ( ( pos = chrNumber.find("chr")) != string::npos ) |
| 3404 | chrNumber.replace( pos, 3, "" ); |
| 3405 | if ( ( pos = chrNumber.find("X")) != string::npos ) //exclude X and Y from the analysis |
| 3406 | continue; |
| 3407 | if ( ( pos = chrNumber.find("Y")) != string::npos ) |
| 3408 | continue; |
| 3409 | int index = findIndex(chrNumber); |
| 3410 | if (index == NA) { |
| 3411 | cerr << "An error occurred in GenomeCopyNumber::evaluateContamination: could not find an index for "<<chrNumber<<"\n"; |
| 3412 | return 0; |
| 3413 | } |
| 3414 | int length = chrCopyNumber_[index].getLength(); |
| 3415 | for (int i = 0; i< length; i++) { |
| 3416 | float observed = chrCopyNumber_[index].getRatioAtBin(i); |
| 3417 | if (observed!=NA) { |
| 3418 | if(isRatioLogged_) |
| 3419 | observed=pow(2,observed); |
| 3420 | float expected = observed; |
| 3421 | if (chrCopyNumber_[index].isMedianCalculated()) { |
| 3422 | expected = chrCopyNumber_[index].getMedianProfileAtI(i) ; |
| 3423 | if (chrCopyNumber_[index].isSmoothed() && WESanalysis == false) |
| 3424 | expected = chrCopyNumber_[index].getSmoothedProfileAtI(i); |
| 3425 | } |
| 3426 | if (!(expected == 1 || expected <= 0 || expected >= 1+2.0/ploidy_ || observed > 3 || observed <= 0) |
| 3427 | && (((1>observed)&&(1>expected))||((1<observed)&&(1<expected)))) {// should it be something related to ploidy_ and not 2 |
| 3428 | float p = (observed-expected)/(observed-expected+2/ploidy_*(1-observed)); |
| 3429 | if (p>-0.5 && p<1.5) { |
| 3430 | values.push_back(p); |
| 3431 | weights.push_back(chrCopyNumber_[index].getFragmentLengths_notNA_At(i)); |
| 3432 | } |
| 3433 | } |
| 3434 | } |
| 3435 | } |
| 3436 | } |
| 3437 | contam = get_median(values); |
| 3438 | contam = get_mean(values); |
| 3439 | contam = get_weighted_mean(values,weights); |
| 3440 | |
| 3441 | if (contam<0) { |
| 3442 | cout << "\t..Evaluation of contamination using contamination = (observedLevel-expectedLevel)/(1-expectedLevel) produced a negative value: "<<contam<<"\n"; |
| 3443 | cout << "\t..This probably means that there is not contamination at all or your \"ploidy\" value is not correct\n Will continue with contamination = 0\n"; |
| 3444 | contam=0; |
| 3445 | } |
| 3446 | |
| 3447 | |
| 3448 | values.clear(); |
| 3449 | return contam; |
| 3450 | } |
| 3451 | |
| 3452 | float GenomeCopyNumber::evaluateContaminationwithLR () { |
no test coverage detected