| 1024 | } |
| 1025 | |
| 1026 | double GenomeCopyNumber::calculateMedianAround (float interval, float around) { |
| 1027 | |
| 1028 | float maxCG = around+interval; |
| 1029 | float minCG = around-interval; |
| 1030 | |
| 1031 | vector <float> myValuesAround; |
| 1032 | vector<ChrCopyNumber>::iterator it; |
| 1033 | |
| 1034 | |
| 1035 | int countAutosomes = 0; |
| 1036 | for ( it=chrCopyNumber_.begin() ; it != chrCopyNumber_.end(); it++ ) { |
| 1037 | if (! (it->getChromosome().find("X")!=string::npos || it->getChromosome().find("Y")!=string::npos)) { |
| 1038 | countAutosomes++; |
| 1039 | } |
| 1040 | } |
| 1041 | |
| 1042 | |
| 1043 | for ( it=chrCopyNumber_.begin() ; it != chrCopyNumber_.end(); it++ ) { |
| 1044 | if (countAutosomes>0 && ! (it->getChromosome().find("X")!=string::npos || it->getChromosome().find("Y")!=string::npos) || countAutosomes==0) |
| 1045 | for (int i = 0; i< it->getLength(); i++) { |
| 1046 | if ((it->getCGprofileAt(i)<=maxCG)&&(it->getCGprofileAt(i)>=minCG)) |
| 1047 | if (it->getValueAt(i) >0) //non-zero values |
| 1048 | myValuesAround.push_back(it->getValueAt(i)); |
| 1049 | } |
| 1050 | } |
| 1051 | if (myValuesAround.size()==0) { |
| 1052 | cerr << "Error: zero reads in windows with the GC-content around "<<around<< " with interval "<< interval<<", will try again with "<< interval*4<<"\n"; |
| 1053 | interval=interval*4; |
| 1054 | maxCG = around+interval; |
| 1055 | minCG = around-interval; |
| 1056 | |
| 1057 | for ( it=chrCopyNumber_.begin() ; it != chrCopyNumber_.end(); it++ ) { |
| 1058 | if (! (it->getChromosome().find("X")!=string::npos || it->getChromosome().find("Y")!=string::npos)) |
| 1059 | for (int i = 0; i< it->getLength(); i++) { |
| 1060 | if ((it->getCGprofileAt(i)<=maxCG)&&(it->getCGprofileAt(i)>=minCG)) |
| 1061 | if (it->getValueAt(i) >0) //non-zero values |
| 1062 | myValuesAround.push_back(it->getValueAt(i)); |
| 1063 | } |
| 1064 | } |
| 1065 | if (myValuesAround.size()==0) { |
| 1066 | cerr << "Error: zero reads in windows with the GC-content around "<<around<< " with interval "<< interval<<"\n"; |
| 1067 | cerr << "Unable to proceed..\n"; |
| 1068 | cerr << "Try to rerun the program with higher number of reads\n"; |
| 1069 | exit(-1); |
| 1070 | } |
| 1071 | } |
| 1072 | float median = get_median(myValuesAround); |
| 1073 | myValuesAround.clear(); |
| 1074 | return median; |
| 1075 | } |
| 1076 | |
| 1077 | double GenomeCopyNumber::calculateMedianAround (GenomeCopyNumber & controlCopyNumber, float interval, float around) { |
| 1078 |
nothing calls this directly
no test coverage detected