From ea21ba73273891261ba2e4d0d85729f308c54d72 Mon Sep 17 00:00:00 2001 From: Pjotr Prins Date: Sat, 26 Aug 2017 07:41:48 +0000 Subject: Tests and enforces added related to https://github.com/genetics-statistics/GEMMA/issues/78 --- src/io.cpp | 44 ++++++++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 14 deletions(-) (limited to 'src/io.cpp') diff --git a/src/io.cpp b/src/io.cpp index fedf69a..f5783cd 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -133,7 +133,7 @@ std::istream &safeGetline(std::istream &is, std::string &t) { } } -// Read SNP file. +// Read SNP file. A single column of SNP names. bool ReadFile_snps(const string file_snps, set &setSnps) { setSnps.clear(); @@ -148,6 +148,7 @@ bool ReadFile_snps(const string file_snps, set &setSnps) { while (getline(infile, line)) { ch_ptr = strtok((char *)line.c_str(), " , \t"); + enforce_msg(ch_ptr,"Problem reading SNP file"); setSnps.insert(ch_ptr); } @@ -157,6 +158,9 @@ bool ReadFile_snps(const string file_snps, set &setSnps) { return true; } +// Read SNP file using a header. The header determines how the +// values for each row are parsed. A valid header can be, for example, +// RS POS CHR bool ReadFile_snps_header(const string &file_snps, set &setSnps) { setSnps.clear(); @@ -175,7 +179,7 @@ bool ReadFile_snps_header(const string &file_snps, set &setSnps) { ReadHeader_io(line, header); if (header.rs_col == 0 && (header.chr_col == 0 || header.pos_col == 0)) { - cout << "missing rs id in the hearder" << endl; + cout << "missing rs id in the header" << endl; } while (!safeGetline(infile, line).eof()) { @@ -185,6 +189,7 @@ bool ReadFile_snps_header(const string &file_snps, set &setSnps) { ch_ptr = strtok((char *)line.c_str(), " , \t"); for (size_t i = 0; i < header.coln; i++) { + enforce_msg(ch_ptr,"Problem reading SNP file"); if (header.rs_col != 0 && header.rs_col == i + 1) { rs = ch_ptr; } @@ -250,7 +255,7 @@ bool ReadFile_log(const string &file_log, double &pheno_mean) { return true; } -// Read bimbam annotation file. +// Read bimbam annotation file which consists of rows of SNP, POS and CHR bool ReadFile_anno(const string &file_anno, map &mapRS2chr, map &mapRS2bp, map &mapRS2cM) { @@ -264,33 +269,39 @@ bool ReadFile_anno(const string &file_anno, map &mapRS2chr, } string line; - char *ch_ptr; - - string rs; - long int b_pos; - string chr; - double cM; while (!safeGetline(infile, line).eof()) { - ch_ptr = strtok((char *)line.c_str(), " , \t"); - rs = ch_ptr; + const char *ch_ptr = strtok((char *)line.c_str(), " , \t"); + enforce_str(ch_ptr, line + " Bad RS format"); + const string rs = ch_ptr; + enforce_str(rs != "", line + " Bad RS format"); + ch_ptr = strtok(NULL, " , \t"); + enforce_str(ch_ptr, line + " Bad format"); + ulong b_pos; if (strcmp(ch_ptr, "NA") == 0) { b_pos = -9; } else { b_pos = atol(ch_ptr); } + enforce_str(b_pos,line + " Bad pos format (is zero)"); + + string chr; ch_ptr = strtok(NULL, " , \t"); if (ch_ptr == NULL || strcmp(ch_ptr, "NA") == 0) { chr = "-9"; } else { chr = ch_ptr; + enforce_str(chr != "", line + " Bad chr format"); } + + double cM; ch_ptr = strtok(NULL, " , \t"); if (ch_ptr == NULL || strcmp(ch_ptr, "NA") == 0) { cM = -9; } else { cM = atof(ch_ptr); + enforce_str(b_pos, line + "Bad cM format (is zero)"); } mapRS2chr[rs] = chr; @@ -377,9 +388,9 @@ bool ReadFile_pheno(const string &file_pheno, while (!safeGetline(infile, line).eof()) { ch_ptr = strtok((char *)line.c_str(), " , \t"); - size_t i = 0; while (i < p_max) { + enforce_msg(ch_ptr,"Wrong number of phenotypes"); if (mapP2c.count(i + 1) != 0) { if (strcmp(ch_ptr, "NA") == 0) { ind_pheno_row[mapP2c[i + 1]] = 0; @@ -618,7 +629,8 @@ bool ReadFile_geno(const string &file_geno, const set &setSnps, gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); int sig; LUDecomp(WtW, pmt, &sig); - LUInvert(WtW, pmt, WtWi); + + LUInvert(WtW, pmt, WtWi); // @@ double v_x, v_w; int c_idv = 0; @@ -667,6 +679,10 @@ bool ReadFile_geno(const string &file_geno, const set &setSnps, } if (mapRS2bp.count(rs) == 0) { + std::string msg = "Can't figure out position for "; + msg += rs; + debug_msg(msg); + chr = "-9"; b_pos = -9; cM = -9; @@ -2979,7 +2995,7 @@ bool bgenKin(const string &file_oxford, vector &indicator_snp, bool ReadHeader_io(const string &line, HEADER &header) { string rs_ptr[] = {"rs", "RS", "snp", "SNP", "snps", "SNPS", "snpid", "SNPID", "rsid", "RSID", "MarkerName"}; - set rs_set(rs_ptr, rs_ptr + 11); + set rs_set(rs_ptr, rs_ptr + 11); // create a set of 11 items string chr_ptr[] = {"chr", "CHR"}; set chr_set(chr_ptr, chr_ptr + 2); string pos_ptr[] = { -- cgit v1.2.3 From be45bcca0ddc1d88066bd18a56ca57d588e41e99 Mon Sep 17 00:00:00 2001 From: Pjotr Prins Date: Sat, 26 Aug 2017 10:00:18 +0000 Subject: Debug: only output check on -debug --- src/io.cpp | 15 ++++++++------- src/io.h | 6 +++--- src/param.cpp | 8 ++++---- test/dev_test_suite.sh | 4 +++- 4 files changed, 18 insertions(+), 15 deletions(-) (limited to 'src/io.cpp') diff --git a/src/io.cpp b/src/io.cpp index f5783cd..80adbe6 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -608,7 +608,7 @@ bool ReadFile_geno(const string &file_geno, const set &setSnps, const double &r2_level, map &mapRS2chr, map &mapRS2bp, map &mapRS2cM, vector &snpInfo, - size_t &ns_test) { + size_t &ns_test, bool debug) { indicator_snp.clear(); snpInfo.clear(); @@ -679,10 +679,11 @@ bool ReadFile_geno(const string &file_geno, const set &setSnps, } if (mapRS2bp.count(rs) == 0) { - std::string msg = "Can't figure out position for "; - msg += rs; - debug_msg(msg); - + if (debug) { + std::string msg = "Can't figure out position for "; + msg += rs; + debug_msg(msg); + } chr = "-9"; b_pos = -9; cM = -9; @@ -1633,7 +1634,7 @@ bool PlinkKin(const string &file_bed, vector &indicator_snp, // genotype and calculate K. bool ReadFile_geno(const string file_geno, vector &indicator_idv, vector &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, - const bool calc_K) { + const bool calc_K, bool debug) { igzstream infile(file_geno.c_str(), igzstream::in); if (!infile) { cout << "error reading genotype file:" << file_geno << endl; @@ -1737,7 +1738,7 @@ bool ReadFile_geno(const string &file_geno, vector &indicator_idv, vector &indicator_snp, vector> &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, - const size_t ns_test) { + const size_t ns_test, bool debug) { igzstream infile(file_geno.c_str(), igzstream::in); if (!infile) { cout << "error reading genotype file:" << file_geno << endl; diff --git a/src/io.h b/src/io.h index 9edc5eb..d9253e3 100644 --- a/src/io.h +++ b/src/io.h @@ -64,7 +64,7 @@ bool ReadFile_geno(const string &file_geno, const set &setSnps, const double &r2_level, map &mapRS2chr, map &mapRS2bp, map &mapRS2cM, vector &snpInfo, - size_t &ns_test); + size_t &ns_test, bool debug); bool ReadFile_bed(const string &file_bed, const set &setSnps, const gsl_matrix *W, vector &indicator_idv, vector &indicator_snp, vector &snpInfo, @@ -94,7 +94,7 @@ bool PlinkKin(const string &file_bed, vector &indicator_snp, bool ReadFile_geno(const string file_geno, vector &indicator_idv, vector &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, - const bool calc_K); + const bool calc_K, bool debug); bool ReadFile_bed(const string &file_bed, vector &indicator_idv, vector &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, const bool calc_K); @@ -102,7 +102,7 @@ bool ReadFile_geno(const string &file_geno, vector &indicator_idv, vector &indicator_snp, vector> &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, - const size_t ns_test); + const size_t ns_test, bool debug); bool ReadFile_bed(const string &file_bed, vector &indicator_idv, vector &indicator_snp, vector> &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, diff --git a/src/param.cpp b/src/param.cpp index c81eaf5..dce55e0 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -337,7 +337,7 @@ void PARAM::ReadFiles(void) { trim_individuals(indicator_cvt, ni_max, mode_debug); if (ReadFile_geno(file_geno, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, - mapRS2bp, mapRS2cM, snpInfo, ns_test) == false) { + mapRS2bp, mapRS2cM, snpInfo, ns_test, mode_debug) == false) { error = true; } gsl_matrix_free(W); @@ -445,7 +445,7 @@ void PARAM::ReadFiles(void) { while (!safeGetline(infile, file_name).eof()) { if (ReadFile_geno(file_name, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, - mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp) == false) { + mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp, mode_debug) == false) { error = true; } @@ -1338,7 +1338,7 @@ void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) { } } else { if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, UtX, K, - calc_K) == false) { + calc_K, mode_debug) == false) { error = true; } } @@ -1358,7 +1358,7 @@ void PARAM::ReadGenotypes(vector> &Xt, gsl_matrix *K, } } else { if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, Xt, K, calc_K, - ni_test, ns_test) == false) { + ni_test, ns_test, mode_debug) == false) { error = true; } } diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh index 67d9cba..9e49251 100755 --- a/test/dev_test_suite.sh +++ b/test/dev_test_suite.sh @@ -24,7 +24,8 @@ testBXDStandardRelatednessMatrixK() { -c ../example/BXD_covariates2.txt \ -a ../example/BXD_snps.txt \ -gk \ - -debug -o $outn + -debug \ + -o $outn assertEquals 0 $? outfn=output/$outn.cXX.txt assertEquals "198" `wc -l < $outfn` @@ -38,6 +39,7 @@ testBXDMultivariateLinearMixedModel() { -a ../example/BXD_snps.txt \ -k ./output/BXD.cXX.txt \ -lmm 2 -maf 0.1 \ + -debug \ -o BXD_mvlmm assertEquals 0 $? -- cgit v1.2.3