| 1136 | } |
| 1137 | |
| 1138 | long double GenomeCopyNumber::calculateRSS(int ploidy) |
| 1139 | { |
| 1140 | string::size_type pos = 0; |
| 1141 | vector<float> observedvalues; |
| 1142 | vector<float> expectedvalues; |
| 1143 | map<string,int>::iterator it; |
| 1144 | for (it=chromosomesInd_.begin() ; it != chromosomesInd_.end(); it++ ) { |
| 1145 | string chrNumber = (*it).first; |
| 1146 | if ( ( pos = chrNumber.find("chr")) != string::npos ) |
| 1147 | chrNumber.replace( pos, 3, "" ); |
| 1148 | if ( ( pos = chrNumber.find("X")) != string::npos ) //exclude X and Y from the analysis |
| 1149 | continue; |
| 1150 | if ( ( pos = chrNumber.find("Y")) != string::npos ) |
| 1151 | continue; |
| 1152 | int index = findIndex(chrNumber); |
| 1153 | int length = chrCopyNumber_[index].getLength(); |
| 1154 | for (int i = 0; i< length; i++) { |
| 1155 | float observed = chrCopyNumber_[index].getRatioAtBin(i); |
| 1156 | if (observed!=NA) { |
| 1157 | if (isRatioLogged_) { |
| 1158 | observed=pow(2,observed); |
| 1159 | } |
| 1160 | float expected = observed; |
| 1161 | if (chrCopyNumber_[index].isMedianCalculated()) { |
| 1162 | expected = chrCopyNumber_[index].getMedianProfileAtI(i); |
| 1163 | if (chrCopyNumber_[index].isSmoothed()) |
| 1164 | expected = chrCopyNumber_[index].getSmoothedProfileAtI(i); |
| 1165 | |
| 1166 | } |
| 1167 | observedvalues.push_back(observed); |
| 1168 | expectedvalues.push_back(expected); |
| 1169 | } |
| 1170 | } |
| 1171 | } |
| 1172 | |
| 1173 | long double RSS = 0; |
| 1174 | for (int i = 0; i < (int)observedvalues.size(); i++) |
| 1175 | { |
| 1176 | if ((observedvalues[i]!=NA) && (expectedvalues[i]!=NA)) |
| 1177 | { |
| 1178 | long double diff = (long double)observedvalues[i] - (long double)round(ploidy*expectedvalues[i])/ploidy; |
| 1179 | RSS = RSS + (long double)pow(diff,2); |
| 1180 | } |
| 1181 | } |
| 1182 | if (observedvalues.size()==0) { |
| 1183 | return 0; |
| 1184 | } |
| 1185 | double normRSS = (RSS/observedvalues.size()); |
| 1186 | observedvalues.clear();expectedvalues.clear(); |
| 1187 | return normRSS; |
| 1188 | } |
| 1189 | |
| 1190 | |
| 1191 | void GenomeCopyNumber::calculateRatioUsingCG( GenomeCopyNumber & controlCopyNumber) { |
no test coverage detected