| 1419 | } |
| 1420 | |
| 1421 | void ChrCopyNumber::fillCGprofile(std::string const& chrFolder) { |
| 1422 | GCprofile_ = vector <float> (length_); |
| 1423 | notNprofile_ = vector <float> (length_); |
| 1424 | ifstream file; |
| 1425 | string filename = chromosome_; |
| 1426 | string possibleFilenames[] = {filename,filename+".fa",filename+".fasta","chr"+filename+".fa","chr"+filename+".fasta"}; |
| 1427 | for (int i = 0; i < 5; i++) { |
| 1428 | |
| 1429 | string myFilename = possibleFilenames[i]; |
| 1430 | string myFullPath = pathAppend(chrFolder,myFilename); |
| 1431 | file.open(myFullPath.c_str()); |
| 1432 | |
| 1433 | if(!file.is_open()) |
| 1434 | file.clear(); |
| 1435 | else |
| 1436 | i = 6; |
| 1437 | } |
| 1438 | |
| 1439 | |
| 1440 | if (!file.is_open()) { |
| 1441 | // throw ("Unable to open fasta file for chr "+chromosome_+" in "+chrFolder+"\n"); |
| 1442 | cerr << "Unable to open fasta file for chr "+chromosome_+" in folder "+chrFolder+"\n\nPlease note, "<< chrFolder << " should be a folder, not a file!\n\n"; |
| 1443 | exit (-1); |
| 1444 | } |
| 1445 | //string myString; |
| 1446 | char letter; // here we read the chromosome letter by letter. Can do it better! |
| 1447 | file >> letter; |
| 1448 | int count = 0; |
| 1449 | int countCG = 0; |
| 1450 | int countN = 0; |
| 1451 | string line; |
| 1452 | string text = ""; |
| 1453 | |
| 1454 | if (letter == '>') |
| 1455 | getline (file,line); |
| 1456 | else { |
| 1457 | count = 1; |
| 1458 | countCG = isCG(letter); |
| 1459 | countN = isN(letter); |
| 1460 | text.push_back(letter); |
| 1461 | } |
| 1462 | |
| 1463 | if (ends_.size()==0) { //all windows have equal length => can use the same windowsize for all windows |
| 1464 | for (int i = 0; i<length_; i++) { |
| 1465 | if (file.eof()) { |
| 1466 | GCprofile_[i] = NA; |
| 1467 | //cout << "End-of-file reached.." << endl; |
| 1468 | } |
| 1469 | while((!file.eof()) && (count < windowSize_)) { |
| 1470 | file>>letter; |
| 1471 | countCG += isCG(letter); |
| 1472 | countN += isN(letter); |
| 1473 | count ++; |
| 1474 | } |
| 1475 | notNprofile_[i] = float(count-countN)/count; |
| 1476 | if (count == countN) |
| 1477 | GCprofile_[i] = NA; |
| 1478 | else |
no test coverage detected