test/vcf.cc

Code
Comments
Other
Rev Date Author Line
3747 14 Aug 18 peter 1 // $Id$
3747 14 Aug 18 peter 2 //
4316 17 Feb 23 peter 3 // Copyright (C) 2018, 2020, 2022, 2023 Peter Johansson
3747 14 Aug 18 peter 4 //
3747 14 Aug 18 peter 5 // This program is free software; you can redistribute it and/or modify
3747 14 Aug 18 peter 6 // it under the terms of the GNU General Public License as published by
3747 14 Aug 18 peter 7 // the Free Software Foundation; either version 3 of the License, or
3747 14 Aug 18 peter 8 // (at your option) any later version.
3747 14 Aug 18 peter 9 //
3747 14 Aug 18 peter 10 // This program is distributed in the hope that it will be useful, but
3747 14 Aug 18 peter 11 // WITHOUT ANY WARRANTY; without even the implied warranty of
3747 14 Aug 18 peter 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
3747 14 Aug 18 peter 13 // General Public License for more details.
3747 14 Aug 18 peter 14 //
3747 14 Aug 18 peter 15 // You should have received a copy of the GNU General Public License
3747 14 Aug 18 peter 16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
3747 14 Aug 18 peter 17
3747 14 Aug 18 peter 18
3747 14 Aug 18 peter 19 // test class omic::VCF
3747 14 Aug 18 peter 20
3747 14 Aug 18 peter 21 #include <config.h>
3747 14 Aug 18 peter 22
3747 14 Aug 18 peter 23 #include "Suite.h"
3747 14 Aug 18 peter 24
3747 14 Aug 18 peter 25 #include "yat/omic/VCF.h"
4145 24 Feb 22 peter 26 #include "yat/utility/split.h"
4021 26 Nov 20 peter 27 #include "yat/utility/utility.h"
3747 14 Aug 18 peter 28
3747 14 Aug 18 peter 29 #include <fstream>
3747 14 Aug 18 peter 30 #include <string>
3747 14 Aug 18 peter 31
3747 14 Aug 18 peter 32 using namespace theplu::yat;
3747 14 Aug 18 peter 33
4316 17 Feb 23 peter 34 template<typename T>
4316 17 Feb 23 peter 35 void equal_vec(test::Suite& suite, std::vector<T> x, std::vector<T> y)
4316 17 Feb 23 peter 36 {
4316 17 Feb 23 peter 37   if (x != y) {
4316 17 Feb 23 peter 38     suite.add(false);
4316 17 Feb 23 peter 39     if (x.size() != y.size())
4316 17 Feb 23 peter 40       suite.err() << "error: vectors not the same\n"
4316 17 Feb 23 peter 41                   << "x size: " << x.size() << "\n"
4316 17 Feb 23 peter 42                   << "y size: " << y.size() << "\n";
4316 17 Feb 23 peter 43     else
4316 17 Feb 23 peter 44       for (size_t i=0; i<x.size(); ++i)
4316 17 Feb 23 peter 45         std::cerr << "'" << x[i] << "' : '" << y[i] << "'\n";
4316 17 Feb 23 peter 46   }
4316 17 Feb 23 peter 47 }
4316 17 Feb 23 peter 48
4316 17 Feb 23 peter 49
4316 17 Feb 23 peter 50 template<typename T1, typename T2>
4316 17 Feb 23 peter 51 void equal(test::Suite& suite, T1 x, T2 y)
4316 17 Feb 23 peter 52 {
4316 17 Feb 23 peter 53   if (x != y) {
4316 17 Feb 23 peter 54     suite.add(false);
4316 17 Feb 23 peter 55     suite.err() << "error: " << x << " != " << y << "\n";
4316 17 Feb 23 peter 56   }
4316 17 Feb 23 peter 57 }
4316 17 Feb 23 peter 58
4316 17 Feb 23 peter 59
4316 17 Feb 23 peter 60 void test_info(test::Suite& suite, omic::VCF::Info info)
4316 17 Feb 23 peter 61 {
4316 17 Feb 23 peter 62   std::string value;
4316 17 Feb 23 peter 63   info.add("FOO", "foo");
4316 17 Feb 23 peter 64   info.get(value, "FOO");
4316 17 Feb 23 peter 65   equal(suite, value, "foo");
4316 17 Feb 23 peter 66
4316 17 Feb 23 peter 67   info.set("FOO", "foo2");
4316 17 Feb 23 peter 68   info.get(value, "FOO");
4316 17 Feb 23 peter 69   equal(suite, value, "foo2");
4316 17 Feb 23 peter 70
4316 17 Feb 23 peter 71   int x;
4316 17 Feb 23 peter 72   info.add("X", "1");
4316 17 Feb 23 peter 73   info.get(x, "X");
4316 17 Feb 23 peter 74   equal(suite, x, 1);
4316 17 Feb 23 peter 75
4316 17 Feb 23 peter 76   info.set("X", 2);
4316 17 Feb 23 peter 77   info.get(x, "X");
4316 17 Feb 23 peter 78   equal(suite, x, 2);
4316 17 Feb 23 peter 79
4316 17 Feb 23 peter 80   std::vector<std::string> values;
4316 17 Feb 23 peter 81   std::vector<std::string> foos(2, "fu");
4316 17 Feb 23 peter 82   info.add("FOOS", foos);
4316 17 Feb 23 peter 83   info.get(values, "FOOS");
4316 17 Feb 23 peter 84   equal_vec(suite, values, foos);
4316 17 Feb 23 peter 85   info.set("FOOS", foos);
4316 17 Feb 23 peter 86   info.get(value, "FOOS");
4316 17 Feb 23 peter 87   equal_vec(suite, values, foos);
4316 17 Feb 23 peter 88
4316 17 Feb 23 peter 89   std::vector<int> y;
4316 17 Feb 23 peter 90   std::vector<int> bars(2, 13);
4316 17 Feb 23 peter 91   info.add("BARS", bars);
4316 17 Feb 23 peter 92   info.get(y, "BARS");
4316 17 Feb 23 peter 93   equal_vec(suite, y, bars);
4316 17 Feb 23 peter 94   info.set("BARS", bars);
4316 17 Feb 23 peter 95   info.get(y, "BARS");
4316 17 Feb 23 peter 96   equal_vec(suite, y, bars);
4316 17 Feb 23 peter 97 }
4316 17 Feb 23 peter 98
4316 17 Feb 23 peter 99
3747 14 Aug 18 peter 100 int main(int argc, char* argv[])
3747 14 Aug 18 peter 101 {
3747 14 Aug 18 peter 102   test::Suite suite(argc, argv);
3747 14 Aug 18 peter 103   try {
3747 14 Aug 18 peter 104     using omic::VCF;
3747 14 Aug 18 peter 105     if (false) {
3747 14 Aug 18 peter 106       VCF vcf;
3747 14 Aug 18 peter 107       test::avoid_compiler_warning(vcf);
3747 14 Aug 18 peter 108     }
3747 14 Aug 18 peter 109     if (false) {
3747 14 Aug 18 peter 110       std::ifstream is("nonsense.vcf");
4021 26 Nov 20 peter 111       VCF vcf(is);
3747 14 Aug 18 peter 112       test::avoid_compiler_warning(vcf);
3747 14 Aug 18 peter 113     }
3747 14 Aug 18 peter 114     VCF vcf("1\t124\t.\tC\tG\t60\tPASS\t\tGT\t0/1\t0/0\t0/0");
3747 14 Aug 18 peter 115
3747 14 Aug 18 peter 116     suite.out() << "test chromosome\n";
3747 14 Aug 18 peter 117     suite.add(vcf.chromosome() == "1");
3747 14 Aug 18 peter 118     vcf.chromosome("chr1");
3747 14 Aug 18 peter 119     suite.add(vcf.chromosome() == "chr1");
3747 14 Aug 18 peter 120
3747 14 Aug 18 peter 121     suite.out() << "test pos\n";
3747 14 Aug 18 peter 122     suite.add(vcf.pos() == 124);
3747 14 Aug 18 peter 123     vcf.pos(52960);
3747 14 Aug 18 peter 124     suite.add(vcf.pos() == 52960);
3747 14 Aug 18 peter 125
3747 14 Aug 18 peter 126     suite.out() << "test id\n";
3747 14 Aug 18 peter 127     suite.add(vcf.id() == ".");
3747 14 Aug 18 peter 128     vcf.id("rs1");
3747 14 Aug 18 peter 129     suite.add(vcf.id() == "rs1");
3747 14 Aug 18 peter 130
3747 14 Aug 18 peter 131     suite.out() << "test ref\n";
3747 14 Aug 18 peter 132     suite.add(vcf.ref() == "C");
3747 14 Aug 18 peter 133     vcf.ref('A');
3747 14 Aug 18 peter 134     vcf.ref("A");
3747 14 Aug 18 peter 135     suite.add(vcf.ref() == "A");
3747 14 Aug 18 peter 136
3747 14 Aug 18 peter 137     suite.out() << "test alt\n";
3747 14 Aug 18 peter 138     suite.add(vcf.alt() == "G");
3747 14 Aug 18 peter 139     vcf.alt('T');
3747 14 Aug 18 peter 140     vcf.alt("T");
3747 14 Aug 18 peter 141     suite.add(vcf.alt() == "T");
3747 14 Aug 18 peter 142     suite.add(vcf.alts().size() == vcf.n_alleles());
3747 14 Aug 18 peter 143     suite.add(vcf.alts()[0] == "T");
3747 14 Aug 18 peter 144
3747 14 Aug 18 peter 145     suite.out() << "test qual\n";
3747 14 Aug 18 peter 146     suite.add(vcf.qual() == "60");
3747 14 Aug 18 peter 147     vcf.qual("89");
3747 14 Aug 18 peter 148     vcf.qual(90);
3747 14 Aug 18 peter 149     suite.add(vcf.qual() == "90");
3747 14 Aug 18 peter 150
3747 14 Aug 18 peter 151     suite.out() << "test filter\n";
3747 14 Aug 18 peter 152     if (!suite.add(vcf.filter() == "PASS"))
3747 14 Aug 18 peter 153       suite.err() << "error: filter is not PASS\n";
3747 14 Aug 18 peter 154     if (!suite.add(vcf.filters().size() == 1))
3747 14 Aug 18 peter 155       suite.err() << "error: filter size is not 1\n";
3747 14 Aug 18 peter 156     vcf.filter("foo");
3747 14 Aug 18 peter 157     if (!suite.add(vcf.filter() == "foo"))
3747 14 Aug 18 peter 158       suite.err() << "error: filter not equal to 'foo'\n";
3747 14 Aug 18 peter 159     vcf.add_filter("bar");
3747 14 Aug 18 peter 160     if (!suite.add(vcf.filter() == "foo;bar"))
3747 14 Aug 18 peter 161       suite.err() << "error: filter not equal to 'foo;bar'\n";
3747 14 Aug 18 peter 162     if (suite.add(vcf.filters().size() == 2)) {
3747 14 Aug 18 peter 163       if (vcf.filters()[0] != "foo")
3747 14 Aug 18 peter 164         suite.add(false);
3747 14 Aug 18 peter 165       if (vcf.filters()[1] != "bar")
3747 14 Aug 18 peter 166         suite.add(false);
3747 14 Aug 18 peter 167     }
3747 14 Aug 18 peter 168     else
3747 14 Aug 18 peter 169       suite.err() << "error: incorrect filters size\n";
3747 14 Aug 18 peter 170
3747 14 Aug 18 peter 171     suite.out() << "test info\n";
4316 17 Feb 23 peter 172     test_info(suite, vcf.info());
3747 14 Aug 18 peter 173     suite.add(vcf.info().str() == "");
3747 14 Aug 18 peter 174     vcf.info().set("SNV");
3747 14 Aug 18 peter 175     suite.add(vcf.info().str() == "SNV");
3747 14 Aug 18 peter 176     vcf.info().add("banana");
3747 14 Aug 18 peter 177     vcf.info().add("FOO", "foo");
3747 14 Aug 18 peter 178     vcf.info().add("BAR", 3.14);
3747 14 Aug 18 peter 179     std::vector<std::string> foos(2, "fu");
3747 14 Aug 18 peter 180     vcf.info().add("FOOS", foos);
3747 14 Aug 18 peter 181     std::vector<int> bars(3, 6);
3747 14 Aug 18 peter 182     vcf.info().add("BARS", bars);
3747 14 Aug 18 peter 183     suite.out() << "info: " << vcf.info() << "\n";
3747 14 Aug 18 peter 184     std::string correct = "SNV;banana;FOO=foo;BAR=3.14;FOOS=fu,fu;BARS=6,6,6";
3747 14 Aug 18 peter 185     if (vcf.info().str() != correct) {
3747 14 Aug 18 peter 186       suite.add(false);
3747 14 Aug 18 peter 187       suite.err() << "error: incorrect info\n";
3747 14 Aug 18 peter 188       suite.err() << "info:     " << vcf.info() << "\n";
3747 14 Aug 18 peter 189       suite.err() << "expected: " << correct << "\n";
3747 14 Aug 18 peter 190     }
4164 13 Mar 22 peter 191     if (vcf.info().count("BAR") != 1) {
4164 13 Mar 22 peter 192       suite.add(false);
4164 13 Mar 22 peter 193       suite.err() << "error: count(\"BAR\") not returning 1\n";
4164 13 Mar 22 peter 194     }
3747 14 Aug 18 peter 195
4021 26 Nov 20 peter 196     // test getting int when value is a String
4021 26 Nov 20 peter 197     try {
4021 26 Nov 20 peter 198       double x;
4021 26 Nov 20 peter 199       vcf.info().get(x, "FOO");
4021 26 Nov 20 peter 200       suite.add(false);
4021 26 Nov 20 peter 201     }
4021 26 Nov 20 peter 202     catch (utility::runtime_error& e) {
4021 26 Nov 20 peter 203       std::ostringstream os;
4021 26 Nov 20 peter 204       utility::print_what(e, os);
4021 26 Nov 20 peter 205       suite.out() << "expected exception: what: " << os.str();
4021 26 Nov 20 peter 206       if (os.str().find("FOO")==std::string::npos) {
4021 26 Nov 20 peter 207         suite.add(false);
4021 26 Nov 20 peter 208         suite.err() << "error: expected 'FOO' mentioned in what()\n";
4021 26 Nov 20 peter 209       }
4021 26 Nov 20 peter 210     }
4021 26 Nov 20 peter 211
3747 14 Aug 18 peter 212     suite.out() << "test data\n";
3747 14 Aug 18 peter 213     std::vector<std::string> GT;
3747 14 Aug 18 peter 214     std::vector<std::string> GT_correct(3, "0/0");
3747 14 Aug 18 peter 215     GT_correct[0] = "0/1";
3747 14 Aug 18 peter 216     vcf.data().get("GT", GT);
3747 14 Aug 18 peter 217     if (GT != GT_correct) {
3747 14 Aug 18 peter 218       suite.err() << "error: incorrect GT\n";
3747 14 Aug 18 peter 219       suite.add(false);
3747 14 Aug 18 peter 220     }
3747 14 Aug 18 peter 221
3747 14 Aug 18 peter 222     std::vector<int> DP;
3747 14 Aug 18 peter 223     DP.push_back(34);
3747 14 Aug 18 peter 224     DP.push_back(87);
3747 14 Aug 18 peter 225     DP.push_back(63);
3747 14 Aug 18 peter 226     vcf.data().add("DP", DP);
3747 14 Aug 18 peter 227
3747 14 Aug 18 peter 228     std::vector<std::vector<int> > PL(3, std::vector<int>(3, 255));
3747 14 Aug 18 peter 229     PL[0][1] = 0;
3747 14 Aug 18 peter 230     PL[1][0] = 0;
3747 14 Aug 18 peter 231     PL[2][0] = 0;
3747 14 Aug 18 peter 232     vcf.data().add("PL", PL);
4170 09 May 22 peter 233     vcf.data().set("PL", PL);
4170 09 May 22 peter 234     vcf.data().set("PL2", PL);
3747 14 Aug 18 peter 235
3747 14 Aug 18 peter 236     std::vector<std::string> ID;
3747 14 Aug 18 peter 237     ID.push_back("foo");
3747 14 Aug 18 peter 238     ID.push_back("bar");
3747 14 Aug 18 peter 239     ID.push_back("baz");
3747 14 Aug 18 peter 240     vcf.data().add("ID", ID);
4021 26 Nov 20 peter 241
4021 26 Nov 20 peter 242     // test unknown key
4021 26 Nov 20 peter 243     try {
4021 26 Nov 20 peter 244       vcf.data().get("banana", GT);
4021 26 Nov 20 peter 245     }
4021 26 Nov 20 peter 246     catch (utility::runtime_error& e) {
4021 26 Nov 20 peter 247       std::ostringstream os;
4021 26 Nov 20 peter 248       utility::print_what(e, os);
4021 26 Nov 20 peter 249       suite.out() << "expected exception: what: " << os.str();
4021 26 Nov 20 peter 250       std::vector<std::string> tags = { "unknown key", "banana"};
4021 26 Nov 20 peter 251       for (size_t i=0; i<tags.size(); ++i)
4021 26 Nov 20 peter 252         if (os.str().find(tags[i]) == std::string::npos) {
4021 26 Nov 20 peter 253           suite.add(false);
4021 26 Nov 20 peter 254           suite.err() << "error: expected '" << tags[i]
4021 26 Nov 20 peter 255                       << "' in error message\n";
4021 26 Nov 20 peter 256         }
4021 26 Nov 20 peter 257     }
4021 26 Nov 20 peter 258
4021 26 Nov 20 peter 259     // test getting int when value is a String
4021 26 Nov 20 peter 260     try {
4021 26 Nov 20 peter 261       std::vector<int> p;
4021 26 Nov 20 peter 262       vcf.data().get("ID", p);
4021 26 Nov 20 peter 263     }
4021 26 Nov 20 peter 264     catch (utility::runtime_error& e) {
4021 26 Nov 20 peter 265       std::ostringstream os;
4021 26 Nov 20 peter 266       utility::print_what(e, os);
4021 26 Nov 20 peter 267       suite.out() << "expected exception: what: " << os.str();
4021 26 Nov 20 peter 268       if (os.str().find("ID")==std::string::npos) {
4021 26 Nov 20 peter 269         suite.add(false);
4021 26 Nov 20 peter 270         suite.err() << "error: expected 'ID' mentioned in what()\n";
4021 26 Nov 20 peter 271       }
4021 26 Nov 20 peter 272     }
4145 24 Feb 22 peter 273
4145 24 Feb 22 peter 274     {
4145 24 Feb 22 peter 275       VCF vcf2(vcf);
4145 24 Feb 22 peter 276       std::vector<size_t> index = {0, 2};
4154 09 Mar 22 peter 277       vcf2.subset(index);
4154 09 Mar 22 peter 278       suite.out() << "vcf: " << vcf << "\n";
4154 09 Mar 22 peter 279       suite.out() << "vcf2:" << vcf2 << "\n";
4145 24 Feb 22 peter 280       std::ostringstream ss;
4145 24 Feb 22 peter 281       ss << vcf;
4145 24 Feb 22 peter 282       std::vector<std::string> vec;
4145 24 Feb 22 peter 283       utility::split(vec, ss.str(), '\t');
4145 24 Feb 22 peter 284       std::ostringstream ss2;
4145 24 Feb 22 peter 285       ss2 << vcf2;
4145 24 Feb 22 peter 286       std::vector<std::string> vec2;
4145 24 Feb 22 peter 287       utility::split(vec2, ss2.str(), '\t');
4145 24 Feb 22 peter 288       for (size_t i=0; i<vec2.size(); ++i) {
4145 24 Feb 22 peter 289         size_t j = i < 9 ? i : 9 + index[i-9];
4154 09 Mar 22 peter 290         if (vec2[i] != vec[j]) {
4145 24 Feb 22 peter 291           suite.add(false);
4154 09 Mar 22 peter 292           suite.err() << "vcf != vcf2: " << vec2[i] << "!=" << vec[j] << "\n";
4145 24 Feb 22 peter 293         }
4145 24 Feb 22 peter 294       }
4145 24 Feb 22 peter 295     }
4145 24 Feb 22 peter 296     {
4145 24 Feb 22 peter 297       VCF vcf2(3);
4145 24 Feb 22 peter 298       if (vcf2.n_samples() != 3) {
4145 24 Feb 22 peter 299         suite.add(false);
4145 24 Feb 22 peter 300         suite.err() << "vcf samples not 3: " << vcf2.n_samples() << "\n";
4145 24 Feb 22 peter 301       }
4145 24 Feb 22 peter 302       vcf2.data() = VCF::Data(4);
4145 24 Feb 22 peter 303       if (vcf2.n_samples() != 4) {
4145 24 Feb 22 peter 304         suite.add(false);
4145 24 Feb 22 peter 305         suite.err() << "vcf samples not 4: " << vcf2.n_samples() << "\n";
4145 24 Feb 22 peter 306       }
4145 24 Feb 22 peter 307     }
3747 14 Aug 18 peter 308   }
3747 14 Aug 18 peter 309   catch (std::exception& e) {
3747 14 Aug 18 peter 310     suite.err() << "error: " << e.what() << "\n";
3747 14 Aug 18 peter 311     suite.add(false);
3747 14 Aug 18 peter 312   }
3747 14 Aug 18 peter 313   return suite.return_value();
3747 14 Aug 18 peter 314 }

Generated on Sat Oct 14 02:21:42 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)