test/bam_header.cc

Code
Comments
Other
Rev Date Author Line
2953 16 Jan 13 peter 1 // $Id$
2953 16 Jan 13 peter 2 //
4359 23 Aug 23 peter 3 // Copyright (C) 2013, 2014, 2015, 2016, 2017, 2020, 2023 Peter Johansson
2953 16 Jan 13 peter 4 //
2953 16 Jan 13 peter 5 // This program is free software; you can redistribute it and/or modify
2953 16 Jan 13 peter 6 // it under the terms of the GNU General Public License as published by
2953 16 Jan 13 peter 7 // the Free Software Foundation; either version 3 of the License, or
2953 16 Jan 13 peter 8 // (at your option) any later version.
2953 16 Jan 13 peter 9 //
2953 16 Jan 13 peter 10 // This program is distributed in the hope that it will be useful, but
2953 16 Jan 13 peter 11 // WITHOUT ANY WARRANTY; without even the implied warranty of
2953 16 Jan 13 peter 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
2953 16 Jan 13 peter 13 // General Public License for more details.
2953 16 Jan 13 peter 14 //
2953 16 Jan 13 peter 15 // You should have received a copy of the GNU General Public License
2953 16 Jan 13 peter 16 // along with this program. If not, see <http://www.gnu.org/licenses/>.
2953 16 Jan 13 peter 17
2953 16 Jan 13 peter 18 #include <config.h>
2953 16 Jan 13 peter 19
2953 16 Jan 13 peter 20 #include "Suite.h"
2953 16 Jan 13 peter 21
3883 24 Mar 20 peter 22 #ifdef YAT_HAVE_HTSLIB
2953 16 Jan 13 peter 23 #include "yat/omic/BamFile.h"
2953 16 Jan 13 peter 24 #include "yat/omic/BamHeader.h"
4290 03 Feb 23 peter 25 #include <set>
2953 16 Jan 13 peter 26 #endif
2953 16 Jan 13 peter 27
4290 03 Feb 23 peter 28
2953 16 Jan 13 peter 29 using namespace theplu::yat;
2953 16 Jan 13 peter 30
4264 13 Jan 23 peter 31 #ifdef YAT_HAVE_HTSLIB
2953 16 Jan 13 peter 32 void test1(test::Suite& suite);
4264 13 Jan 23 peter 33 #else
4264 13 Jan 23 peter 34 void test1(test::Suite& suite) {}
4264 13 Jan 23 peter 35 #endif
2953 16 Jan 13 peter 36
2953 16 Jan 13 peter 37 int main(int argc, char* argv[])
2953 16 Jan 13 peter 38 {
3201 03 May 14 peter 39   test::Suite suite(argc, argv, true);
2953 16 Jan 13 peter 40   try {
2953 16 Jan 13 peter 41     test1(suite);
2953 16 Jan 13 peter 42   }
2953 16 Jan 13 peter 43   catch (std::runtime_error& e) {
2953 16 Jan 13 peter 44     suite.err() << "what: " << e.what() << "\n";
2953 16 Jan 13 peter 45     suite.add(false);
2953 16 Jan 13 peter 46   }
2953 16 Jan 13 peter 47
2953 16 Jan 13 peter 48   return suite.return_value();
2953 16 Jan 13 peter 49 }
2953 16 Jan 13 peter 50
4264 13 Jan 23 peter 51
4264 13 Jan 23 peter 52 #ifdef YAT_HAVE_HTSLIB
4264 13 Jan 23 peter 53 using namespace omic;
4264 13 Jan 23 peter 54
4264 13 Jan 23 peter 55 void test_copy_and_move(test::Suite& suite, const BamHeader& header)
4264 13 Jan 23 peter 56 {
4264 13 Jan 23 peter 57   suite.out() << "test copy constructor\n";
4264 13 Jan 23 peter 58   {
4264 13 Jan 23 peter 59     if (true) {
4264 13 Jan 23 peter 60       BamHeader copy(header);
4264 13 Jan 23 peter 61       if (copy.text() != header.text()) {
4264 13 Jan 23 peter 62         suite.add(false);
4264 13 Jan 23 peter 63         suite.err() << "copy not identical\n";
4264 13 Jan 23 peter 64       }
4264 13 Jan 23 peter 65     }
4264 13 Jan 23 peter 66     // test that copy destructor doesn't destroy header data
4264 13 Jan 23 peter 67     std::string str = header.text();
4264 13 Jan 23 peter 68     suite.add(str.size());
4264 13 Jan 23 peter 69   }
4264 13 Jan 23 peter 70
4264 13 Jan 23 peter 71   suite.out() << "test move constructor\n";
4264 13 Jan 23 peter 72   {
4264 13 Jan 23 peter 73     if (true) {
4264 13 Jan 23 peter 74       BamHeader copy(header);
4264 13 Jan 23 peter 75       BamHeader hdr(std::move(copy));
4264 13 Jan 23 peter 76       if (hdr.text() != header.text()) {
4264 13 Jan 23 peter 77         suite.add(false);
4264 13 Jan 23 peter 78         suite.err() << "copy constructed with move ctor not identical\n";
4264 13 Jan 23 peter 79       }
4264 13 Jan 23 peter 80     }
4264 13 Jan 23 peter 81     std::string str = header.text();
4264 13 Jan 23 peter 82     suite.add(str.size());
4264 13 Jan 23 peter 83   }
4264 13 Jan 23 peter 84
4264 13 Jan 23 peter 85
4264 13 Jan 23 peter 86   suite.out() << "test copy assignment\n";
4264 13 Jan 23 peter 87   {
4264 13 Jan 23 peter 88     {
4264 13 Jan 23 peter 89       BamHeader copy;
4264 13 Jan 23 peter 90       copy = header;
4264 13 Jan 23 peter 91       if (copy.text() != header.text()) {
4264 13 Jan 23 peter 92         suite.add(false);
4264 13 Jan 23 peter 93         suite.err() << "assigned copy not identical\n";
4264 13 Jan 23 peter 94       }
4264 13 Jan 23 peter 95     }
4264 13 Jan 23 peter 96     std::string str = header.text();
4264 13 Jan 23 peter 97     suite.add(str.size());
4264 13 Jan 23 peter 98   }
4264 13 Jan 23 peter 99
4264 13 Jan 23 peter 100
4264 13 Jan 23 peter 101   suite.out() << "test move assignment\n";
4264 13 Jan 23 peter 102   {
4264 13 Jan 23 peter 103     {
4264 13 Jan 23 peter 104       BamHeader copy(header);
4264 13 Jan 23 peter 105       BamHeader hdr;
4264 13 Jan 23 peter 106       hdr = std::move(copy);
4264 13 Jan 23 peter 107       if (hdr.text() != header.text()) {
4264 13 Jan 23 peter 108         suite.add(false);
4264 13 Jan 23 peter 109         suite.err() << "copy (from move assignment) not identical\n";
4264 13 Jan 23 peter 110       }
4264 13 Jan 23 peter 111     }
4264 13 Jan 23 peter 112     std::string str = header.text();
4264 13 Jan 23 peter 113     suite.add(str.size());
4264 13 Jan 23 peter 114   }
4264 13 Jan 23 peter 115 }
4264 13 Jan 23 peter 116
4264 13 Jan 23 peter 117
2953 16 Jan 13 peter 118 void test1(test::Suite& suite)
2953 16 Jan 13 peter 119 {
2953 16 Jan 13 peter 120   std::string file = "../../data/foo.sorted.bam";
2953 16 Jan 13 peter 121
2953 16 Jan 13 peter 122   InBamFile bam_stream(file);
2953 16 Jan 13 peter 123   BamHeader hdr = bam_stream.header();
4264 13 Jan 23 peter 124   test_copy_and_move(suite, hdr);
2953 16 Jan 13 peter 125   suite.out() << hdr.n_targets() << "\n";
2953 16 Jan 13 peter 126   suite.add(hdr.n_targets()==84);
2953 16 Jan 13 peter 127   suite.out() << hdr.target_name(0) << "\n";
2953 16 Jan 13 peter 128   if (!suite.add(std::string("1")==hdr.target_name(0)))
2953 16 Jan 13 peter 129     suite.err() << "error\n";
2999 14 Mar 13 peter 130   int32_t tid = hdr.tid("1");
2999 14 Mar 13 peter 131   suite.out() << "tid: " << tid << "\n";
2999 14 Mar 13 peter 132   if (!suite.add(tid==0))
2999 14 Mar 13 peter 133     suite.err() << "error\n";
2999 14 Mar 13 peter 134
2999 14 Mar 13 peter 135   int begin, end;
2999 14 Mar 13 peter 136   hdr.parse_region("GL000197.1", tid, begin, end);
2999 14 Mar 13 peter 137   suite.out() << tid << " " << begin << " " << end << "\n";
2999 14 Mar 13 peter 138   suite.add(tid==35);
2999 14 Mar 13 peter 139   hdr.parse_region("GL000197.1:1000-2000", tid, begin, end);
2999 14 Mar 13 peter 140   suite.out() << tid << " " << begin << " " << end << "\n";
2999 14 Mar 13 peter 141   if (!suite.add(tid==35))
2999 14 Mar 13 peter 142     suite.err() << "incorrect tid\n";
2999 14 Mar 13 peter 143   if (!suite.add(begin==999))
2999 14 Mar 13 peter 144     suite.err() << "incorrect begin\n";
2999 14 Mar 13 peter 145   if (!suite.add(end==2000))
2999 14 Mar 13 peter 146     suite.err() << "incorrect end\n";
2999 14 Mar 13 peter 147
2953 16 Jan 13 peter 148   suite.out() << hdr.target_name(35) << "\n";
2953 16 Jan 13 peter 149   if (!suite.add(std::string("GL000197.1")==hdr.target_name(35)))
2953 16 Jan 13 peter 150     suite.err() << "error\n";
2953 16 Jan 13 peter 151   suite.out() << hdr.target_length(0) << "\n";
2953 16 Jan 13 peter 152   suite.add(hdr.target_length(0)==249250621);
3408 16 Apr 15 peter 153
3408 16 Apr 15 peter 154   std::string str = hdr.text();
3408 16 Apr 15 peter 155   std::string str1(str);
3408 16 Apr 15 peter 156   str1 += "@PG\tID:prog\tVN:1.0\n";
3408 16 Apr 15 peter 157   omic::BamHeader hdr2(hdr);
3408 16 Apr 15 peter 158   hdr.text(str1);
3408 16 Apr 15 peter 159   std::string str2 = hdr.text();
3408 16 Apr 15 peter 160   if (str2.substr(0, str.size()) != str) {
3408 16 Apr 15 peter 161     suite.add(false);
3408 16 Apr 15 peter 162     suite.err() << "incorrect text:\n" << str2 << "\nexpected:\n"
3408 16 Apr 15 peter 163                 << str1 << "\n";
3408 16 Apr 15 peter 164   }
3408 16 Apr 15 peter 165   else if (str2 != str1) {
3408 16 Apr 15 peter 166     suite.add(false);
3408 16 Apr 15 peter 167     suite.err() << "error: line was not added to header\n";
3408 16 Apr 15 peter 168   }
3408 16 Apr 15 peter 169   if (hdr2.text() != str) {
3408 16 Apr 15 peter 170     suite.add(false);
3408 16 Apr 15 peter 171     suite.err() << "function text(1) affects copies\n";
3408 16 Apr 15 peter 172   }
3408 16 Apr 15 peter 173
3409 17 Apr 15 peter 174   int32_t n_chr = hdr.n_targets();
3409 17 Apr 15 peter 175
4290 03 Feb 23 peter 176   // check that we can add another chromosome via a @SG field
3409 17 Apr 15 peter 177   std::stringstream ss(hdr.text());
3409 17 Apr 15 peter 178   std::string str3;
3409 17 Apr 15 peter 179   std::string line;
3409 17 Apr 15 peter 180   int32_t n_targ = 0;
3409 17 Apr 15 peter 181   while (getline(ss, line)) {
3409 17 Apr 15 peter 182     str3 += line + "\n";
3409 17 Apr 15 peter 183     if (line.substr(0,3) == "@SQ") {
3409 17 Apr 15 peter 184       ++n_targ;
3409 17 Apr 15 peter 185       if (n_targ == n_chr)
3409 17 Apr 15 peter 186         str3 += "@SQ\tSN:Zbababababa\tLN:123\n";
3409 17 Apr 15 peter 187     }
3409 17 Apr 15 peter 188   }
3409 17 Apr 15 peter 189   hdr.text(str3);
3409 17 Apr 15 peter 190
3409 17 Apr 15 peter 191   int32_t n = hdr.n_targets();
3409 17 Apr 15 peter 192   if (n != n_chr+1) {
3409 17 Apr 15 peter 193     suite.add(false);
3409 17 Apr 15 peter 194     suite.err() << "error: incorrect number of targets: " << n
3409 17 Apr 15 peter 195                 << " expected " << (n_chr+1) << "\n";
3409 17 Apr 15 peter 196     suite.err() << "header text:\n" << hdr.text() << "===\n";
3409 17 Apr 15 peter 197   }
3409 17 Apr 15 peter 198
3410 21 Apr 15 peter 199   suite.err() << hdr.text() << "\n";
3410 21 Apr 15 peter 200   try {
3410 21 Apr 15 peter 201     std::string sample_name = hdr.read_group("foo", "SM");
3410 21 Apr 15 peter 202     if (sample_name!="Tumour") {
3410 21 Apr 15 peter 203       suite.add(false);
3410 21 Apr 15 peter 204       suite.err() << "error: incorrect sample name: '" << sample_name << "'\n";
3410 21 Apr 15 peter 205     }
3410 21 Apr 15 peter 206
3477 09 Mar 16 peter 207     typedef BamHeader::read_group_iterator iterator;
3477 09 Mar 16 peter 208     iterator it = hdr.read_group_begin();
3599 22 Jan 17 peter 209     test::test_bidirectional_traversal_iterator(it);
3599 22 Jan 17 peter 210     test::test_readable_iterator(it);
3477 09 Mar 16 peter 211     if (it == hdr.read_group_end()) {
3477 09 Mar 16 peter 212       suite.add(false);
3477 09 Mar 16 peter 213       suite.err() << "it == hdr.read_group_end()\n";
3477 09 Mar 16 peter 214     }
3477 09 Mar 16 peter 215     else {
3477 09 Mar 16 peter 216       if (!suite.add(*it == "foo"))
3477 09 Mar 16 peter 217         suite.err() << "*it != foo\n";
3477 09 Mar 16 peter 218       std::map<std::string, std::string> m = hdr.read_group(*it);
3477 09 Mar 16 peter 219       std::string val = m["SM"];
3477 09 Mar 16 peter 220       if (!suite.add(val == "Tumour"))
3477 09 Mar 16 peter 221         suite.err() << "incorrect val: " << val << " expected Tumour\n";
3477 09 Mar 16 peter 222       val = m["PL"];
3477 09 Mar 16 peter 223       if (!suite.add(val == "ILLUMINA"))
3477 09 Mar 16 peter 224         suite.err() << "incorrect val: " << val << " expected ILLUMINA\n";
3477 09 Mar 16 peter 225     }
3477 09 Mar 16 peter 226
3410 21 Apr 15 peter 227     std::string bwa_version = hdr.program_group("bwa", "VN");
3410 21 Apr 15 peter 228     if (bwa_version!="0.6.1-r104") {
3410 21 Apr 15 peter 229       suite.add(false);
3410 21 Apr 15 peter 230       suite.err() << "error: incorrect bwa_version: '" << bwa_version << "'\n";
3410 21 Apr 15 peter 231     }
3410 21 Apr 15 peter 232   }
3410 21 Apr 15 peter 233   catch (std::runtime_error& e) {
3410 21 Apr 15 peter 234     suite.add(false);
3410 21 Apr 15 peter 235     suite.err() << "error: exception: " << e.what() << "\n";
3410 21 Apr 15 peter 236   }
3412 22 Apr 15 peter 237
3412 22 Apr 15 peter 238   try {
3412 22 Apr 15 peter 239     // try calling unknown ID
3412 22 Apr 15 peter 240     hdr.read_group("banana", "SM");
3412 22 Apr 15 peter 241     suite.add(false);
3412 22 Apr 15 peter 242     suite.err() << "read_group(\"banana\", \"SM\") didn't throw\n";
3412 22 Apr 15 peter 243   }
3412 22 Apr 15 peter 244   catch (utility::runtime_error& e) {
3412 22 Apr 15 peter 245     suite.out() << "expected exception with what(): " << e.what() << "\n";
3412 22 Apr 15 peter 246   }
3412 22 Apr 15 peter 247
3412 22 Apr 15 peter 248   try {
3412 22 Apr 15 peter 249     // try calling unknown key
3412 22 Apr 15 peter 250     hdr.read_group("foo", "HM");
3412 22 Apr 15 peter 251     suite.add(false);
3412 22 Apr 15 peter 252     suite.err() << "read_group(\"foo\", \"HM\") didn't throw\n";
3412 22 Apr 15 peter 253   }
3412 22 Apr 15 peter 254   catch (utility::runtime_error& e) {
3412 22 Apr 15 peter 255     suite.out() << "expected exception with what(): " << e.what() << "\n";
3412 22 Apr 15 peter 256   }
4290 03 Feb 23 peter 257
4290 03 Feb 23 peter 258
4290 03 Feb 23 peter 259   std::set<std::string>
4290 03 Feb 23 peter 260     programs(hdr.program_group_begin(), hdr.program_group_end());
4290 03 Feb 23 peter 261   for (const auto& prog : programs)
4290 03 Feb 23 peter 262     suite.out() << prog << "\n";
4290 03 Feb 23 peter 263   if (!programs.count("bwa") || !programs.count("samtools") ||
4290 03 Feb 23 peter 264       programs.size()!=3) {
4290 03 Feb 23 peter 265     suite.add(false);
4290 03 Feb 23 peter 266     suite.err() << "incorrect programs\n";
4290 03 Feb 23 peter 267   }
4290 03 Feb 23 peter 268   hdr.program_group("bwa");
4264 13 Jan 23 peter 269 }
2953 16 Jan 13 peter 270 #endif

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