From cedf5b4448388ad1593143b9f52bea8a3194c930 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Mon, 3 Oct 2016 00:18:30 +0100 Subject: [PATCH 01/29] tmp change to make the r version to compile --- src/dEploid.cpp | 2 +- src/dEploidIO.cpp | 31 +++++++++++++++++++++++++++---- src/dEploidIO.hpp | 8 ++++++-- 3 files changed, 34 insertions(+), 7 deletions(-) diff --git a/src/dEploid.cpp b/src/dEploid.cpp index dd3648e6..feeedef8 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -32,7 +32,7 @@ int main( int argc, char *argv[] ){ try { DEploidIO dEploidIO; - (void)dEploidIO.core( argc, argv ); + (void)dEploidIO.core(); if ( dEploidIO.version() ){ dEploidIO.printVersion(); diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index ac21b40a..00ce7173 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -26,6 +26,7 @@ #include // assert #include // std::setw #include +#include DEploidIO::DEploidIO(){ @@ -33,6 +34,25 @@ DEploidIO::DEploidIO(){ } +DEploidIO::DEploidIO(const std::string &arg) { + this->init(); + std::istringstream iss(arg); + copy(std::istream_iterator(iss), + std::istream_iterator(), + std::back_inserter(argv_)); + this->core(); +} + + + +DEploidIO::DEploidIO(int argc, char *argv[]) { + this->init(); + argv_ = std::vector(argv + 1, argv + argc); + this->argv_i = argv_.begin(); + this->core(); +} + + DEploidIO::~DEploidIO(){ if ( this->excludedMarkers != NULL ){ delete this->excludedMarkers; @@ -45,10 +65,7 @@ DEploidIO::~DEploidIO(){ } -void DEploidIO::core(int argc, char *argv[]) { - argv_ = std::vector(argv + 1, argv + argc); - this->argv_i = argv_.begin(); - +void DEploidIO::core() { if ( argv_.size() == 0 ) { this->set_help(true); return; @@ -432,3 +449,9 @@ void DEploidIO::printHelp(){ cout << "./dEploid -vcf data/testData/PG0390-C.test.vcf -plaf data/testData/labStrains.test.PLAF.txt -o PG0390-CNopanel -noPanel"<< endl; cout << "./dEploid -vcf data/testData/PG0390-C.test.vcf -exclude data/testData/labStrains.test.exclude.txt -plaf data/testData/labStrains.test.PLAF.txt -o PG0390-CNopanelExclude -noPanel"<< endl; } + + +std::ostream& operator<< (std::ostream& stream, const DEploidIO& dEploidIO) { + for (std::string arg : dEploidIO.argv_) stream << " " << arg; + return stream; +} diff --git a/src/dEploidIO.hpp b/src/dEploidIO.hpp index 1b9955f3..923694a3 100644 --- a/src/dEploidIO.hpp +++ b/src/dEploidIO.hpp @@ -50,9 +50,11 @@ class DEploidIO{ friend class McmcMachinery; public: DEploidIO(); + DEploidIO(const std::string &arg); + DEploidIO(int argc, char *argv[]); ~DEploidIO (); - void core(int argc, char *argv[]); + void core(); void printHelp(); bool help() const { return help_; } void printVersion(); @@ -80,6 +82,9 @@ class DEploidIO{ // Log void write (McmcSample * mcmcSample, Panel * panel ); + bool randomSeedWasSet() const {return this->randomSeedWasSet_; } + + friend std::ostream& operator<< (std::ostream& stream, const DEploidIO& dEploidIO); private: @@ -94,7 +99,6 @@ class DEploidIO{ size_t randomSeed() const { return randomSeed_;} bool randomSeedWasSet_; void setRandomSeedWasSet(const bool random){ this->randomSeedWasSet_ = random; } - bool randomSeedWasSet() const {return this->randomSeedWasSet_; } bool initialPropWasGiven_; From 701f3102f27d0ae4335a9cc7d083cc5c0f004ee2 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Mon, 3 Oct 2016 09:16:07 +0100 Subject: [PATCH 02/29] update unit test --- .gitignore | 7 + src/dEploid.cpp | 4 +- src/dEploidIO.cpp | 1 - tests/unittest/test_dEploidIO.cpp | 331 +++++++++++++++--------------- 4 files changed, 177 insertions(+), 166 deletions(-) diff --git a/.gitignore b/.gitignore index fabba440..4414ba3e 100644 --- a/.gitignore +++ b/.gitignore @@ -22,6 +22,10 @@ *.a *.lib +# Backup file generated while rewriting r scripts +utilities/dataExplore.r.bak +utilities/interpretDEploid.r.bak + # Executables *.exe *.out @@ -79,3 +83,6 @@ vcf.* *.gcov libbash_test.info test_coverage/ + +# Doxyfile +docs/doxygen/Doxyfile diff --git a/src/dEploid.cpp b/src/dEploid.cpp index feeedef8..887eb800 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -30,9 +30,7 @@ using namespace std; int main( int argc, char *argv[] ){ try { - - DEploidIO dEploidIO; - (void)dEploidIO.core(); + DEploidIO dEploidIO(argc, argv); if ( dEploidIO.version() ){ dEploidIO.printVersion(); diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index 00ce7173..7dc5bdad 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -44,7 +44,6 @@ DEploidIO::DEploidIO(const std::string &arg) { } - DEploidIO::DEploidIO(int argc, char *argv[]) { this->init(); argv_ = std::vector(argv + 1, argv + argc); diff --git a/tests/unittest/test_dEploidIO.cpp b/tests/unittest/test_dEploidIO.cpp index 9c7dbad8..8d048406 100644 --- a/tests/unittest/test_dEploidIO.cpp +++ b/tests/unittest/test_dEploidIO.cpp @@ -30,46 +30,45 @@ class TestIO : public CppUnit::TestCase { CPPUNIT_TEST_SUITE_END(); private: - DEploidIO* input_; double epsilon3; public: void setUp() { this->epsilon3 = 0.000000000001; - this->input_ = new DEploidIO (); } void tearDown() { - delete input_; } void testInitialization(){ - CPPUNIT_ASSERT_EQUAL( this->input_->randomSeedWasSet(), false ); - CPPUNIT_ASSERT_EQUAL( this->input_->doExportRecombProb(), false ); - CPPUNIT_ASSERT_EQUAL( this->input_->compressVcf(), false ); - CPPUNIT_ASSERT_EQUAL( this->input_->initialPropWasGiven(), false ); - CPPUNIT_ASSERT_EQUAL( this->input_->excludeSites() , false ); - CPPUNIT_ASSERT( this->input_->excludedMarkers == NULL); - CPPUNIT_ASSERT_EQUAL( this->input_->randomSeed(), (size_t)0 ); - CPPUNIT_ASSERT_EQUAL( this->input_->help(), false); - CPPUNIT_ASSERT_EQUAL( this->input_->usePanel(), true); - CPPUNIT_ASSERT_EQUAL( this->input_->precision_ , (size_t)8); - CPPUNIT_ASSERT( this->input_->prefix_ == "pf3k-dEploid"); - CPPUNIT_ASSERT_EQUAL( this->input_->kStrain_ , (size_t)5); - CPPUNIT_ASSERT_EQUAL( this->input_->nMcmcSample_ , (size_t)800); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateProp() , true); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdatePair() , true); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateSingle() , true); - CPPUNIT_ASSERT_EQUAL( this->input_->doExportPostProb() , false); - CPPUNIT_ASSERT_DOUBLES_EQUAL( this->input_->mcmcBurn_, 0.5, epsilon3); - CPPUNIT_ASSERT_EQUAL( this->input_->mcmcMachineryRate_ , (size_t)5); - CPPUNIT_ASSERT_DOUBLES_EQUAL( this->input_->missCopyProb_ , 0.01, epsilon3); - CPPUNIT_ASSERT_EQUAL( this->input_->useConstRecomb() , false ); - CPPUNIT_ASSERT_EQUAL( this->input_->forbidCopyFromSame() , false ); - CPPUNIT_ASSERT_DOUBLES_EQUAL( this->input_->constRecombProb_ , 1.0, epsilon3); - CPPUNIT_ASSERT_DOUBLES_EQUAL( this->input_->averageCentimorganDistance_, 15000.0, epsilon3 ); - CPPUNIT_ASSERT_DOUBLES_EQUAL( this->input_->Ne_ , 10.0, epsilon3 ); - CPPUNIT_ASSERT_EQUAL( this->input_->doExportVcf(), false ); + DEploidIO* dEploidIOptr = new DEploidIO (); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->randomSeedWasSet(), false ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->doExportRecombProb(), false ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->compressVcf(), false ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->initialPropWasGiven(), false ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->excludeSites() , false ); + CPPUNIT_ASSERT( dEploidIOptr->excludedMarkers == NULL); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->randomSeed(), (size_t)0 ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->help(), false); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->usePanel(), true); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->precision_ , (size_t)8); + CPPUNIT_ASSERT( dEploidIOptr->prefix_ == "pf3k-dEploid"); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->kStrain_ , (size_t)5); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->nMcmcSample_ , (size_t)800); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->doUpdateProp() , true); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->doUpdatePair() , true); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->doUpdateSingle() , true); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->doExportPostProb() , false); + CPPUNIT_ASSERT_DOUBLES_EQUAL( dEploidIOptr->mcmcBurn_, 0.5, epsilon3); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->mcmcMachineryRate_ , (size_t)5); + CPPUNIT_ASSERT_DOUBLES_EQUAL( dEploidIOptr->missCopyProb_ , 0.01, epsilon3); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->useConstRecomb() , false ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->forbidCopyFromSame() , false ); + CPPUNIT_ASSERT_DOUBLES_EQUAL( dEploidIOptr->constRecombProb_ , 1.0, epsilon3); + CPPUNIT_ASSERT_DOUBLES_EQUAL( dEploidIOptr->averageCentimorganDistance_, 15000.0, epsilon3 ); + CPPUNIT_ASSERT_DOUBLES_EQUAL( dEploidIOptr->Ne_ , 10.0, epsilon3 ); + CPPUNIT_ASSERT_EQUAL( dEploidIOptr->doExportVcf(), false ); + delete dEploidIOptr; } @@ -82,10 +81,10 @@ class TestIO : public CppUnit::TestCase { "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1", "-exclude", "data/testData/labStrains.test.exclude.txt", "-vcfOut", "-z"}; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(11, argv1) ); - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(13, argv1) ); - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(14, argv1) ); - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(15, argv1) ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(11, argv1) ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(13, argv1) ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(14, argv1) ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(15, argv1) ); char *argv2[] = { "./dEploid", "-vcf", "data/testData/PG0390-C.test.vcf", @@ -93,8 +92,8 @@ class TestIO : public CppUnit::TestCase { "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1", "-exclude", "data/testData/labStrains.test.exclude.txt"}; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(9, argv2) ); - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(11, argv2) ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(9, argv2) ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(11, argv2) ); } @@ -105,8 +104,8 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-initialP", "0.1", "0.2", "0.3", "0.4", "-k", "4" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(16, argv1) ); - CPPUNIT_ASSERT_THROW ( this->input_->core(14, argv1), NumOfPropNotMatchNumStrain ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(16, argv1) ); + CPPUNIT_ASSERT_THROW ( DEploidIO(14, argv1), NumOfPropNotMatchNumStrain ); char *argv2[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -114,7 +113,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-initialP", "0.1", "0.2", "0.3", "0.3", "0.2"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(15, argv2), SumOfPropNotOne ); + CPPUNIT_ASSERT_THROW ( DEploidIO(15, argv2), SumOfPropNotOne ); char *argv3[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -122,7 +121,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-initialP", "0.1", "0.2", "0.3", "0.3", "0.02"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(15, argv3), SumOfPropNotOne ); + CPPUNIT_ASSERT_THROW ( DEploidIO(15, argv3), SumOfPropNotOne ); char *argv4[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -130,7 +129,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-initialP", "-o", "tmp"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(12, argv4), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(12, argv4), NotEnoughArg ); char *argv5[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -138,55 +137,57 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-initialP"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv5), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv5), NotEnoughArg ); } void testPrintHelp(){ char *argv[] = { "./dEploid" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(1, argv) ); - - CPPUNIT_ASSERT_EQUAL((size_t)0, input_->argv_.size()); - CPPUNIT_ASSERT_EQUAL( true, input_->help()); - CPPUNIT_ASSERT_NO_THROW ( this->input_->printHelp() ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(1, argv) ); + DEploidIO dEploidIO1(1, argv); + CPPUNIT_ASSERT_EQUAL((size_t)0, dEploidIO1.argv_.size()); + CPPUNIT_ASSERT_EQUAL( true, dEploidIO1.help()); + CPPUNIT_ASSERT_NO_THROW ( dEploidIO1.printHelp() ); char *argv1[] = { "./dEploid", "-h" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(2, argv1) ); - - CPPUNIT_ASSERT_EQUAL((size_t)1, input_->argv_.size()); - CPPUNIT_ASSERT_EQUAL( true, input_->help()); - CPPUNIT_ASSERT_NO_THROW ( this->input_->printHelp() ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(2, argv1) ); + DEploidIO dEploidIO2(2, argv1); + CPPUNIT_ASSERT_EQUAL((size_t)1, dEploidIO2.argv_.size()); + CPPUNIT_ASSERT_EQUAL( true, dEploidIO2.help()); + CPPUNIT_ASSERT_NO_THROW ( dEploidIO2.printHelp() ); char *argv2[] = { "./dEploid", "-help" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(2, argv2) ); - - CPPUNIT_ASSERT_EQUAL((size_t)1, input_->argv_.size()); - CPPUNIT_ASSERT_EQUAL( true, input_->help()); - CPPUNIT_ASSERT_NO_THROW ( this->input_->printHelp() ); + CPPUNIT_ASSERT_NO_THROW(DEploidIO(2, argv2)); + DEploidIO dEploidIO3(2, argv2); + CPPUNIT_ASSERT_EQUAL((size_t)1, dEploidIO3.argv_.size()); + CPPUNIT_ASSERT_EQUAL(true, dEploidIO3.help()); + CPPUNIT_ASSERT_NO_THROW(dEploidIO3.printHelp()); } void testPrintVersion(){ char *argv1[] = { "./dEploid", "-v" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(2, argv1) ); - CPPUNIT_ASSERT_EQUAL( true, input_->version()); - CPPUNIT_ASSERT_NO_THROW ( this->input_->printVersion() ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(2, argv1) ); + DEploidIO dEploidIO1(2, argv1); + CPPUNIT_ASSERT_EQUAL( true, dEploidIO1.version()); + CPPUNIT_ASSERT_NO_THROW ( dEploidIO1.printVersion() ); char *argv2[] = { "./dEploid", "-v" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(2, argv2) ); - CPPUNIT_ASSERT_EQUAL( true, input_->version()); - CPPUNIT_ASSERT_NO_THROW ( this->input_->printVersion() ); + CPPUNIT_ASSERT_NO_THROW ( DEploidIO(2, argv2) ); + DEploidIO dEploidIO2(2, argv2); + CPPUNIT_ASSERT_EQUAL( true, dEploidIO2.version()); + CPPUNIT_ASSERT_NO_THROW ( dEploidIO2.printVersion() ); } void testOutOfRange(){ - char *argv1[] = { "./dEploid", + char *argv1[] = {"./dEploid", "-ref", "data/testData/PG0390-C.test.ref", "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-burn", "1.1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv1), OutOfRange ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv1), OutOfRange ); char *argv2[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -194,7 +195,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-miss", "5.1"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv2), OutOfRange ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv2), OutOfRange ); char *argv3[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -202,7 +203,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-recomb", "2.1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv3), OutOfRange ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv3), OutOfRange ); } @@ -213,7 +214,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv1), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv1), NotEnoughArg ); char *argv2[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -221,7 +222,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-seed"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv2), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv2), NotEnoughArg ); char *argv3[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -229,7 +230,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-p"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv3), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv3), NotEnoughArg ); char *argv4[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -237,35 +238,35 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-k"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv4), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv4), NotEnoughArg ); char *argv5[] = { "./dEploid", "-ref", "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(8, argv5), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(8, argv5), NotEnoughArg ); char *argv6[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", "-alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(8, argv6), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(8, argv6), NotEnoughArg ); char *argv7[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "-panel", "data/testData/labStrains.test.panel.txt" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(8, argv7), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(8, argv7), NotEnoughArg ); char *argv8[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(8, argv8), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(8, argv8), NotEnoughArg ); char *argv9[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -273,7 +274,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-exclude" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv9), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv9), NotEnoughArg ); } @@ -285,8 +286,8 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-seed", "asdf"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv2), NotEnoughArg ); - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv2), WrongType ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv2), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv2), WrongType ); char *argv3[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -294,8 +295,8 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-p", "asdf"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv3), NotEnoughArg ); - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv3), WrongType ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv3), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv3), WrongType ); char *argv4[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -303,8 +304,8 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-k", "asdf"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv4), NotEnoughArg ); - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv4), WrongType ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv4), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv4), WrongType ); char *argv5[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -312,8 +313,8 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-nSample", "asdf"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv5), NotEnoughArg ); - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv5), WrongType ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv5), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv5), WrongType ); char *argv6[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", @@ -321,8 +322,8 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-rate", "asdf"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv6), NotEnoughArg ); - CPPUNIT_ASSERT_THROW ( this->input_->core(11, argv6), WrongType ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv6), NotEnoughArg ); + CPPUNIT_ASSERT_THROW ( DEploidIO(11, argv6), WrongType ); } @@ -332,10 +333,10 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv1), FileNameMissing ); - CPPUNIT_ASSERT_THROW_MESSAGE( "\033[1;31mRef count\033[0m file path missing!", this->input_->core(9, argv1), FileNameMissing ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv1), FileNameMissing ); + CPPUNIT_ASSERT_THROW_MESSAGE( "\033[1;31mRef count\033[0m file path missing!", DEploidIO(9, argv1), FileNameMissing ); try { - this->input_->core(9, argv1); + DEploidIO(9, argv1); } catch (const exception &e) { CPPUNIT_ASSERT_EQUAL( string("\033[1;31mRef count\033[0m file path missing!"), string(e.what()) ); @@ -346,10 +347,10 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv2), FileNameMissing ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv2), FileNameMissing ); try { - this->input_->core(9, argv2); + DEploidIO(9, argv2); } catch (const exception &e) { CPPUNIT_ASSERT_EQUAL( string("\033[1;31mAlt count\033[0m file path missing!"), string(e.what()) ); @@ -361,7 +362,7 @@ class TestIO : public CppUnit::TestCase { "-alt", "data/testData/PG0390-C.test.alt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv3), FileNameMissing ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv3), FileNameMissing ); // panel missing char *argv4[] = { "./dEploid", @@ -369,7 +370,7 @@ class TestIO : public CppUnit::TestCase { "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv4), FileNameMissing ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv4), FileNameMissing ); } @@ -380,9 +381,9 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW( this->input_->core(11, argv), InvalidInputFile ); + CPPUNIT_ASSERT_THROW( DEploidIO(11, argv), InvalidInputFile ); try { - this->input_->core(11, argv); + DEploidIO(11, argv); } catch (const exception &e) { CPPUNIT_ASSERT_EQUAL( string("Invalid input file: \033[1;31mPG0390_first100ref.txt\033[0m"), string(e.what()) ); @@ -392,7 +393,7 @@ class TestIO : public CppUnit::TestCase { void testUnknowArg(){ char *argv[] = { "./dEploid", "-unknow"}; - CPPUNIT_ASSERT_THROW( this->input_->core(2, argv), UnknowArg ); + CPPUNIT_ASSERT_THROW( DEploidIO(2, argv), UnknowArg ); } @@ -404,7 +405,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-noPanel" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv1), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv1), FlagsConflict ); // noPanel conflict with panel char *argv2[] = { "./dEploid", @@ -413,7 +414,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel", "-panel", "data/testData/labStrains.test.panel.txt"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(10, argv2), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(10, argv2), FlagsConflict ); // ref, alt conflict with vcf char *argv3[] = { "./dEploid", @@ -423,7 +424,7 @@ class TestIO : public CppUnit::TestCase { "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1", "-vcf", "data/testData/PG0390-C.test.vcf"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(13, argv3), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(13, argv3), FlagsConflict ); // vcf conflict with ref char *argv4[] = { "./dEploid", @@ -433,7 +434,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(13, argv4), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(13, argv4), FlagsConflict ); // vcf conflict with alt char *argv5[] = { "./dEploid", @@ -443,7 +444,7 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(13, argv5), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(13, argv5), FlagsConflict ); // exportPostProb conflict with noPanel char *argv6[] = { "./dEploid", @@ -451,7 +452,7 @@ class TestIO : public CppUnit::TestCase { "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel", "-exportPostProb" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv6), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv6), FlagsConflict ); // noPanel conflict with exportPostProb char *argv7[] = { "./dEploid", @@ -459,7 +460,7 @@ class TestIO : public CppUnit::TestCase { "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-exportPostProb", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv7), FlagsConflict ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv7), FlagsConflict ); } @@ -469,14 +470,14 @@ class TestIO : public CppUnit::TestCase { "-alt", "data/testData/PG0390.C_alt.txt", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv1), LociNumberUnequal ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv1), LociNumberUnequal ); char *argv2[] = { "./dEploid", "-ref", "data/testData/PG0390-C.test.ref", "-alt", "data/testData/PG0390-C.test.alt", "-plaf", "data/testData/labStrains_samples_PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt" }; - CPPUNIT_ASSERT_THROW ( this->input_->core(9, argv2), LociNumberUnequal ); + CPPUNIT_ASSERT_THROW ( DEploidIO(9, argv2), LociNumberUnequal ); } @@ -488,36 +489,38 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(11, argv1) ); - CPPUNIT_ASSERT_EQUAL((size_t)594, this->input_->nLoci_); + + CPPUNIT_ASSERT_NO_THROW (DEploidIO(11, argv1)); + DEploidIO dEploidIO1(11, argv1); + CPPUNIT_ASSERT_EQUAL((size_t)594, dEploidIO1.nLoci_); //Pf3D7_01_v3 93157 0.35824742268032 - CPPUNIT_ASSERT_DOUBLES_EQUAL(0.35824742268032, this->input_->plaf_[0], 0.000000000001); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.35824742268032, dEploidIO1.plaf_[0], 0.000000000001); //Pf3D7_01_v3 93157 85 - CPPUNIT_ASSERT_EQUAL(85.0, this->input_->refCount_[0]); + CPPUNIT_ASSERT_EQUAL(85.0, dEploidIO1.refCount_[0]); //Pf3D7_01_v3 93157 0 - CPPUNIT_ASSERT_EQUAL(0.0, this->input_->altCount_[0]); + CPPUNIT_ASSERT_EQUAL(0.0, dEploidIO1.altCount_[0]); //Pf3D7_01_v3 95518 0.553038882597102 - CPPUNIT_ASSERT_DOUBLES_EQUAL(0.553038882597102, this->input_->plaf_[4], 0.000000000001); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.553038882597102, dEploidIO1.plaf_[4], 0.000000000001); //Pf3D7_01_v3 95518 156 - CPPUNIT_ASSERT_EQUAL(156.0, this->input_->refCount_[4]); + CPPUNIT_ASSERT_EQUAL(156.0, dEploidIO1.refCount_[4]); //Pf3D7_01_v3 95518 46 - CPPUNIT_ASSERT_EQUAL(46.0, this->input_->altCount_[4]); + CPPUNIT_ASSERT_EQUAL(46.0, dEploidIO1.altCount_[4]); //Pf3D7_01_v3 112038 0.178471474703944 - CPPUNIT_ASSERT_DOUBLES_EQUAL(0.178471474703944, this->input_->plaf_[20], 0.000000000001); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.178471474703944, dEploidIO1.plaf_[20], 0.000000000001); //Pf3D7_01_v3 112038 157 - CPPUNIT_ASSERT_EQUAL(157.0, this->input_->refCount_[20]); + CPPUNIT_ASSERT_EQUAL(157.0, dEploidIO1.refCount_[20]); //Pf3D7_01_v3 112038 28 - CPPUNIT_ASSERT_EQUAL(28.0, this->input_->altCount_[20]); + CPPUNIT_ASSERT_EQUAL(28.0, dEploidIO1.altCount_[20]); //Pf3D7_01_v3 180192 0.747258396161626 - CPPUNIT_ASSERT_DOUBLES_EQUAL(0.747258396161626, this->input_->plaf_[99], 0.000000000001); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.747258396161626, dEploidIO1.plaf_[99], 0.000000000001); //Pf3D7_01_v3 180192 199 - CPPUNIT_ASSERT_EQUAL(199.0, this->input_->refCount_[99]); + CPPUNIT_ASSERT_EQUAL(199.0, dEploidIO1.refCount_[99]); //Pf3D7_01_v3 180192 0 - CPPUNIT_ASSERT_EQUAL(0.0, this->input_->altCount_[99]); + CPPUNIT_ASSERT_EQUAL(0.0, dEploidIO1.altCount_[99]); } @@ -528,25 +531,29 @@ class TestIO : public CppUnit::TestCase { "-plaf", "data/testData/labStrains.test.PLAF.txt", "-panel", "data/testData/labStrains.test.panel.txt", "-o", "tmp1", "-forbidUpdateProp", "-forbidUpdateSingle", "-forbidUpdatePair" }; - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(14, argv1) ); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateProp(), false); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateSingle(), false); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdatePair(), false); - - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(13, argv1) ); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateProp(), false); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateSingle(), false); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdatePair(), true); - - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(12, argv1) ); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateProp(), false); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateSingle(), true); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdatePair(), true); - - CPPUNIT_ASSERT_NO_THROW ( this->input_->core(11, argv1) ); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateProp() , true); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdateSingle() , true); - CPPUNIT_ASSERT_EQUAL( this->input_->doUpdatePair() , true); + CPPUNIT_ASSERT_NO_THROW (DEploidIO(14, argv1)); + DEploidIO dEploidIO1(14, argv1); + CPPUNIT_ASSERT_EQUAL( dEploidIO1.doUpdateProp(), false); + CPPUNIT_ASSERT_EQUAL( dEploidIO1.doUpdateSingle(), false); + CPPUNIT_ASSERT_EQUAL( dEploidIO1.doUpdatePair(), false); + + CPPUNIT_ASSERT_NO_THROW (DEploidIO(13, argv1)); + DEploidIO dEploidIO2(13, argv1); + CPPUNIT_ASSERT_EQUAL( dEploidIO2.doUpdateProp(), false); + CPPUNIT_ASSERT_EQUAL( dEploidIO2.doUpdateSingle(), false); + CPPUNIT_ASSERT_EQUAL( dEploidIO2.doUpdatePair(), true); + + CPPUNIT_ASSERT_NO_THROW (DEploidIO(12, argv1)); + DEploidIO dEploidIO3(12, argv1); + CPPUNIT_ASSERT_EQUAL( dEploidIO3.doUpdateProp(), false); + CPPUNIT_ASSERT_EQUAL( dEploidIO3.doUpdateSingle(), true); + CPPUNIT_ASSERT_EQUAL( dEploidIO3.doUpdatePair(), true); + + CPPUNIT_ASSERT_NO_THROW (DEploidIO(11, argv1)); + DEploidIO dEploidIO4(11, argv1); + CPPUNIT_ASSERT_EQUAL( dEploidIO4.doUpdateProp() , true); + CPPUNIT_ASSERT_EQUAL( dEploidIO4.doUpdateSingle() , true); + CPPUNIT_ASSERT_EQUAL( dEploidIO4.doUpdatePair() , true); } @@ -554,7 +561,7 @@ class TestIO : public CppUnit::TestCase { char *argv1[] = { "./dEploid", "-vcf", "data/testData/PG0389-C.test.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv1), InvalidInputFile ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv1), InvalidInputFile ); } @@ -562,66 +569,66 @@ class TestIO : public CppUnit::TestCase { char *argv1[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.alt.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv1), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv1), VcfInvalidHeaderFieldNames ); char *argv2[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.chrom2.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv2), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv2), VcfInvalidHeaderFieldNames ); char *argv3[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.chrom.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv3), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv3), VcfInvalidHeaderFieldNames ); char *argv4[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.filter.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv4), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv4), VcfInvalidHeaderFieldNames ); char *argv5[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.format.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv5), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv5), VcfInvalidHeaderFieldNames ); char *argv6[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.id.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv6), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv6), VcfInvalidHeaderFieldNames ); char *argv7[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.info.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv7), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv7), VcfInvalidHeaderFieldNames ); char *argv8[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.pos.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv8), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv8), VcfInvalidHeaderFieldNames ); char *argv9[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.qual.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv9), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv9), VcfInvalidHeaderFieldNames ); char *argv10[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badHeaderFieldNames.ref.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv10), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv10), VcfInvalidHeaderFieldNames ); } void testVcfNoAD(){ char *argv1[] = { "./dEploid", "-vcf", "data/testData/crappyVcf/badVariant.noAD.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv1), VcfCoverageFieldNotFound ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv1), VcfCoverageFieldNotFound ); } void testInvalidVcfGz(){ char *argv1[] = { "./dEploid", "-vcf", "data/testData/PG0389-C.test.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv1), InvalidInputFile ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv1), InvalidInputFile ); } @@ -629,66 +636,66 @@ class TestIO : public CppUnit::TestCase { char *argv1[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.alt.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv1), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv1), VcfInvalidHeaderFieldNames ); char *argv2[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.chrom2.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv2), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv2), VcfInvalidHeaderFieldNames ); char *argv3[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.chrom.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv3), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv3), VcfInvalidHeaderFieldNames ); char *argv4[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.filter.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv4), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv4), VcfInvalidHeaderFieldNames ); char *argv5[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.format.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv5), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv5), VcfInvalidHeaderFieldNames ); char *argv6[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.id.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv6), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv6), VcfInvalidHeaderFieldNames ); char *argv7[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.info.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv7), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv7), VcfInvalidHeaderFieldNames ); char *argv8[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.pos.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv8), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv8), VcfInvalidHeaderFieldNames ); char *argv9[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.qual.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv9), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv9), VcfInvalidHeaderFieldNames ); char *argv10[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badHeaderFieldNames.ref.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv10), VcfInvalidHeaderFieldNames ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv10), VcfInvalidHeaderFieldNames ); } void testVcfGzNoAD(){ char *argv1[] = { "./dEploid", "-vcf", "data/testData/crappyVcfGz/badVariant.noAD.vcf.gz", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(6, argv1), VcfCoverageFieldNotFound ); + CPPUNIT_ASSERT_THROW ( DEploidIO(6, argv1), VcfCoverageFieldNotFound ); } void testVcfOutUnSpecified(){ char *argv1[] = { "./dEploid", "-vcf", "data/testData/PG0390-C.test.vcf", "-plaf", "data/testData/labStrains.test.PLAF.txt", "-noPanel", "-z"}; - CPPUNIT_ASSERT_THROW ( this->input_->core(7, argv1), VcfOutUnSpecified ); + CPPUNIT_ASSERT_THROW ( DEploidIO(7, argv1), VcfOutUnSpecified ); } }; From 8e624a58019af0bc7185328c6eb6e95ec5c0f4bf Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Tue, 4 Oct 2016 02:10:37 +0100 Subject: [PATCH 03/29] progress --- src/dEploid.cpp | 5 ++- src/dEploidIO.cpp | 70 +++++++++++++++---------------- src/dEploidIO.hpp | 4 +- tests/unittest/test_dEploidIO.cpp | 12 +++--- 4 files changed, 47 insertions(+), 44 deletions(-) diff --git a/src/dEploid.cpp b/src/dEploid.cpp index 887eb800..57b17d0b 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -31,14 +31,15 @@ using namespace std; int main( int argc, char *argv[] ){ try { DEploidIO dEploidIO(argc, argv); + std::ostream *output = &std::cout; if ( dEploidIO.version() ){ - dEploidIO.printVersion(); + dEploidIO.printVersion(*output); return EXIT_SUCCESS; } if ( dEploidIO.help() ){ - dEploidIO.printHelp(); + dEploidIO.printHelp(*output); return EXIT_SUCCESS; } diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index 7dc5bdad..3e68a68d 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -410,43 +410,43 @@ void DEploidIO::readNextStringto( string &readto ){ } -void DEploidIO::printVersion(){ - cout << endl - << "dEploid " << VERSION - << endl - << "Git commit: " << dEploidGitVersion_ << endl; +void DEploidIO::printVersion(std::ostream& out){ + out << endl + << "dEploid " << VERSION + << endl + << "Git commit: " << dEploidGitVersion_ << endl; } -void DEploidIO::printHelp(){ - cout << endl - << "dEploid " << VERSION - << endl - << endl; - cout << "Usage:" - << endl; - cout << setw(20) << "-h or -help" << " -- " << "Help. List the following content."< Date: Tue, 4 Oct 2016 13:34:25 +0100 Subject: [PATCH 04/29] set argument to the first value --- src/dEploidIO.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index 3e68a68d..eaeef839 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -40,6 +40,7 @@ DEploidIO::DEploidIO(const std::string &arg) { copy(std::istream_iterator(iss), std::istream_iterator(), std::back_inserter(argv_)); + this->argv_i = argv_.begin(); this->core(); } From e5e9517c106c5323f3356333703cbaef74a1e312 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 13:32:00 +0100 Subject: [PATCH 05/29] update copy right --- src/dEploid.cpp | 6 +++--- src/dEploidIO.cpp | 3 +-- src/dEploidIO.hpp | 6 +++--- src/debug/mcmcDebug.cpp | 6 +++--- src/debug/vcfReaderDebug.cpp | 6 +++--- src/exceptions.hpp | 6 +++--- src/export/dEploidIOExport.cpp | 6 +++--- src/export/dEploidIOExportPosteriorProb.cpp | 6 +++--- src/mcmc.cpp | 10 +++++----- src/mcmc.hpp | 8 ++++---- src/panel.cpp | 6 +++--- src/panel.hpp | 6 +++--- src/random/fastfunc.cpp | 6 +++--- src/random/fastfunc.hpp | 10 +++++----- src/random/mersenne_twister.cpp | 6 +++--- src/random/mersenne_twister.hpp | 10 +++++----- src/random/random_generator.cpp | 6 +++--- src/random/random_generator.hpp | 12 ++++++------ src/txtReader.cpp | 6 +++--- src/txtReader.hpp | 6 +++--- src/updateHap.cpp | 6 +++--- src/updateHap.hpp | 6 +++--- src/utility.cpp | 6 +++--- src/utility.hpp | 6 +++--- src/variantIndex.cpp | 6 +++--- src/variantIndex.hpp | 6 +++--- src/vcfReader.cpp | 6 +++--- src/vcfReader.hpp | 6 +++--- utilities/dEploidTools.r | 19 ++++++++++++++++++- 29 files changed, 110 insertions(+), 94 deletions(-) diff --git a/src/dEploid.cpp b/src/dEploid.cpp index 57b17d0b..8526d88d 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include // std::cout #include "mcmc.hpp" diff --git a/src/dEploidIO.cpp b/src/dEploidIO.cpp index eaeef839..ef5b7f4d 100644 --- a/src/dEploidIO.cpp +++ b/src/dEploidIO.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,7 +18,6 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - */ #include "dEploidIO.hpp" diff --git a/src/dEploidIO.hpp b/src/dEploidIO.hpp index ba23b37c..69f93d91 100644 --- a/src/dEploidIO.hpp +++ b/src/dEploidIO.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include /* strtol, strtod */ #include diff --git a/src/debug/mcmcDebug.cpp b/src/debug/mcmcDebug.cpp index a159909f..2a6fb6f3 100644 --- a/src/debug/mcmcDebug.cpp +++ b/src/debug/mcmcDebug.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "mcmc.hpp" //#include "utility.hpp" diff --git a/src/debug/vcfReaderDebug.cpp b/src/debug/vcfReaderDebug.cpp index 87203417..8770abc8 100644 --- a/src/debug/vcfReaderDebug.cpp +++ b/src/debug/vcfReaderDebug.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "global.h" #include // std::cout diff --git a/src/exceptions.hpp b/src/exceptions.hpp index 69f8a6fc..8de98464 100644 --- a/src/exceptions.hpp +++ b/src/exceptions.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include /* string */ diff --git a/src/export/dEploidIOExport.cpp b/src/export/dEploidIOExport.cpp index 52b1db17..9bf94ef4 100644 --- a/src/export/dEploidIOExport.cpp +++ b/src/export/dEploidIOExport.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "dEploidIO.hpp" #include "mcmc.hpp" diff --git a/src/export/dEploidIOExportPosteriorProb.cpp b/src/export/dEploidIOExportPosteriorProb.cpp index 6d502593..48dd91f7 100644 --- a/src/export/dEploidIOExportPosteriorProb.cpp +++ b/src/export/dEploidIOExportPosteriorProb.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "dEploidIO.hpp" #include "updateHap.hpp" diff --git a/src/mcmc.cpp b/src/mcmc.cpp index 6afc5793..f83eff76 100644 --- a/src/mcmc.cpp +++ b/src/mcmc.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "mcmc.hpp" #include "utility.hpp" @@ -212,11 +212,11 @@ vector McmcMachinery::titre2prop(vector & tmpTitre){ } -void McmcMachinery::runMcmcChain( ){ +void McmcMachinery::runMcmcChain( bool showProgress ){ for ( this->currentMcmcIteration_ = 0 ; currentMcmcIteration_ < this->maxIteration_ ; currentMcmcIteration_++){ dout << endl; dout << "MCMC iteration: " << this->currentMcmcIteration_ << endl; - if ( this->currentMcmcIteration_ > 0 && this->currentMcmcIteration_%100 == 0 ){ + if ( this->currentMcmcIteration_ > 0 && this->currentMcmcIteration_%100 == 0 && showProgress ){ clog << "\r" << " MCMC step" << setw(4) << int(currentMcmcIteration_ * 100 / this->maxIteration_) << "% completed."<sampleMcmcEvent(); diff --git a/src/mcmc.hpp b/src/mcmc.hpp index 8335cf23..2af5dfc4 100644 --- a/src/mcmc.hpp +++ b/src/mcmc.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include #include @@ -68,7 +68,7 @@ class McmcMachinery { //McmcMachinery(); McmcMachinery( DEploidIO* pdfDeconfIO, Panel *panel, McmcSample *mcmcSample ); ~McmcMachinery(); - void runMcmcChain( ); + void runMcmcChain( bool showProgress = true ); private: McmcSample* mcmcSample_; diff --git a/src/panel.cpp b/src/panel.cpp index 0e30b986..4ded441b 100644 --- a/src/panel.cpp +++ b/src/panel.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "panel.hpp" #include diff --git a/src/panel.hpp b/src/panel.hpp index fd222a94..0a35d136 100644 --- a/src/panel.hpp +++ b/src/panel.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #ifndef PANEL #define PANEL diff --git a/src/random/fastfunc.cpp b/src/random/fastfunc.cpp index f294f12d..2c22ec1c 100644 --- a/src/random/fastfunc.cpp +++ b/src/random/fastfunc.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "fastfunc.hpp" diff --git a/src/random/fastfunc.hpp b/src/random/fastfunc.hpp index 94eb74bc..15fa3da6 100644 --- a/src/random/fastfunc.hpp +++ b/src/random/fastfunc.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,11 +18,11 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . + * + */ -*/ - -#ifndef scrm_src_random_fastfunc -#define scrm_src_random_fastfunc +#ifndef dEploid_src_random_fastfunc +#define dEploid_src_random_fastfunc #include #include diff --git a/src/random/mersenne_twister.cpp b/src/random/mersenne_twister.cpp index 2f0f09aa..ce9a1c92 100644 --- a/src/random/mersenne_twister.cpp +++ b/src/random/mersenne_twister.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "mersenne_twister.hpp" diff --git a/src/random/mersenne_twister.hpp b/src/random/mersenne_twister.hpp index c5374ceb..cacdd3b5 100644 --- a/src/random/mersenne_twister.hpp +++ b/src/random/mersenne_twister.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,11 +18,11 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . + * + */ -*/ - -#ifndef scrm_src_random_mersenne_twister -#define scrm_src_random_mersenne_twister +#ifndef dEploid_src_random_mersenne_twister +#define dEploid_src_random_mersenne_twister #include #include "random_generator.hpp" diff --git a/src/random/random_generator.cpp b/src/random/random_generator.cpp index 7d683fb9..32e2e83e 100644 --- a/src/random/random_generator.cpp +++ b/src/random/random_generator.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "random_generator.hpp" #include diff --git a/src/random/random_generator.hpp b/src/random/random_generator.hpp index 1b67fb0d..c492b710 100644 --- a/src/random/random_generator.hpp +++ b/src/random/random_generator.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,11 +18,11 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . + * + */ -*/ - -#ifndef scrm_src_random_random_generator -#define scrm_src_random_random_generator +#ifndef dEploid_src_random_random_generator +#define dEploid_src_random_random_generator //#include "../macros.h" // Needs to be before cassert @@ -88,7 +88,7 @@ class RandomGenerator protected: // Sample from a unit exponential distribution // Unit tested - // fastlog can return 0, which causes a bug in scrm. + // fastlog can return 0, which causes a bug in dEploid. // log or fastlog does seems to have an influence on the runtime. virtual double sampleUnitExponential() { return -(this->ff()->fastlog(sample())); diff --git a/src/txtReader.cpp b/src/txtReader.cpp index a8a51b79..53fbc8be 100644 --- a/src/txtReader.cpp +++ b/src/txtReader.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include #include diff --git a/src/txtReader.hpp b/src/txtReader.hpp index 23d6ea7d..43e0b4b9 100644 --- a/src/txtReader.hpp +++ b/src/txtReader.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #ifndef TXTREADER diff --git a/src/updateHap.cpp b/src/updateHap.cpp index be1d51d4..2d451180 100644 --- a/src/updateHap.cpp +++ b/src/updateHap.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "updateHap.hpp" #include // std::reverse diff --git a/src/updateHap.hpp b/src/updateHap.hpp index df23b855..949599c3 100644 --- a/src/updateHap.hpp +++ b/src/updateHap.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include #include diff --git a/src/utility.cpp b/src/utility.cpp index 3f135f5b..77d2737b 100644 --- a/src/utility.cpp +++ b/src/utility.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "logbeta.h" #include "utility.hpp" diff --git a/src/utility.hpp b/src/utility.hpp index f6669a82..8d055cf1 100644 --- a/src/utility.hpp +++ b/src/utility.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include #include diff --git a/src/variantIndex.cpp b/src/variantIndex.cpp index d43646c6..a35d0680 100644 --- a/src/variantIndex.cpp +++ b/src/variantIndex.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "txtReader.hpp" diff --git a/src/variantIndex.hpp b/src/variantIndex.hpp index 378c0e99..7384dde5 100644 --- a/src/variantIndex.hpp +++ b/src/variantIndex.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #ifndef VARIANTINDEX #define VARIANTINDEX diff --git a/src/vcfReader.cpp b/src/vcfReader.cpp index 9d3ea296..cdc98005 100644 --- a/src/vcfReader.cpp +++ b/src/vcfReader.cpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include "global.h" #include // assert diff --git a/src/vcfReader.hpp b/src/vcfReader.hpp index 2b19f810..a897287a 100644 --- a/src/vcfReader.hpp +++ b/src/vcfReader.hpp @@ -6,7 +6,7 @@ * * This file is part of dEploid. * - * scrm is free software: you can redistribute it and/or modify + * dEploid is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. @@ -18,8 +18,8 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see . - -*/ + * + */ #include /* string */ #include /* vector */ diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index a7f32b04..5cae7c1c 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -368,7 +368,24 @@ fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclu } - +#' Deconvolute mixed haplotypes, and reporting the mixture proportions from each sample +#' This function provieds an interface for calling \emph{dEploid} from R. +#' The command line options are passed via the \code{args} argument and \code{file} +#' +#' @section blahblah: +#' Blah blah +#' +#' @param tmpProp blah blah blah +#' +#' @param title blah blah blah +#' +#' @return A named list of something something ... +#' +#' @export +#' +#' @examples +#' set.seed(1234) +#' plot.prob <-function (tmpProp, title){ rainbowColorBin = 16 barplot(t(tmpProp), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="SNP index", ylab="Posterior probabilities", main=title) From 6579b4fdd07f9f54c50fcdac4e001e60d63258d9 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 13:51:02 +0100 Subject: [PATCH 06/29] update --- src/dEploid.cpp | 2 +- utilities/dataExplore.r | 2 +- utilities/interpretDEploid.r | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dEploid.cpp b/src/dEploid.cpp index 8526d88d..d4624b96 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -59,7 +59,7 @@ int main( int argc, char *argv[] ){ McmcSample * mcmcSample = new McmcSample(); McmcMachinery mcmcMachinery(&dEploidIO, panel, mcmcSample); - mcmcMachinery.runMcmcChain(); + mcmcMachinery.runMcmcChain(true); dEploidIO.write(mcmcSample, panel); diff --git a/utilities/dataExplore.r b/utilities/dataExplore.r index 27feae94..6a09bc05 100644 --- a/utilities/dataExplore.r +++ b/utilities/dataExplore.r @@ -1,4 +1,4 @@ -rm(list=ls()); dEploidRootDir="/home/joezhu/DEploid" +rm(list=ls()); dEploidRootDir="/home/joezhu/DEploid-r/.DEploid" # DESCRIPTION: # # USAGE: diff --git a/utilities/interpretDEploid.r b/utilities/interpretDEploid.r index 51d8ab5c..f7587935 100644 --- a/utilities/interpretDEploid.r +++ b/utilities/interpretDEploid.r @@ -1,4 +1,4 @@ -rm(list=ls()); dEploidRootDir="/home/joezhu/DEploid" +rm(list=ls()); dEploidRootDir="/home/joezhu/DEploid-r/.DEploid" # DESCRIPTION: # # USAGE: From 05374bc1df4b9b81f9b02de4fede98745affbd4d Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 14:31:22 +0100 Subject: [PATCH 07/29] use valgrind-devel instead of valgrind on mac --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 40a22f93..77b458a4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,7 +11,7 @@ before_install: - echo $LC_ALL - if [ $TRAVIS_OS_NAME == linux ]; then sudo apt-get update -qq; sudo apt-get install -qq libcppunit-dev valgrind r-base-core lcov python-pip doxygen graphviz; pip install --user cpp-coveralls; fi - if [ $TRAVIS_OS_NAME == linux ]; then apt-cache policy zlib*; fi - - if [ $TRAVIS_OS_NAME == osx ]; then brew tap homebrew/dupes && brew update && brew install llvm cppunit valgrind zlib && brew link --force llvm; fi + - if [ $TRAVIS_OS_NAME == osx ]; then brew tap homebrew/dupes && brew update && brew install llvm cppunit valgrind-devel zlib && brew link --force llvm; fi before_script: - ./bootstrap From a0d59cc2a33d990e84b51bee19e4b972aed2cc01 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 14:35:50 +0100 Subject: [PATCH 08/29] valgrind-devel don't work --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 77b458a4..40a22f93 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,7 +11,7 @@ before_install: - echo $LC_ALL - if [ $TRAVIS_OS_NAME == linux ]; then sudo apt-get update -qq; sudo apt-get install -qq libcppunit-dev valgrind r-base-core lcov python-pip doxygen graphviz; pip install --user cpp-coveralls; fi - if [ $TRAVIS_OS_NAME == linux ]; then apt-cache policy zlib*; fi - - if [ $TRAVIS_OS_NAME == osx ]; then brew tap homebrew/dupes && brew update && brew install llvm cppunit valgrind-devel zlib && brew link --force llvm; fi + - if [ $TRAVIS_OS_NAME == osx ]; then brew tap homebrew/dupes && brew update && brew install llvm cppunit valgrind zlib && brew link --force llvm; fi before_script: - ./bootstrap From fb987086564a4dd7d9145fcd74712eaf2c6d533d Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 14:45:07 +0100 Subject: [PATCH 09/29] valgrind suppress on --- tests/test_binary.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_binary.sh b/tests/test_binary.sh index 83b55bef..3cdf45b8 100755 --- a/tests/test_binary.sh +++ b/tests/test_binary.sh @@ -16,7 +16,7 @@ function test_dEploid { echo -n "." # Test for memory leaks - valgrind --error-exitcode=1 --leak-check=full -q ./dEploid $@ -seed $i > /dev/null + valgrind --error-exitcode=1 --leak-check=full -q --gen-suppressions=yes ./dEploid $@ -seed $i > /dev/null if [ $? -ne 0 ]; then echo "" echo "Valgrind check of \"./dEploid $@ -seed $i\" failed." @@ -39,7 +39,7 @@ function test_noRepeat { fi # Test for memory leaks - valgrind --error-exitcode=1 --leak-check=full -q ./dEploid $@ > /dev/null + valgrind --error-exitcode=1 --leak-check=full -q --gen-suppressions=yes ./dEploid $@ > /dev/null if [ $? -ne 0 ]; then echo "" echo "Valgrind check of \"./dEploid $@ \" failed." From 7d34f760470bc3c331e7e1277eb8380595a30ebb Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 16:46:14 +0100 Subject: [PATCH 10/29] rework on random number generator, so r version will work --- src/dEploid.cpp | 5 ++++- src/dEploidIO.hpp | 2 +- src/mcmc.cpp | 14 +++++++------- src/mcmc.hpp | 10 +++++----- src/random/random_generator.hpp | 1 + src/updateHap.cpp | 6 +++--- src/updateHap.hpp | 12 ++++++------ src/utility.cpp | 2 +- src/utility.hpp | 2 +- tests/unittest/test_mcmc.cpp | 8 ++++++-- 10 files changed, 35 insertions(+), 27 deletions(-) diff --git a/src/dEploid.cpp b/src/dEploid.cpp index d4624b96..977fa92d 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -57,11 +57,14 @@ int main( int argc, char *argv[] ){ } McmcSample * mcmcSample = new McmcSample(); + MersenneTwister rg(dEploidIO.randomSeed()); - McmcMachinery mcmcMachinery(&dEploidIO, panel, mcmcSample); + McmcMachinery mcmcMachinery(&dEploidIO, panel, mcmcSample, &rg); mcmcMachinery.runMcmcChain(true); dEploidIO.write(mcmcSample, panel); + //rg->clearFastFunc(); + //delete rg; if ( panel ){ delete panel; diff --git a/src/dEploidIO.hpp b/src/dEploidIO.hpp index 69f93d91..bba9afb5 100644 --- a/src/dEploidIO.hpp +++ b/src/dEploidIO.hpp @@ -85,6 +85,7 @@ class DEploidIO{ bool randomSeedWasSet() const {return this->randomSeedWasSet_; } friend std::ostream& operator<< (std::ostream& stream, const DEploidIO& dEploidIO); + size_t randomSeed() const { return randomSeed_;} private: @@ -96,7 +97,6 @@ class DEploidIO{ string excludeFileName_; string prefix_; size_t randomSeed_; - size_t randomSeed() const { return randomSeed_;} bool randomSeedWasSet_; void setRandomSeedWasSet(const bool random){ this->randomSeedWasSet_ = random; } diff --git a/src/mcmc.cpp b/src/mcmc.cpp index f83eff76..820c7a38 100644 --- a/src/mcmc.cpp +++ b/src/mcmc.cpp @@ -31,15 +31,15 @@ McmcSample::McmcSample(){}; McmcSample::~McmcSample(){}; -McmcMachinery::McmcMachinery(DEploidIO* pfDeconfIO, Panel *panel, McmcSample *mcmcSample ){ // initialiseMCMCmachinery +McmcMachinery::McmcMachinery(DEploidIO* pfDeconfIO, Panel *panel, McmcSample *mcmcSample, RandomGenerator* rg_ ){ // initialiseMCMCmachinery this->dEploidIO_ = pfDeconfIO; this->panel_ = panel; this->mcmcSample_ = mcmcSample; - this->seed_ = this->dEploidIO_->randomSeed(); - this->hapRg_ = new MersenneTwister(this->seed_); + //this->hapRg_ = new MersenneTwister(this->seed_); + this->hapRg_ = rg_ ; this->mcmcEventRg_ = this->hapRg_; this->propRg_ = this->hapRg_; this->initialHapRg_ = this->hapRg_; @@ -65,10 +65,10 @@ McmcMachinery::McmcMachinery(DEploidIO* pfDeconfIO, Panel *panel, McmcSample *mc McmcMachinery::~McmcMachinery(){ - if ( this->hapRg_ ){ - this->hapRg_->clearFastFunc(); - delete hapRg_; - } + //if ( this->hapRg_ ){ + //this->hapRg_->clearFastFunc(); + //delete hapRg_; + //} //if ( this->mcmcEventRg_ ){ //this->mcmcEventRg_->clearFastFunc(); diff --git a/src/mcmc.hpp b/src/mcmc.hpp index 2af5dfc4..d22e6741 100644 --- a/src/mcmc.hpp +++ b/src/mcmc.hpp @@ -66,7 +66,7 @@ class McmcMachinery { #endif public: //McmcMachinery(); - McmcMachinery( DEploidIO* pdfDeconfIO, Panel *panel, McmcSample *mcmcSample ); + McmcMachinery( DEploidIO* pdfDeconfIO, Panel *panel, McmcSample *mcmcSample, RandomGenerator* rg_ ); ~McmcMachinery(); void runMcmcChain( bool showProgress = true ); @@ -89,10 +89,10 @@ class McmcMachinery { size_t strainIndex2_; size_t seed_; - MersenneTwister* hapRg_; - MersenneTwister* mcmcEventRg_; - MersenneTwister* propRg_; - MersenneTwister* initialHapRg_; + RandomGenerator* hapRg_; + RandomGenerator* mcmcEventRg_; + RandomGenerator* propRg_; + RandomGenerator* initialHapRg_; std::default_random_engine* std_generator_;// (this->seed_); std::normal_distribution* initialTitre_normal_distribution_;// (MN_LOG_TITRE, SD_LOG_TITRE); diff --git a/src/random/random_generator.hpp b/src/random/random_generator.hpp index c492b710..408de933 100644 --- a/src/random/random_generator.hpp +++ b/src/random/random_generator.hpp @@ -35,6 +35,7 @@ class RandomGenerator { friend class MersenneTwister; + friend class RRandomGenerator; public: RandomGenerator() { this->ff_ = new FastFunc(); }; RandomGenerator( FastFunc* ff ) { this->ff_ = ff; }; diff --git a/src/updateHap.cpp b/src/updateHap.cpp index 2d451180..642590ce 100644 --- a/src/updateHap.cpp +++ b/src/updateHap.cpp @@ -33,7 +33,7 @@ UpdateHap::UpdateHap( vector &refCount, vector &expectedWsaf, vector &proportion, vector < vector > &haplotypes, - MersenneTwister* rg, + RandomGenerator* rg, size_t segmentStartIndex, size_t nLoci, Panel* panel, @@ -86,7 +86,7 @@ UpdateSingleHap::UpdateSingleHap( vector &refCount, vector &expectedWsaf, vector &proportion, vector < vector > &haplotypes, - MersenneTwister* rg, + RandomGenerator* rg, size_t segmentStartIndex, size_t nLoci, Panel* panel, double missCopyProb, @@ -307,7 +307,7 @@ UpdatePairHap::UpdatePairHap( vector &refCount, vector &expectedWsaf, vector &proportion, vector < vector > &haplotypes, - MersenneTwister* rg, + RandomGenerator* rg, size_t segmentStartIndex, size_t nLoci, Panel* panel, double missCopyProb, bool forbidCopyFromSame, diff --git a/src/updateHap.hpp b/src/updateHap.hpp index 949599c3..9afc33a6 100644 --- a/src/updateHap.hpp +++ b/src/updateHap.hpp @@ -51,7 +51,7 @@ class UpdateHap{ vector &expectedWsaf, vector &proportion, vector < vector > &haplotypes, - MersenneTwister* rg, + RandomGenerator* rg, size_t segmentStartIndex, size_t nLoci, Panel* panel, @@ -60,9 +60,9 @@ class UpdateHap{ Panel* panel_; double missCopyProb_; - MersenneTwister* recombRg_; - MersenneTwister* recombLevel2Rg_; - MersenneTwister* missCopyRg_; + RandomGenerator* recombRg_; + RandomGenerator* recombLevel2Rg_; + RandomGenerator* missCopyRg_; size_t kStrain_; size_t nPanel_; @@ -106,7 +106,7 @@ class UpdateSingleHap : public UpdateHap{ vector &expectedWsaf, vector &proportion, vector < vector > &haplotypes, - MersenneTwister* rg, + RandomGenerator* rg, size_t segmentStartIndex, size_t nLoci, Panel* panel, double missCopyProb, @@ -160,7 +160,7 @@ class UpdatePairHap : public UpdateHap{ vector &expectedWsaf, vector &proportion, vector < vector > &haplotypes, - MersenneTwister* rg, + RandomGenerator* rg, size_t segmentStartIndex, size_t nLoci, Panel* panel, double missCopyProb, bool forbidCopyFromSame, diff --git a/src/utility.cpp b/src/utility.cpp index 77d2737b..ce4e35d1 100644 --- a/src/utility.cpp +++ b/src/utility.cpp @@ -123,7 +123,7 @@ double calcLLK( double ref, double alt, double unadjustedWsaf, double err, doubl } -size_t sampleIndexGivenProp ( MersenneTwister* rg, vector proportion ){ +size_t sampleIndexGivenProp ( RandomGenerator* rg, vector proportion ){ #ifndef NDEBUG auto biggest = std::max_element(std::begin(proportion), std::end(proportion)); return std::distance(proportion.begin(), biggest); diff --git a/src/utility.hpp b/src/utility.hpp index 8d055cf1..adc680fc 100644 --- a/src/utility.hpp +++ b/src/utility.hpp @@ -76,7 +76,7 @@ void normalizeBySum ( vector & array ); void normalizeBySumMat ( vector > & matrix ); vector calcLLKs( vector &refCount, vector &altCount, vector &expectedWsaf, size_t firstIndex, size_t length ); double calcLLK( double ref, double alt, double unadjustedWsaf, double err = 0.01, double fac=100 ) ; -size_t sampleIndexGivenProp ( MersenneTwister* rg, vector proportion ); +size_t sampleIndexGivenProp ( RandomGenerator* rg, vector proportion ); vector reshapeMatToVec ( vector < vector > &Mat ); diff --git a/tests/unittest/test_mcmc.cpp b/tests/unittest/test_mcmc.cpp index abf59562..a3ca0267 100644 --- a/tests/unittest/test_mcmc.cpp +++ b/tests/unittest/test_mcmc.cpp @@ -51,6 +51,7 @@ class TestMcmcMachinery: public CppUnit::TestCase { DEploidIO* dEploidIO_; Panel* panel_; McmcMachinery* mcmcMachinery_; + MersenneTwister* rg_; size_t nRepeat; double epsilon1; double epsilon2; @@ -95,8 +96,9 @@ class TestMcmcMachinery: public CppUnit::TestCase { void setUp() { mcmcSample_ = new McmcSample(); dEploidIO_ = new DEploidIO(); + rg_ = new MersenneTwister(dEploidIO_->randomSeed()); panel_ = new Panel(); - mcmcMachinery_ = new McmcMachinery(this->dEploidIO_, this->panel_, this->mcmcSample_ ); + mcmcMachinery_ = new McmcMachinery(this->dEploidIO_, this->panel_, this->mcmcSample_, this->rg_ ); nRepeat = 1000000; epsilon1 = 0.01; epsilon2 = 0.001; @@ -105,6 +107,8 @@ class TestMcmcMachinery: public CppUnit::TestCase { void tearDown() { + rg_->clearFastFunc(); + delete rg_; delete mcmcMachinery_; delete mcmcSample_; delete dEploidIO_; @@ -113,7 +117,7 @@ class TestMcmcMachinery: public CppUnit::TestCase { void testMainConstructor(){ - McmcMachinery tmpMcmcMachinery(this->dEploidIO_, this->panel_, this->mcmcSample_ ); + McmcMachinery tmpMcmcMachinery(this->dEploidIO_, this->panel_, this->mcmcSample_, this->rg_ ); } From 4374e93b4e067b3049bdcfe46fa04e1cee720d4d Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sat, 8 Oct 2016 18:13:07 +0100 Subject: [PATCH 11/29] fix memory leak --- src/dEploid.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/dEploid.cpp b/src/dEploid.cpp index 977fa92d..91a873cc 100644 --- a/src/dEploid.cpp +++ b/src/dEploid.cpp @@ -63,9 +63,8 @@ int main( int argc, char *argv[] ){ mcmcMachinery.runMcmcChain(true); dEploidIO.write(mcmcSample, panel); - //rg->clearFastFunc(); - //delete rg; + rg.clearFastFunc(); if ( panel ){ delete panel; } From c78bc87e5d7e2229325c8765ca265a3422ac1bd6 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 00:32:08 +0100 Subject: [PATCH 12/29] progress --- src/dEploidIO.hpp | 3 +++ src/mcmc.cpp | 6 +++--- src/mcmc.hpp | 5 +++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/dEploidIO.hpp b/src/dEploidIO.hpp index bba9afb5..ee32a778 100644 --- a/src/dEploidIO.hpp +++ b/src/dEploidIO.hpp @@ -48,6 +48,7 @@ class DEploidIO{ friend class TestMcmcMachinery; #endif friend class McmcMachinery; + friend class RMcmcSample; public: DEploidIO(); DEploidIO(const std::string &arg); @@ -67,6 +68,8 @@ class DEploidIO{ string panelFileName_; size_t nLoci() const { return this->nLoci_; } + size_t kStrain() const { return this->kStrain_;} + size_t nMcmcSample() const { return this->nMcmcSample_; } double averageCentimorganDistance() const { return this->averageCentimorganDistance_; } double Ne() const { return this->Ne_; } double constRecombProb() const { return this->constRecombProb_; } diff --git a/src/mcmc.cpp b/src/mcmc.cpp index 820c7a38..fdc6ea52 100644 --- a/src/mcmc.cpp +++ b/src/mcmc.cpp @@ -31,9 +31,9 @@ McmcSample::McmcSample(){}; McmcSample::~McmcSample(){}; -McmcMachinery::McmcMachinery(DEploidIO* pfDeconfIO, Panel *panel, McmcSample *mcmcSample, RandomGenerator* rg_ ){ // initialiseMCMCmachinery +McmcMachinery::McmcMachinery(DEploidIO* dEploidIO, Panel *panel, McmcSample *mcmcSample, RandomGenerator* rg_ ){ // initialiseMCMCmachinery - this->dEploidIO_ = pfDeconfIO; + this->dEploidIO_ = dEploidIO; this->panel_ = panel; this->mcmcSample_ = mcmcSample; this->seed_ = this->dEploidIO_->randomSeed(); @@ -47,7 +47,7 @@ McmcMachinery::McmcMachinery(DEploidIO* pfDeconfIO, Panel *panel, McmcSample *mc //this->propRg_ = new MersenneTwister(this->seed_); //this->initialHapRg_ = new MersenneTwister(this->seed_); - this->calcMaxIteration( dEploidIO_->nMcmcSample_ , dEploidIO_->mcmcMachineryRate_, pfDeconfIO->mcmcBurn_ ); + this->calcMaxIteration( dEploidIO_->nMcmcSample_ , dEploidIO_->mcmcMachineryRate_, dEploidIO_->mcmcBurn_ ); this->MN_LOG_TITRE = 0.0; this->SD_LOG_TITRE = 3.0; diff --git a/src/mcmc.hpp b/src/mcmc.hpp index d22e6741..a0eaa2c5 100644 --- a/src/mcmc.hpp +++ b/src/mcmc.hpp @@ -42,6 +42,7 @@ class McmcSample { #endif friend class McmcMachinery; friend class DEploidIO; + friend class RMcmcSample; public: McmcSample(); ~McmcSample(); @@ -52,11 +53,11 @@ class McmcSample { } - private: - vector < int > moves; vector < vector > proportion; vector < vector > hap; vector < double > sumLLKs; + private: + vector < int > moves; }; From ca320b91f477e62909a63cae144bed272c7d2ec9 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 14:27:41 +0100 Subject: [PATCH 13/29] update version tag --- version | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/version b/version index 13210269..35258725 100644 --- a/version +++ b/version @@ -1 +1 @@ -v0.3-release +v0.3-beta.1 From 3790f6d312a949667c637d7881929605fa165fcd Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 14:35:55 +0100 Subject: [PATCH 14/29] plotProportions --- utilities/dEploidTools.r | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 5cae7c1c..11491ce6 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -347,9 +347,9 @@ plot.wsaf.vs.index <- function ( coverage, expWSAF = c(), expWSAFChrom = c(), ex points(plotIndex, expWSAF[expWSAFChrom == chromI], col="blue") } } - } + fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclude ){ dEploidOutput = fun.dEploidPrefix ( dEploidPrefix ) tmpProp = read.table(dEploidOutput$propFileName, header=F) @@ -368,27 +368,28 @@ fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclu } -#' Deconvolute mixed haplotypes, and reporting the mixture proportions from each sample -#' This function provieds an interface for calling \emph{dEploid} from R. -#' The command line options are passed via the \code{args} argument and \code{file} +#' @title Plot proportions #' -#' @section blahblah: -#' Blah blah +#' @description Plot the MCMC samples of the proportion, indexed by the MCMC chain. #' -#' @param tmpProp blah blah blah +#' @param proportions Matrix of the MCMC proportion samples. The matrix size is number of the MCMC samples by the number of strains. #' -#' @param title blah blah blah +#' @param title Figure title. #' -#' @return A named list of something something ... +#' @return #' #' @export #' #' @examples -#' set.seed(1234) +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' panelFile = system.file("extdata", "labStrains.test.panel.txt", package = "DEploid") +#' PG0390 = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) +#' plotProportions( PG0390$Proportions, "PG0390-C proportions" ) #' -plot.prob <-function (tmpProp, title){ +plotProportions <-function (proportions, title = ""){ rainbowColorBin = 16 - barplot(t(tmpProp), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="SNP index", ylab="Posterior probabilities", main=title) + barplot(t(proportions), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="SNP index", ylab="Posterior probabilities", main=title) } @@ -399,7 +400,7 @@ plot.postProb.ofCase <- function ( inPrefix, outPrefix, case ){ ncol = ceiling(length(chromName)/2) par(mfrow = c(ncol,length(chromName)/ncol)) for ( chromI in chromName ){ - plot.prob ( obj[which( chromI == obj$CHROM),c(3:dim(obj)[2])], "") + plotProportions ( obj[which( chromI == obj$CHROM),c(3:dim(obj)[2])], "") } dev.off() } From 1bdf4660083127677260dbce9b2e90277fb83b35 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 15:08:15 +0100 Subject: [PATCH 15/29] documentate extractCoverageFromVcf --- utilities/dEploidTools.r | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 11491ce6..de0fc678 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -80,7 +80,7 @@ fun.dEploidPrefix <- function ( prefix ){ fun.extract.coverage <- function ( inputs ){ if ( inputs$vcfFileName != "" ){ - return (fun.extract.vcf (inputs$vcfFileName)) + return (extractCoverageFromVcf (inputs$vcfFileName)) } else { return (fun.extract.coverage.from.txt (inputs$refFileName, inputs$altFileName)) } @@ -108,8 +108,25 @@ fun.extract.exclude <- function (excludeFileName, excludeBool){ } } - -fun.extract.vcf <- function ( vcfName, ADFieldIndex = 2 ){ +#' @title Extract read counts from VCF +#' +#' @description Extract read counts from VCF file of a single sample. +#' +#' @note The VCF file should only contain one sample. If more samples present in the VCF, it only returns coverage for of the first sample. +#' +#' @param vcfName Path of the VCF file. +#' +#' @param ADFieldIndex Index of the AD field of the sample field. For example, if the format is "GT:AD:DP:GQ:PL", the AD index is 2 (by default). +#' +#' @return A data.frame contains four columns: chromosomes, positions, reference allele count, alternative allele count. +#' +#' @export +#' +#' @examples +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390 = extractCoverageFromVcf(vcfFile) +#' +extractCoverageFromVcf <- function ( vcfName, ADFieldIndex = 2 ){ # Assume that AD is the second field catCmd = "cat" if ( grepl("gzip", system(paste("file --mime-type", vcfName), T) ) == TRUE ){ From 31b12ad07149b5b2cd58894a700b47860168bc89 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 16:11:45 +0100 Subject: [PATCH 16/29] document extractCoverageFromTxt --- utilities/dEploidTools.r | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index de0fc678..79f2bf6c 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -82,13 +82,31 @@ fun.extract.coverage <- function ( inputs ){ if ( inputs$vcfFileName != "" ){ return (extractCoverageFromVcf (inputs$vcfFileName)) } else { - return (fun.extract.coverage.from.txt (inputs$refFileName, inputs$altFileName)) + return (extractCoverageFromTxt (inputs$refFileName, inputs$altFileName)) } } - -fun.extract.coverage.from.txt <- function ( refFileName, altFileName ){ +#' @title Extract read counts from plain text file +#' +#' @description Extract read counts from tab-delimited text files of a single sample. +#' +#' @note The allele count files must be tab-delimited. The allele count files contain three columns: chromosomes, positions and allele count. +#' +#' @param refFileName Path of the reference allele count file. +#' +#' @param altFileName Path of the alternative allele count file. +#' +#' @return A data.frame contains four columns: chromosomes, positions, reference allele count, alternative allele count. +#' +#' @export +#' +#' @examples +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390 = extractCoverageFromTxt(refFile, altFile) +#' +extractCoverageFromTxt <- function ( refFileName, altFileName ){ ref = read.table(refFileName, header = TRUE, comment.char = "") alt = read.table(altFileName, header = TRUE, comment.char = "") return ( data.frame( CHROM = ref[,1], From 4ea531c215368632c9ba9bb99321abd267e03d74 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 17:59:40 +0100 Subject: [PATCH 17/29] document plotAltVsRef --- utilities/dEploidTools.r | 75 ++++++++++++++++++++++++++++++---------- 1 file changed, 56 insertions(+), 19 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 79f2bf6c..55f6293d 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -200,15 +200,56 @@ fun.dic.by.theta <- function ( tmpllk, thetallk ){ } -plot.prop <-function (propMat, title = "Components"){ - rainbowColorBin = 5 - barplot(t(propMat), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="Iteration", ylab="Component Freq", main=title) + +#' @title Plot proportions +#' +#' @description Plot the MCMC samples of the proportion, indexed by the MCMC chain. +#' +#' @param proportions Matrix of the MCMC proportion samples. The matrix size is number of the MCMC samples by the number of strains. +#' +#' @param title Figure title. +#' +#' @return +#' +#' @export +#' +#' @examples +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' panelFile = system.file("extdata", "labStrains.test.panel.txt", package = "DEploid") +#' PG0390 = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) +#' plotProportions( PG0390$Proportions, "PG0390-C proportions" ) +#' +plotProportions <-function (proportions, title = "Components"){ + rainbowColorBin = 16 + barplot(t(proportions), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="Iteration", ylab="Component proportion", main=title) } -plot.altVsRef <- function ( ref, alt, title = "Alt vs Ref", exclude.ref = c(), exclude.alt = c() ){ +#' @title Plot coverage +#' +#' @description Plot alternative allele count vs reference allele count at each site. +#' +#' @param ref Numeric array of reference allele count. +#' +#' @param alt Numeric array of alternative allele count. +#' +#' @param title Figure title, "Alt vs Ref" by default +#' +#' @param exclude.ref @param exclude.alt +#' +#' @return +#' +#' @export +#' +#' @examples +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390 = extractCoverageFromVcf(vcfFile) +#' plotAltVsRef( PG0390$refCount, PG0390$altCount ) +#' +plotAltVsRef <- function ( ref, alt, title = "Alt vs Ref", exclude.ref = c(), exclude.alt = c() ){ tmp.range = 1.1*mean(max(alt), max(ref)) - plot ( ref, alt, xlim=c(0,tmp.range), ylim=c(0,tmp.range), cex = 0.5, xlab = "REF", ylab = "ALT", main = title) + plot ( ref, alt, xlim=c(0, tmp.range), ylim=c(0,tmp.range), cex = 0.5, xlab = "REF", ylab = "ALT", main = title) points (exclude.ref, exclude.alt, col = "red") abline(v =50, untf = FALSE, lty = 2) abline(h =50, untf = FALSE, lty = 2) @@ -303,7 +344,7 @@ fun.dataExplore <- function (coverage, plafInfo, prefix = "") { png ( paste ( prefix, "altVsRefAndWSAFvsPLAF.png", sep = "" ), width = 1800, height = 600) par( mfrow = c(1,3) ) - plot.altVsRef ( ref, alt ) + plotAltVsRef ( ref, alt ) WSAF = fun.calc.obsWSAF ( alt, ref ) @@ -331,7 +372,7 @@ fun.interpretDEploid.1 <- function (coverage, plafInfo, dEploidPrefix, prefix = # png ( paste ( prefix, ".interpretDEploidFigure.1.png", sep = "" ), width = 1500, height = 500, bg="transparent") png ( paste ( prefix, ".interpretDEploidFigure.1.png", sep = "" ), width = 1500, height = 1000) par( mfrow = c(2,3) ) - plot.altVsRef ( ref, alt ) + plotAltVsRef ( ref, alt ) obsWSAF = fun.calc.obsWSAF ( alt, ref ) plot.wsaf.hist ( obsWSAF ) @@ -347,7 +388,7 @@ fun.interpretDEploid.1 <- function (coverage, plafInfo, dEploidPrefix, prefix = } plot.plaf.vs.wsaf ( PLAF, obsWSAF, expWSAF ) - plot.prop( tmpProp ) + plotProportions( tmpProp ) tmpTitle = fun.getWSAF.corr (obsWSAF, expWSAF, dEploidOutput$dicLogFileName) plot.wsaf ( obsWSAF, expWSAF, tmpTitle ) @@ -403,11 +444,11 @@ fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclu } -#' @title Plot proportions +#' @title Painting haplotype according the reference panel #' -#' @description Plot the MCMC samples of the proportion, indexed by the MCMC chain. +#' @description Plot the posterior probabilities of a haplotype given the refernece panel. #' -#' @param proportions Matrix of the MCMC proportion samples. The matrix size is number of the MCMC samples by the number of strains. +#' @param posteriorProbabilities Posterior probabilities matrix with the size of number of loci by the number of reference strain. #' #' @param title Figure title. #' @@ -416,15 +457,11 @@ fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclu #' @export #' #' @examples -#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") -#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' panelFile = system.file("extdata", "labStrains.test.panel.txt", package = "DEploid") -#' PG0390 = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) -#' plotProportions( PG0390$Proportions, "PG0390-C proportions" ) +#' to do ... #' -plotProportions <-function (proportions, title = ""){ +haplotypePainter <-function (posteriorProbabilities, title = ""){ rainbowColorBin = 16 - barplot(t(proportions), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="SNP index", ylab="Posterior probabilities", main=title) + barplot(t(posteriorProbabilities), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="SNP index", ylab="Posterior probabilities", main=title) } @@ -435,7 +472,7 @@ plot.postProb.ofCase <- function ( inPrefix, outPrefix, case ){ ncol = ceiling(length(chromName)/2) par(mfrow = c(ncol,length(chromName)/ncol)) for ( chromI in chromName ){ - plotProportions ( obj[which( chromI == obj$CHROM),c(3:dim(obj)[2])], "") + haplotypePainter ( obj[which( chromI == obj$CHROM),c(3:dim(obj)[2])], "") } dev.off() } From 24ebe59f5cd6b4c3ca021de02c7833a9202e6376 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 18:03:45 +0100 Subject: [PATCH 18/29] document plotAltVsRef --- utilities/dEploidTools.r | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 55f6293d..5d8273df 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -236,7 +236,9 @@ plotProportions <-function (proportions, title = "Components"){ #' #' @param title Figure title, "Alt vs Ref" by default #' -#' @param exclude.ref @param exclude.alt +#' @param exclude.ref Numeric array of reference allele count at sites that are not deconvoluted. +#' +#' @param exclude.alt Numeric array of alternative allele count at sites that are not deconvoluted #' #' @return #' From 24e30de527418e4f092d66b8c513ff24647c1b61 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 18:49:23 +0100 Subject: [PATCH 19/29] document computeObsWSAF --- utilities/dEploidTools.r | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 5d8273df..8c535d40 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -333,7 +333,26 @@ plot.wsaf <- function (obsWSAF, expWSAF, title = ""){ } -fun.calc.obsWSAF <- function (alt, ref) { +#' @title Compute observed WSAF +#' +#' @description Compute observed allele frequency within sample from the allele counts. +#' +#' @param ref Numeric array of reference allele count. +#' +#' @param alt Numeric array of alternative allele count. +#' +#' @return Numeric array of observed allele frequency within sample. +#' +#' @seealso \code{\link{histWSAF}} for histogram, \code{\link{plotWSAF}} for WSAFs plotted against indices. +#' +#' @export +#' +#' @examples +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390 = extractCoverageFromVcf(vcfFile) +#' computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' +computeObsWSAF <- function (alt, ref) { return ( alt / (ref + alt + 0.00000001) ) } @@ -348,11 +367,11 @@ fun.dataExplore <- function (coverage, plafInfo, prefix = "") { plotAltVsRef ( ref, alt ) - WSAF = fun.calc.obsWSAF ( alt, ref ) + obsWSAF = computeObsWSAF ( alt, ref ) - plot.wsaf.hist ( WSAF ) + plot.wsaf.hist ( obsWSAF ) - plot.plaf.vs.wsaf ( PLAF, WSAF ) + plot.plaf.vs.wsaf ( PLAF, obsWSAF ) dev.off() } @@ -376,7 +395,7 @@ fun.interpretDEploid.1 <- function (coverage, plafInfo, dEploidPrefix, prefix = par( mfrow = c(2,3) ) plotAltVsRef ( ref, alt ) - obsWSAF = fun.calc.obsWSAF ( alt, ref ) + obsWSAF = computeObsWSAF ( alt, ref ) plot.wsaf.hist ( obsWSAF ) if (exclude$excludeBool){ @@ -406,7 +425,7 @@ plot.wsaf.vs.index <- function ( coverage, expWSAF = c(), expWSAFChrom = c(), ex chromList = levels(coverage$CHROM) ref = coverage$refCount alt = coverage$altCount - obsWSAF = fun.calc.obsWSAF ( alt, ref ) + obsWSAF = computeObsWSAF ( alt, ref ) for ( chromI in chromList ){ plot( obsWSAF[coverage$CHROM==chromI], col="red", ylim=c(0,1), main = paste(titlePrefix, chromI, "WSAF")) From 72f37c800333148915cbff78664d2a774432f013 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 19:29:15 +0100 Subject: [PATCH 20/29] document histWSAF --- utilities/dEploidTools.r | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 8c535d40..c0f3146d 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -261,12 +261,33 @@ plotAltVsRef <- function ( ref, alt, title = "Alt vs Ref", exclude.ref = c(), ex } -plot.wsaf.hist <- function ( obsWSAF, indexed = TRUE, title ="Histogram 00) ) == 1) } - hist(obsWSAF, main=title, breaks = seq(0, 1, by =0.1), xlab = "WSAF") + return (hist(obsWSAF[tmpWSAF_index], main=title, breaks = seq(0, 1, by =0.1), xlab = "WSAF")) } @@ -369,7 +390,7 @@ fun.dataExplore <- function (coverage, plafInfo, prefix = "") { obsWSAF = computeObsWSAF ( alt, ref ) - plot.wsaf.hist ( obsWSAF ) + histWSAF ( obsWSAF ) plot.plaf.vs.wsaf ( PLAF, obsWSAF ) @@ -396,7 +417,7 @@ fun.interpretDEploid.1 <- function (coverage, plafInfo, dEploidPrefix, prefix = plotAltVsRef ( ref, alt ) obsWSAF = computeObsWSAF ( alt, ref ) - plot.wsaf.hist ( obsWSAF ) + histWSAF ( obsWSAF ) if (exclude$excludeBool){ excludeLogic = ( paste(coverage$CHROM, coverage$POS) %in% paste(exclude$excludeTable$CHROM, exclude$excludeTable$POS) ) From db3819dfab3af6b776854bb2e1aadb41b2e1d724 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 20:17:09 +0100 Subject: [PATCH 21/29] extract plaf and plot plaf vs wsaf --- utilities/dEploidTools.r | 58 ++++++++++++++++++++++++++++++------ utilities/dataExplore.r | 2 +- utilities/interpretDEploid.r | 2 +- 3 files changed, 51 insertions(+), 11 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index c0f3146d..8c966ef7 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -173,8 +173,24 @@ extractCoverageFromVcf <- function ( vcfName, ADFieldIndex = 2 ){ } -fun.extract.plaf <- function ( plafName ){ - return ( read.table(plafName, header=T) ) +#' @title Extract PLAF +#' +#' @description Extract population level allele frequency (PLAF) from text file. +#' +#' @note The text file must have header, and population level allele frequency recorded in the "PLAF" field. +#' +#' @param plafName Path of the PLAF text file. +#' +#' @return A numeric array of PLAF +#' +#' @export +#' +#' @examples +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") +#' plaf = extractPLAF(plafFile) +#' +extractPLAF<- function ( plafName ){ + return ( read.table(plafName, header=T)$PLAF ) } @@ -328,7 +344,31 @@ fun.getllk.dic <- function ( llkTable, ref, alt, expWSAF, logFileName ){ } -plot.plaf.vs.wsaf <- function ( plaf, obsWSAF, expWSAF = c(), title = "PLAF vs WSAF" ){ +#' @title Plot WSAF vs PLAF +#' +#' @description Plot allele frequencies within sample against population level. +#' +#' @param plaf Numeric array of population level allele frequency. +#' +#' @param obsWSAF Numeric array of observed altenative allele frequencies within sample. +#' +#' @param expWSAF Numeric array of expected WSAF from model. +#' +#' @param title Figure title, "WSAF vs PLAF" by default +#' +#' @return +#' +#' @export +#' +#' @examples +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390 = extractCoverageFromVcf(vcfFile) +#' obsWSAF = computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") +#' plaf = extractPLAF(plafFile) +#' plotWSAFvsPLAF(plaf, obsWSAF) +#' +plotWSAFvsPLAF <- function ( plaf, obsWSAF, expWSAF = c(), title = "WSAF vs PLAF" ){ plot ( plaf, obsWSAF, cex = 0.5, xlim = c(0,1), ylim = c(0,1), col = "red", main = title, xlab = "PLAF", ylab = "WSAF" ) if ( length(expWSAF) > 0 ){ points ( plaf, expWSAF, cex = 0.5, col = "blue") @@ -378,8 +418,8 @@ computeObsWSAF <- function (alt, ref) { } -fun.dataExplore <- function (coverage, plafInfo, prefix = "") { - PLAF = plafInfo$PLAF +fun.dataExplore <- function (coverage, PLAF, prefix = "") { +# PLAF = plafInfo$PLAF ref = coverage$refCount alt = coverage$altCount @@ -392,15 +432,15 @@ fun.dataExplore <- function (coverage, plafInfo, prefix = "") { histWSAF ( obsWSAF ) - plot.plaf.vs.wsaf ( PLAF, obsWSAF ) + plotWSAFvsPLAF ( PLAF, obsWSAF ) dev.off() } -fun.interpretDEploid.1 <- function (coverage, plafInfo, dEploidPrefix, prefix = "", exclude ) { +fun.interpretDEploid.1 <- function (coverage, PLAF, dEploidPrefix, prefix = "", exclude ) { - PLAF = plafInfo$PLAF +# PLAF = plafInfo$PLAF ref = coverage$refCount alt = coverage$altCount @@ -428,7 +468,7 @@ fun.interpretDEploid.1 <- function (coverage, plafInfo, dEploidPrefix, prefix = ref = ref[includeindex] alt = alt[includeindex] } - plot.plaf.vs.wsaf ( PLAF, obsWSAF, expWSAF ) + plotWSAFvsPLAF ( PLAF, obsWSAF, expWSAF ) plotProportions( tmpProp ) diff --git a/utilities/dataExplore.r b/utilities/dataExplore.r index 6a09bc05..3c2b8a4b 100644 --- a/utilities/dataExplore.r +++ b/utilities/dataExplore.r @@ -20,7 +20,7 @@ myInput = fun.parse ( args ) myCoverageInfo = fun.extract.coverage ( myInput ) -myPlafInfo = fun.extract.plaf ( myInput$plafFileName ) +myPlafInfo = extractPLAF( myInput$plafFileName ) fun.dataExplore (myCoverageInfo, myPlafInfo, myInput$outPrefix) diff --git a/utilities/interpretDEploid.r b/utilities/interpretDEploid.r index f7587935..b62149f7 100644 --- a/utilities/interpretDEploid.r +++ b/utilities/interpretDEploid.r @@ -20,7 +20,7 @@ myInput = fun.parse ( args ) myCoverageInfo = fun.extract.coverage ( myInput ) -myPlafInfo = fun.extract.plaf ( myInput$plafFileName ) +myPlafInfo = extractPLAF( myInput$plafFileName ) myExcludeInfo = fun.extract.exclude (myInput$excludeFileName, myInput$excludeBool) From f60c63d5834609323e939c7f7285501a36b0d97d Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 21:44:01 +0100 Subject: [PATCH 22/29] plot WSAF obs vs exp --- utilities/dEploidTools.r | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 8c966ef7..3b2d2b5b 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -385,7 +385,30 @@ fun.getWSAF.corr <- function( obsWSAF, expWSAF, dicLogFileName ){ } -plot.wsaf <- function (obsWSAF, expWSAF, title = ""){ +#' @title Plot WSAF +#' +#' @description Plot observed alternative allele frequency within sample against expected WSAF. +#' +#' @param obsWSAF Numeric array of observed WSAF. +#' +#' @param expWSAF Numeric array of expected WSAF. +#' +#' @param title Figure title. +#' +#' @return +#' +#' @export +#' +#' @examples +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390 = extractCoverageFromVcf(vcfFile) +#' obsWSAF = computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' PG0390.deconv = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) +#' prop = PG0390.deconv$Proportions[dim(PG0390.deconv$Proportions)[1],] +#' expWSAF = t(PG0390.deconv$Haps) %*% prop +#' plotObsExpWSAF(obsWSAF, expWSAF) +#' +plotObsExpWSAF <- function (obsWSAF, expWSAF, title = "WSAF(observed vs expected)"){ plot(obsWSAF, expWSAF, pch=19, col="blue", xlab="Observed WSAF (ALT/(ALT+REF))", ylab="Expected WSAF (h%*%p)", main=title, xlim = c(-0.05, 1.05), cex = 0.5, ylim = c(-0.05, 1.05)); @@ -473,7 +496,7 @@ fun.interpretDEploid.1 <- function (coverage, PLAF, dEploidPrefix, prefix = "", plotProportions( tmpProp ) tmpTitle = fun.getWSAF.corr (obsWSAF, expWSAF, dEploidOutput$dicLogFileName) - plot.wsaf ( obsWSAF, expWSAF, tmpTitle ) + plotObsExpWSAF ( obsWSAF, expWSAF, tmpTitle ) tmpTitle = fun.getllk.dic (llkTable, ref, alt, expWSAF, dEploidOutput$dicLogFileName ) plot.llk( llkTable, ref, alt, expWSAF, tmpTitle ) From f8919e33b493c7df956ac054df0389545d759ee8 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 22:17:28 +0100 Subject: [PATCH 23/29] remove return --- utilities/dEploidTools.r | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 3b2d2b5b..cacc0c1a 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -225,8 +225,6 @@ fun.dic.by.theta <- function ( tmpllk, thetallk ){ #' #' @param title Figure title. #' -#' @return -#' #' @export #' #' @examples @@ -256,8 +254,6 @@ plotProportions <-function (proportions, title = "Components"){ #' #' @param exclude.alt Numeric array of alternative allele count at sites that are not deconvoluted #' -#' @return -#' #' @export #' #' @examples @@ -356,8 +352,6 @@ fun.getllk.dic <- function ( llkTable, ref, alt, expWSAF, logFileName ){ #' #' @param title Figure title, "WSAF vs PLAF" by default #' -#' @return -#' #' @export #' #' @examples @@ -395,8 +389,6 @@ fun.getWSAF.corr <- function( obsWSAF, expWSAF, dicLogFileName ){ #' #' @param title Figure title. #' -#' @return -#' #' @export #' #' @examples @@ -557,12 +549,9 @@ fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclu #' #' @param title Figure title. #' -#' @return -#' #' @export #' #' @examples -#' to do ... #' haplotypePainter <-function (posteriorProbabilities, title = ""){ rainbowColorBin = 16 From 604c4d0d78f91d31dd21b2c8915a9971f13fba2e Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 22:31:56 +0100 Subject: [PATCH 24/29] update --- utilities/dEploidTools.r | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index cacc0c1a..16afdd1c 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -419,7 +419,7 @@ plotObsExpWSAF <- function (obsWSAF, expWSAF, title = "WSAF(observed vs expected #' #' @return Numeric array of observed allele frequency within sample. #' -#' @seealso \code{\link{histWSAF}} for histogram, \code{\link{plotWSAF}} for WSAFs plotted against indices. +#' @seealso \code{\link{histWSAF}} for histogram. #' #' @export #' From cb9c9a13e33941013fa6bf6a4282e120f37a15cf Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Sun, 9 Oct 2016 22:38:23 +0100 Subject: [PATCH 25/29] update --- utilities/dEploidTools.r | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 16afdd1c..7c0f6398 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -395,6 +395,7 @@ fun.getWSAF.corr <- function( obsWSAF, expWSAF, dicLogFileName ){ #' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") #' PG0390 = extractCoverageFromVcf(vcfFile) #' obsWSAF = computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") #' PG0390.deconv = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) #' prop = PG0390.deconv$Proportions[dim(PG0390.deconv$Proportions)[1],] #' expWSAF = t(PG0390.deconv$Haps) %*% prop @@ -551,8 +552,6 @@ fun.interpretDEploid.2 <- function ( coverage, dEploidPrefix, prefix = "", exclu #' #' @export #' -#' @examples -#' haplotypePainter <-function (posteriorProbabilities, title = ""){ rainbowColorBin = 16 barplot(t(posteriorProbabilities), beside=F, border=NA, col=rainbow(rainbowColorBin), space=0, xlab="SNP index", ylab="Posterior probabilities", main=title) From 4b15ddeae18b85dc4bceabc96bb2c5a68e81b052 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Mon, 10 Oct 2016 13:15:03 +0100 Subject: [PATCH 26/29] change zcat to zcat <, so it will work on mac --- utilities/dEploidTools.r | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 7c0f6398..e382b820 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -148,7 +148,7 @@ extractCoverageFromVcf <- function ( vcfName, ADFieldIndex = 2 ){ # Assume that AD is the second field catCmd = "cat" if ( grepl("gzip", system(paste("file --mime-type", vcfName), T) ) == TRUE ){ - catCmd = "zcat" + catCmd = "zcat <" } skipNum = as.numeric(system(paste(catCmd, vcfName, " | head -5000 | grep \"##\" | wc -l"), T)) From c6322ae312abb6955546b8c576d70c04358ee4ce Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Thu, 13 Oct 2016 19:21:19 +0100 Subject: [PATCH 27/29] update examples --- utilities/dEploidTools.r | 87 ++++++++++++++++++++++++++++++++-------- 1 file changed, 71 insertions(+), 16 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index e382b820..8a040941 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -229,10 +229,20 @@ fun.dic.by.theta <- function ( tmpllk, thetallk ){ #' #' @examples #' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") -#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") #' panelFile = system.file("extdata", "labStrains.test.panel.txt", package = "DEploid") -#' PG0390 = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) -#' plotProportions( PG0390$Proportions, "PG0390-C proportions" ) +#' +#' # Example 1 +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390CoverageTxt = extractCoverageFromTxt(refFile, altFile) +#' PG0390CoverageTxt.deconv = dEploid(paste("-ref", refFile, "-alt", altFile, "-plaf", plafFile, "-noPanel")) +#' plotProportions( PG0390CoverageTxt.deconv$Proportions, "PG0390-C proportions" ) +#' +#' # Example 2 +#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390CoverageVcf = extractCoverageFromVcf(vcfFile) +#' PG0390CoverageVcf.deconv = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) +#' plotProportions( PG0390CoverageVcf.deconv$Proportions, "PG0390-C proportions" ) #' plotProportions <-function (proportions, title = "Components"){ rainbowColorBin = 16 @@ -257,9 +267,16 @@ plotProportions <-function (proportions, title = "Components"){ #' @export #' #' @examples +#' # Example 1 +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390CoverageTxt = extractCoverageFromTxt(refFile, altFile) +#' plotAltVsRef( PG0390CoverageTxt$refCount, PG0390CoverageTxt$altCount ) +#' +#' # Example 2 #' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' PG0390 = extractCoverageFromVcf(vcfFile) -#' plotAltVsRef( PG0390$refCount, PG0390$altCount ) +#' PG0390CoverageVcf = extractCoverageFromVcf(vcfFile) +#' plotAltVsRef( PG0390CoverageVcf$refCount, PG0390CoverageVcf$altCount ) #' plotAltVsRef <- function ( ref, alt, title = "Alt vs Ref", exclude.ref = c(), exclude.alt = c() ){ tmp.range = 1.1*mean(max(alt), max(ref)) @@ -288,9 +305,18 @@ plotAltVsRef <- function ( ref, alt, title = "Alt vs Ref", exclude.ref = c(), ex #' @export #' #' @examples +#' # Example 1 +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390CoverageTxt = extractCoverageFromTxt(refFile, altFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageTxt$altCount, PG0390CoverageTxt$refCount ) +#' histWSAF(obsWSAF) +#' myhist = histWSAF(obsWSAF, FALSE) +#' +#' # Example 2 #' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' PG0390 = extractCoverageFromVcf(vcfFile) -#' obsWSAF = computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' PG0390CoverageVcf = extractCoverageFromVcf(vcfFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount ) #' histWSAF(obsWSAF) #' myhist = histWSAF(obsWSAF, FALSE) #' @@ -355,9 +381,19 @@ fun.getllk.dic <- function ( llkTable, ref, alt, expWSAF, logFileName ){ #' @export #' #' @examples +#' # Example 1 +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390CoverageTxt = extractCoverageFromTxt(refFile, altFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageTxt$altCount, PG0390CoverageTxt$refCount ) +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") +#' plaf = extractPLAF(plafFile) +#' plotWSAFvsPLAF(plaf, obsWSAF) +#' +#' # Example 2 #' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' PG0390 = extractCoverageFromVcf(vcfFile) -#' obsWSAF = computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' PG0390CoverageVcf = extractCoverageFromVcf(vcfFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount ) #' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") #' plaf = extractPLAF(plafFile) #' plotWSAFvsPLAF(plaf, obsWSAF) @@ -392,13 +428,25 @@ fun.getWSAF.corr <- function( obsWSAF, expWSAF, dicLogFileName ){ #' @export #' #' @examples +#' # Example 1 +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390CoverageTxt = extractCoverageFromTxt(refFile, altFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageTxt$altCount, PG0390CoverageTxt$refCount ) +#' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") +#' PG0390CoverageTxt.deconv = dEploid(paste("-ref", refFile, "-alt", altFile, "-plaf", plafFile, "-noPanel")) +#' prop = PG0390CoverageTxt.deconv$Proportions[dim(PG0390CoverageTxt.deconv$Proportions)[1],] +#' expWSAF = t(PG0390CoverageTxt.deconv$Haps) %*% prop +#' plotObsExpWSAF(obsWSAF, expWSAF) +#' +#' # Example 2 #' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' PG0390 = extractCoverageFromVcf(vcfFile) -#' obsWSAF = computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' PG0390CoverageVcf = extractCoverageFromVcf(vcfFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount ) #' plafFile = system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid") -#' PG0390.deconv = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) -#' prop = PG0390.deconv$Proportions[dim(PG0390.deconv$Proportions)[1],] -#' expWSAF = t(PG0390.deconv$Haps) %*% prop +#' PG0390CoverageVcf.deconv = dEploid(paste("-vcf", vcfFile, "-plaf", plafFile, "-noPanel")) +#' prop = PG0390CoverageVcf.deconv$Proportions[dim(PG0390CoverageVcf.deconv$Proportions)[1],] +#' expWSAF = t(PG0390CoverageVcf.deconv$Haps) %*% prop #' plotObsExpWSAF(obsWSAF, expWSAF) #' plotObsExpWSAF <- function (obsWSAF, expWSAF, title = "WSAF(observed vs expected)"){ @@ -425,9 +473,16 @@ plotObsExpWSAF <- function (obsWSAF, expWSAF, title = "WSAF(observed vs expected #' @export #' #' @examples +#' # Example 1 +#' refFile = system.file("extdata", "PG0390-C.test.ref", package = "DEploid") +#' altFile = system.file("extdata", "PG0390-C.test.alt", package = "DEploid") +#' PG0390CoverageTxt = extractCoverageFromTxt(refFile, altFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageTxt$altCount, PG0390CoverageTxt$refCount ) +#' +#' # Example 2 #' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' PG0390 = extractCoverageFromVcf(vcfFile) -#' computeObsWSAF( PG0390$altCount, PG0390$refCount ) +#' PG0390CoverageVcf = extractCoverageFromVcf(vcfFile) +#' obsWSAF = computeObsWSAF( PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount ) #' computeObsWSAF <- function (alt, ref) { return ( alt / (ref + alt + 0.00000001) ) From c2f3236c53713def2681db07c43b3db3d248d914 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Thu, 13 Oct 2016 20:34:55 +0100 Subject: [PATCH 28/29] do not use hard coded gz --- utilities/dEploidTools.r | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/utilities/dEploidTools.r b/utilities/dEploidTools.r index 8a040941..e37c7e74 100644 --- a/utilities/dEploidTools.r +++ b/utilities/dEploidTools.r @@ -146,12 +146,24 @@ fun.extract.exclude <- function (excludeFileName, excludeBool){ #' extractCoverageFromVcf <- function ( vcfName, ADFieldIndex = 2 ){ # Assume that AD is the second field - catCmd = "cat" - if ( grepl("gzip", system(paste("file --mime-type", vcfName), T) ) == TRUE ){ - catCmd = "zcat <" + h <- function(w){ + if( any( grepl( "gzfile connection", w) ) ) + invokeRestart( "muffleWarning" ) } - skipNum = as.numeric(system(paste(catCmd, vcfName, " | head -5000 | grep \"##\" | wc -l"), T)) + gzf = gzfile(vcfFile, open = "rb") + skipNum = 0 + line = withCallingHandlers( readLines(gzf, n=1), warning=h) + while ( length(line) > 0 ){ + if (grepl("##", line )){ + skipNum = skipNum+1 + } else { + break + } + line = withCallingHandlers( readLines(gzf, n=1), warning=h) + } + close(gzf) + vcf = read.table( gzfile(vcfName), skip=skipNum, header=T, comment.char="", stringsAsFactors = FALSE, check.names=FALSE) sampleName = names(vcf)[10] From 57feebfddd07b3dbd66773e2a099e87f1a39a6a3 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Fri, 14 Oct 2016 12:14:02 +0100 Subject: [PATCH 29/29] only use valgrind on linux --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 40a22f93..5817d067 100644 --- a/.travis.yml +++ b/.travis.yml @@ -18,7 +18,7 @@ before_script: script: - make - make check - - ./tests/test_binary.sh + - if [ $TRAVIS_OS_NAME == linux ]; then ./tests/test_binary.sh; fi - ./tests/testPOS.sh - ./tests/test_binaryVcfVsTxt.sh - ./tests/test-against-previous-version.sh