Add some additional check

Author: valeu

src/main.cpp

@@ -520,6 +520,16 @@ int main(int argc, char *argv[])
     std::string chrLenFile = (std::string)cf.Value("general","chrLenFile","");
     cout << "..File with chromosome lengths:\t" << chrLenFile << "\n";
 
+    if (targetBed!="" & chrLenFile!="") {
+        //make a check that in the chromosome length file there is no chromosomes absent in the captureRegions:
+        bool chrLenFileIsOK = checkChrLen(chrLenFile,targetBed);
+        if (!chrLenFileIsOK) {
+            cerr <<"Will exit\n";
+            exit (1);
+        }
+    }
+
+
     bool isMinMappabilitySet = cf.hasValue("general","minMappabilityPerWindow");
     minMappabilityPerWindow = double(cf.Value("general","minMappabilityPerWindow",0.85));
     if (isMinMappabilitySet && isUseGC) {

src/myFunc.cpp

@@ -71,6 +71,36 @@ std::vector<std::string> &split(const std::string &s, char delim, std::vector<st
     return elems;
 }
 
+bool checkChrLen(const std::string &chrLenFile,const std::string &targetBed) {
+    std::vector<std::string> names;
+	std::vector<int> lengths;
+    std::vector<std::string> names_bed;
+	readFileWithGenomeInfo(chrLenFile, names, lengths);
+    readChrNamesInBed(targetBed, names_bed);
+    bool toReturn = true;
+
+    if(names.empty()){
+        cerr << "Error:Cound not read "<< chrLenFile<<"\n";
+        exit(1);
+    }
+    if(names_bed.empty()){
+        cerr << "Error:Cound not read "<< targetBed<<"\n";
+        exit(1);
+    }
+
+    for (int i=0; i<names.size();i++) {
+        if(std::find(names_bed.begin(), names_bed.end(), names[i]) != names_bed.end()) {
+                /* names_bed contains names[i]; everything is OK */
+        } else {
+                /* names_bed does not contain names[i] */
+                toReturn = false;
+                cerr << "Error: chromosome "<< names[i]<< " present in your "<<chrLenFile << " file was not detected in your file with capture regions " <<targetBed<<"\n";
+                cerr << "Please solve this issue and rerun Control-FREEC\n";
+                cerr << "For example, you can remove chromosome "<< names[i]<<" from your "<<chrLenFile<<"\n";
+        }
+    }
+    return toReturn;
+}
 
 std::vector<std::string> split(const std::string &s, char delim) {
     std::vector<std::string> elems;
@@ -234,6 +264,38 @@ float get_iqr(const std::vector<float>& myvector) {
   return upper_quartile - lower_quartile;
 }
 
+void readChrNamesInBed(const std::string &targetBed, std::vector<std::string>&names_bed){
+    ifstream file(targetBed.c_str());
+	if (!file.is_open()) {
+        cerr << "Error: unable to open "+targetBed+"\n" ;
+        exit(-1);
+	}
+	string line;
+	string name;
+	bool isFai=0;
+	while (std::getline(file,line)) {
+		if (! line.length()) continue;
+		if (line[0] == '#') continue;
+
+		std::vector<std::string> strs = split(line, '\t');
+		if (strs.size()<2) {
+		    continue;
+		}
+        name  = strs[0];
+		strs.clear();
+		myReplace(name, " ", "");
+
+		//delete "Chr"
+		string::size_type pos = 0;
+		if ( (pos = name.find("chr", pos)) != string::npos )
+			name.replace( pos, 3, "" );
+
+        names_bed.push_back(name);
+	}
+	file.close();
+}
+
+
 void readFileWithGenomeInfo(const std::string &chrLenFileName, std::vector<std::string>& names, std::vector<int>& lengths) {
 	//reading the file with genome information
 	ifstream file(chrLenFileName.c_str());

src/myFunc.h

@@ -83,6 +83,7 @@ float get_weighted_mean(const std::vector<float>& data, const std::vector<float>
 float get_sum(const std::vector<float>& data) ;
 float get_iqr(const std::vector<float>& data);
 void readFileWithGenomeInfo(const std::string &chrLenFileName, std::vector<std::string>& names, std::vector<int>& lengths);
+void readChrNamesInBed(const std::string &targetBed, std::vector<std::string>&names_bed);
 unsigned long sum(const std::vector<int>& data);
 long getLineNumber(std::string const& file, const std::string& pathToSamtools, const std::string& pathToSambamba, const std::string& SambambaThreads);
 long getReadNumberFromPileup(std::string const& file);
@@ -150,6 +151,7 @@ void advance_to(const std::string& haystack, size_t& offset, char needle) ;
 std::vector<float> get_quartiles(std::vector<float> vect);
 
 int calculateTotalLength(std::vector <int> lefts,std::vector <int> rights);
+bool checkChrLen(const std::string &chrLenFile,const std::string &targetBed) ;
 
 #ifdef _WIN32
 double expm1(double x);