test/pileup.cc

Code
Comments
Other
Rev Date Author Line
3317 19 Sep 14 peter 1 // $Id$
3317 19 Sep 14 peter 2
3317 19 Sep 14 peter 3 /*
3999 08 Oct 20 peter 4   Copyright (C) 2014, 2015, 2019, 2020 Peter Johansson
3317 19 Sep 14 peter 5
3317 19 Sep 14 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3317 19 Sep 14 peter 7
3317 19 Sep 14 peter 8   The yat library is free software; you can redistribute it and/or
3317 19 Sep 14 peter 9   modify it under the terms of the GNU General Public License as
3317 19 Sep 14 peter 10   published by the Free Software Foundation; either version 3 of the
3317 19 Sep 14 peter 11   License, or (at your option) any later version.
3317 19 Sep 14 peter 12
3317 19 Sep 14 peter 13   The yat library is distributed in the hope that it will be useful,
3317 19 Sep 14 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3317 19 Sep 14 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3317 19 Sep 14 peter 16   General Public License for more details.
3317 19 Sep 14 peter 17
3317 19 Sep 14 peter 18   You should have received a copy of the GNU General Public License
3317 19 Sep 14 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3317 19 Sep 14 peter 20 */
3317 19 Sep 14 peter 21
3317 19 Sep 14 peter 22 #include <config.h>
3317 19 Sep 14 peter 23
3317 19 Sep 14 peter 24 #include "Suite.h"
3317 19 Sep 14 peter 25
3439 20 Nov 15 peter 26 #include <yat/utility/config_public.h>
3439 20 Nov 15 peter 27
3883 24 Mar 20 peter 28 #ifdef YAT_HAVE_HTSLIB
3317 19 Sep 14 peter 29 #include <yat/omic/BamFile.h>
3317 19 Sep 14 peter 30 #include <yat/omic/BamReadIterator.h>
3317 19 Sep 14 peter 31 #include "yat/omic/Pileup.h"
3369 10 Feb 15 peter 32 #endif
3317 19 Sep 14 peter 33
3317 19 Sep 14 peter 34 #include <yat/utility/Aligner.h>
3317 19 Sep 14 peter 35
3317 19 Sep 14 peter 36 #include <algorithm>
3317 19 Sep 14 peter 37 #include <cassert>
3317 19 Sep 14 peter 38 #include <string>
3317 19 Sep 14 peter 39 #include <vector>
3317 19 Sep 14 peter 40
3317 19 Sep 14 peter 41 using namespace theplu::yat;
3317 19 Sep 14 peter 42
3317 19 Sep 14 peter 43 // The simplest of genotypers; reporting the first and second "thirdiles"
3317 19 Sep 14 peter 44 std::string genotype(std::string str)
3317 19 Sep 14 peter 45 {
3317 19 Sep 14 peter 46   std::string res("  ");
3317 19 Sep 14 peter 47   std::sort(str.begin(), str.end());
3317 19 Sep 14 peter 48   res[0] = str[str.size()/3];
3317 19 Sep 14 peter 49   res[1] = str[2*str.size()/3];
3317 19 Sep 14 peter 50   return res;
3317 19 Sep 14 peter 51 }
3317 19 Sep 14 peter 52
3317 19 Sep 14 peter 53
3888 25 Mar 20 peter 54 #ifdef YAT_HAVE_HTSLIB
3356 22 Nov 14 peter 55 char seq_nt16(size_t x)
3356 22 Nov 14 peter 56 {
3356 22 Nov 14 peter 57   return seq_nt16_str[x];
3356 22 Nov 14 peter 58 }
3356 22 Nov 14 peter 59
3888 25 Mar 20 peter 60
3317 19 Sep 14 peter 61 void test1(test::Suite& suite, const std::string& fn)
3317 19 Sep 14 peter 62 {
3317 19 Sep 14 peter 63   omic::InBamFile in;
3317 19 Sep 14 peter 64   in.open(fn);
3317 19 Sep 14 peter 65
3317 19 Sep 14 peter 66   using omic::BamReadIterator;
3317 19 Sep 14 peter 67   BamReadIterator first(in);
3317 19 Sep 14 peter 68   BamReadIterator last;
3317 19 Sep 14 peter 69
3317 19 Sep 14 peter 70   size_t count = 0;
3317 19 Sep 14 peter 71   std::vector<std::string> correct_gt;
3317 19 Sep 14 peter 72   correct_gt.push_back("AA");
3317 19 Sep 14 peter 73   correct_gt.push_back("CC");
3317 19 Sep 14 peter 74   correct_gt.push_back("CC");
3317 19 Sep 14 peter 75   correct_gt.push_back("TT");
3317 19 Sep 14 peter 76   correct_gt.push_back("=A");
3317 19 Sep 14 peter 77   correct_gt.push_back("AA");
3317 19 Sep 14 peter 78   correct_gt.push_back("AA");
3317 19 Sep 14 peter 79   correct_gt.push_back("GG");
3317 19 Sep 14 peter 80   correct_gt.push_back("AA");
3317 19 Sep 14 peter 81   correct_gt.push_back("AA");
3317 19 Sep 14 peter 82   correct_gt.push_back("AA");
3317 19 Sep 14 peter 83
3317 19 Sep 14 peter 84   correct_gt.push_back("GG");
3317 19 Sep 14 peter 85   correct_gt.push_back("GG");
3317 19 Sep 14 peter 86   correct_gt.push_back("CC");
3317 19 Sep 14 peter 87   correct_gt.push_back("CC");
3317 19 Sep 14 peter 88   correct_gt.push_back("CC");
3317 19 Sep 14 peter 89
3317 19 Sep 14 peter 90   correct_gt.push_back("TT");
3317 19 Sep 14 peter 91   correct_gt.push_back("TT");
3317 19 Sep 14 peter 92   correct_gt.push_back("=G");
3317 19 Sep 14 peter 93   correct_gt.push_back("=C");
3317 19 Sep 14 peter 94
3317 19 Sep 14 peter 95   using omic::Pileup;
3317 19 Sep 14 peter 96   Pileup<BamReadIterator> pileup(first, last);
3317 19 Sep 14 peter 97   while (pileup.good()) {
3317 19 Sep 14 peter 98     if (pileup.pos()+1 < 24341) {
3317 19 Sep 14 peter 99       pileup.increment();
3317 19 Sep 14 peter 100       continue;
3317 19 Sep 14 peter 101     }
3317 19 Sep 14 peter 102     Pileup<BamReadIterator>::const_iterator iter = pileup.begin();
3317 19 Sep 14 peter 103     std::string str;
3317 19 Sep 14 peter 104     for (; iter!=pileup.end(); ++iter) {
3327 10 Oct 14 peter 105       if (iter->bam().end() <= pileup.pos()) {
3327 10 Oct 14 peter 106         suite.err() << "error: " << pileup.pos() << " "
3327 10 Oct 14 peter 107                     << iter->bam().pos() << " "
3327 10 Oct 14 peter 108                     << iter->bam().end() << " "
3327 10 Oct 14 peter 109                     << iter->bam().cigar_str()
3327 10 Oct 14 peter 110                     << "\n";
3327 10 Oct 14 peter 111         suite.add(false);
3317 19 Sep 14 peter 112         continue;
3327 10 Oct 14 peter 113       }
3356 22 Nov 14 peter 114       str += seq_nt16(iter->sequence());
3317 19 Sep 14 peter 115     }
3317 19 Sep 14 peter 116
3317 19 Sep 14 peter 117     if (pileup.pos()+1 == 24341 && count!=0) {
3317 19 Sep 14 peter 118       suite.err() << "incorrect count\n";
3317 19 Sep 14 peter 119       suite.add(false);
3317 19 Sep 14 peter 120     }
3317 19 Sep 14 peter 121
3317 19 Sep 14 peter 122     if (count < correct_gt.size())
3317 19 Sep 14 peter 123       if (genotype(str) != correct_gt[count]) {
3317 19 Sep 14 peter 124         suite.err() << "incorrect genotype: " << genotype(str) << "\n";
3317 19 Sep 14 peter 125         suite.err() << "expected: " << correct_gt[count] << "\n";
3317 19 Sep 14 peter 126         suite.add(false);
3317 19 Sep 14 peter 127       }
3317 19 Sep 14 peter 128
3317 19 Sep 14 peter 129     if (count <= 3) {
3317 19 Sep 14 peter 130       if (pileup.pos()+1 != static_cast<int>(count+24341)) {
3317 19 Sep 14 peter 131         suite.err() << "incorrect pos or count\n";
3317 19 Sep 14 peter 132         suite.err() << "pos: " << pileup.pos()+1 << "\n";
3317 19 Sep 14 peter 133         suite.err() << "count: " << count << "\n";
3317 19 Sep 14 peter 134         suite.add(false);
3317 19 Sep 14 peter 135       }
3317 19 Sep 14 peter 136     }
3317 19 Sep 14 peter 137     else if (count == 4 && !pileup.skip_ref() ) {
3317 19 Sep 14 peter 138       suite.out() << "expected skip ref\n";
3317 19 Sep 14 peter 139       suite.add(false);
3317 19 Sep 14 peter 140     }
3317 19 Sep 14 peter 141     else if (count==5) {
3317 19 Sep 14 peter 142       if (pileup.pos()+1 != static_cast<int>(count+24340)) {
3317 19 Sep 14 peter 143         suite.err() << "incorrect pos or count\n";
3317 19 Sep 14 peter 144         suite.err() << "pos: " << pileup.pos()+1 << "\n";
3317 19 Sep 14 peter 145         suite.err() << "count: " << count << "\n";
3317 19 Sep 14 peter 146         suite.add(false);
3317 19 Sep 14 peter 147       }
3317 19 Sep 14 peter 148     }
3317 19 Sep 14 peter 149     pileup.increment();
3317 19 Sep 14 peter 150     ++count;
3317 19 Sep 14 peter 151   }
3801 04 May 19 peter 152   Pileup<BamReadIterator>::const_reverse_iterator rend = pileup.rend();
3801 04 May 19 peter 153   std::distance(pileup.rbegin(), rend);
3802 12 May 19 peter 154   pileup.size();
3317 19 Sep 14 peter 155 }
3317 19 Sep 14 peter 156
3317 19 Sep 14 peter 157
3317 19 Sep 14 peter 158 void test2(test::Suite& suite, const std::string& fn)
3317 19 Sep 14 peter 159 {
3317 19 Sep 14 peter 160   std::ostringstream error;
3317 19 Sep 14 peter 161   theplu::yat::omic::InBamFile in(fn);
3317 19 Sep 14 peter 162   using omic::BamRead;
3317 19 Sep 14 peter 163   using omic::BamReadIterator;
3317 19 Sep 14 peter 164   using omic::Pileup;
3317 19 Sep 14 peter 165   BamReadIterator first(in, 1, 24150, 25000);
3317 19 Sep 14 peter 166   BamReadIterator last;
3317 19 Sep 14 peter 167
3317 19 Sep 14 peter 168   std::vector<BamRead> reads(2);
3317 19 Sep 14 peter 169   reads[0] = *first;
3317 19 Sep 14 peter 170   ++first;
3317 19 Sep 14 peter 171   reads[1] = *first;
3317 19 Sep 14 peter 172   in.close();
3317 19 Sep 14 peter 173   std::string ref = reads[0].sequence();
3317 19 Sep 14 peter 174
3317 19 Sep 14 peter 175   suite.out() << "no modifications\n";
3317 19 Sep 14 peter 176   // no modification
3317 19 Sep 14 peter 177   typedef std::vector<BamRead>::iterator iterator;
3317 19 Sep 14 peter 178   typedef Pileup<iterator>::const_iterator plp_iterator;
3317 19 Sep 14 peter 179   for (Pileup<iterator> plp(reads.begin(), reads.end());plp.pos()<24070;
3317 19 Sep 14 peter 180        plp.increment()) {
3317 19 Sep 14 peter 181     for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
3317 19 Sep 14 peter 182       if (i->bam().pos()!=reads[1].pos())
3317 19 Sep 14 peter 183         continue;
3356 22 Nov 14 peter 184       char nt = seq_nt16(i->sequence());
3317 19 Sep 14 peter 185       if (nt != ref[plp.pos()-reads[0].pos()]) {
3317 19 Sep 14 peter 186         error << "'" << nt << "' not equal '"
3317 19 Sep 14 peter 187               << ref[plp.pos()-reads[0].pos()]
3317 19 Sep 14 peter 188               << "' at pos " << reads[0].pos() << " for bam with start pos: "
3317 19 Sep 14 peter 189               << i->bam().pos();
3317 19 Sep 14 peter 190         throw std::runtime_error(error.str());
3317 19 Sep 14 peter 191       }
3317 19 Sep 14 peter 192     }
3317 19 Sep 14 peter 193   }
3317 19 Sep 14 peter 194
3317 19 Sep 14 peter 195   // insertion
3317 19 Sep 14 peter 196   suite.out() << "insertion\n";
3317 19 Sep 14 peter 197   assert(reads[1].sequence_length()==100);
3317 19 Sep 14 peter 198   utility::Aligner::Cigar cig;
3317 19 Sep 14 peter 199   cig.push_back(BAM_CMATCH, 10);
3317 19 Sep 14 peter 200   cig.push_back(BAM_CINS, 1);
3317 19 Sep 14 peter 201   cig.push_back(BAM_CMATCH, 89);
3317 19 Sep 14 peter 202   reads[1].cigar(cig);
3317 19 Sep 14 peter 203   int count = 0;
3317 19 Sep 14 peter 204   for (Pileup<iterator> plp(reads.begin(), reads.end());plp.good();
3317 19 Sep 14 peter 205        plp.increment()) {
3317 19 Sep 14 peter 206     for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
3317 19 Sep 14 peter 207       if (!same_query_name(i->bam(), reads[1]))
3317 19 Sep 14 peter 208         continue;
3317 19 Sep 14 peter 209       ++count;
3317 19 Sep 14 peter 210       suite.err() << count << "\n";
3356 22 Nov 14 peter 211       char nt = seq_nt16(i->sequence());
3317 19 Sep 14 peter 212       if (nt != reads[1].sequence()[count-1])
3317 19 Sep 14 peter 213         error << "nt: " << nt << " not as expected. count: " << count << "\n";
3317 19 Sep 14 peter 214       // first 10 matches
3317 19 Sep 14 peter 215       if (count <= 10) {
3317 19 Sep 14 peter 216         if (plp.skip_ref())
3317 19 Sep 14 peter 217           error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 218         if (plp.pos() != reads[1].pos()+count-1)
3317 19 Sep 14 peter 219           error << "count: " << count << " with pos: " << plp.pos() << "\n";
3317 19 Sep 14 peter 220       }
3317 19 Sep 14 peter 221       else if (count==11) { // insertion
3317 19 Sep 14 peter 222         if (!plp.skip_ref())
3317 19 Sep 14 peter 223           error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 224       }
3317 19 Sep 14 peter 225       else {
3317 19 Sep 14 peter 226         if (plp.skip_ref())
3317 19 Sep 14 peter 227           error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 228         if (plp.pos() != reads[1].pos()+count-2)
3317 19 Sep 14 peter 229           error << "count: " << count << " with pos: " << plp.pos() << "\n";
3317 19 Sep 14 peter 230       }
3317 19 Sep 14 peter 231
3317 19 Sep 14 peter 232       if (error.str().size())
3317 19 Sep 14 peter 233         throw std::runtime_error(error.str());
3317 19 Sep 14 peter 234     }
3317 19 Sep 14 peter 235   }
3317 19 Sep 14 peter 236
3317 19 Sep 14 peter 237
3317 19 Sep 14 peter 238   // deletion
3317 19 Sep 14 peter 239   suite.out() << "deletion\n";
3317 19 Sep 14 peter 240   assert(reads[1].sequence_length()==100);
3317 19 Sep 14 peter 241   cig.clear();
3317 19 Sep 14 peter 242   cig.push_back(BAM_CMATCH, 10);
3317 19 Sep 14 peter 243   cig.push_back(BAM_CDEL, 1);
3317 19 Sep 14 peter 244   cig.push_back(BAM_CMATCH, 90);
3317 19 Sep 14 peter 245   assert(cig.size()==3);
3317 19 Sep 14 peter 246   reads[1].cigar(cig);
3317 19 Sep 14 peter 247   count = 0;
3317 19 Sep 14 peter 248   for (Pileup<iterator> plp(reads.begin(), reads.end()); plp.good();
3317 19 Sep 14 peter 249        plp.increment()) {
3317 19 Sep 14 peter 250     for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
3317 19 Sep 14 peter 251       if (!same_query_name(i->bam(), reads[1]))
3317 19 Sep 14 peter 252         continue;
3317 19 Sep 14 peter 253       ++count;
3356 22 Nov 14 peter 254       char nt = seq_nt16(i->sequence());
3317 19 Sep 14 peter 255       if (plp.pos() != reads[1].pos()+count-1)
3317 19 Sep 14 peter 256         error << "count: " << count << " with pos: " << plp.pos() << "\n";
3317 19 Sep 14 peter 257       // first 10 matches
3317 19 Sep 14 peter 258       if (count <= 10) {
3317 19 Sep 14 peter 259         if (plp.skip_ref())
3317 19 Sep 14 peter 260           error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 261         if (nt != reads[1].sequence()[count-1])
3317 19 Sep 14 peter 262           error << "nt: " << nt << " not as expected\n";
3317 19 Sep 14 peter 263       }
3317 19 Sep 14 peter 264       else if (count==11) {
3317 19 Sep 14 peter 265         if (plp.skip_ref())
3317 19 Sep 14 peter 266           error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 267         if (!i->is_del())
3317 19 Sep 14 peter 268           error << count << " expected deletion\n";
3317 19 Sep 14 peter 269       }
3317 19 Sep 14 peter 270       else {
3317 19 Sep 14 peter 271         if (plp.skip_ref())
3317 19 Sep 14 peter 272           error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 273         if (nt != reads[1].sequence()[count-2])
3317 19 Sep 14 peter 274           error << "nt: " << nt << " not as expected\n";
3317 19 Sep 14 peter 275       }
3317 19 Sep 14 peter 276
3317 19 Sep 14 peter 277       if (error.str().size())
3317 19 Sep 14 peter 278         throw std::runtime_error(error.str());
3317 19 Sep 14 peter 279     }
3317 19 Sep 14 peter 280   }
3317 19 Sep 14 peter 281
3317 19 Sep 14 peter 282   // soft clipped
3317 19 Sep 14 peter 283   suite.out() << "soft clipped\n";
3317 19 Sep 14 peter 284   assert(reads[1].sequence_length()==100);
3317 19 Sep 14 peter 285   cig.clear();
3317 19 Sep 14 peter 286   cig.push_back(BAM_CSOFT_CLIP, 10);
3317 19 Sep 14 peter 287   cig.push_back(BAM_CMATCH, 90);
3317 19 Sep 14 peter 288   reads[1].cigar(cig);
3317 19 Sep 14 peter 289   count = 0;
3317 19 Sep 14 peter 290   for (Pileup<iterator> plp(reads.begin(), reads.end()); plp.good();
3317 19 Sep 14 peter 291        plp.increment()) {
3317 19 Sep 14 peter 292     for (plp_iterator i = plp.begin(); i!=plp.end(); ++i) {
3317 19 Sep 14 peter 293       if (!same_query_name(i->bam(), reads[1]))
3317 19 Sep 14 peter 294         continue;
3317 19 Sep 14 peter 295       ++count;
3356 22 Nov 14 peter 296       char nt = seq_nt16(i->sequence());
3317 19 Sep 14 peter 297       if (plp.skip_ref())
3317 19 Sep 14 peter 298         error << count << " unexpected skip_ref\n";
3317 19 Sep 14 peter 299       if (plp.pos() != reads[1].pos()+count-1)
3317 19 Sep 14 peter 300         error << "count: " << count << " with pos: " << plp.pos() << "\n";
3317 19 Sep 14 peter 301       if (nt != reads[1].sequence()[count+9])
3317 19 Sep 14 peter 302         error << "count: " << count << " nt: " << nt << " not as expected\n";
3317 19 Sep 14 peter 303
3317 19 Sep 14 peter 304       if (error.str().size())
3317 19 Sep 14 peter 305         throw std::runtime_error(error.str());
3317 19 Sep 14 peter 306     }
3317 19 Sep 14 peter 307   }
3317 19 Sep 14 peter 308 }
3317 19 Sep 14 peter 309 #else
3317 19 Sep 14 peter 310 // tests if libbam is not available, should never be run
3317 19 Sep 14 peter 311 void test1(test::Suite& suite, const std::string& fn) {assert(0);}
3317 19 Sep 14 peter 312 void test2(test::Suite& suite, const std::string& fn) {assert(0);}
3317 19 Sep 14 peter 313 #endif
3317 19 Sep 14 peter 314
3317 19 Sep 14 peter 315 int main(int argc, char* argv[])
3317 19 Sep 14 peter 316 {
3317 19 Sep 14 peter 317   test::Suite suite(argc, argv, true);
3317 19 Sep 14 peter 318   std::string fn = "../../data/foo.sorted.bam";
3317 19 Sep 14 peter 319   try {
3317 19 Sep 14 peter 320     test1(suite, fn);
3317 19 Sep 14 peter 321     test2(suite, fn);
3317 19 Sep 14 peter 322   }
3317 19 Sep 14 peter 323   catch (std::runtime_error& e) {
3317 19 Sep 14 peter 324     suite.err() << "test failed: " << e.what() << "\n";
3317 19 Sep 14 peter 325     suite.add(false);
3317 19 Sep 14 peter 326   }
3317 19 Sep 14 peter 327   return suite.return_value();
3317 19 Sep 14 peter 328 }

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