From 480f83193a5524a474a8c52a97880e5639a19245 Mon Sep 17 00:00:00 2001 From: Jeremiah Wala Date: Thu, 22 Sep 2016 23:41:32 -0400 Subject: [PATCH] fixed unit tests --- .travis.scripts/publish-doxygen.sh | 4 +- SeqLib/GenomicRegionCollection.cpp | 188 ++++++++++++++--------------- bwa | 2 +- seq_test/seq_test.cpp | 64 +++------- seq_test/test_data/test.bed | 3 - seq_test/test_data/test.vcf | 151 ----------------------- src/BamWriter.cpp | 10 +- 7 files changed, 118 insertions(+), 304 deletions(-) delete mode 100644 seq_test/test_data/test.bed delete mode 100644 seq_test/test_data/test.vcf diff --git a/.travis.scripts/publish-doxygen.sh b/.travis.scripts/publish-doxygen.sh index 61e2a613e..716748206 100755 --- a/.travis.scripts/publish-doxygen.sh +++ b/.travis.scripts/publish-doxygen.sh @@ -7,13 +7,13 @@ cd ${HOME}; wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/travis/doxygen_1.8.8-1_amd64.deb sudo dpkg --install doxygen_1.8.8-1_amd64.deb - cd ${HOME}/build/jwalabroad/SeqLib; + cd ${HOME}/build/walaj/SeqLib; doxygen echo -e "Publishing doxygen...\n"; git config --global user.email "travis@travis-ci.org"; git config --global user.name "travis-ci"; - git clone --branch=gh-pages https://${GH_TOKEN}@github.com/jwalabroad/SeqLib gh-pages; + git clone --branch=gh-pages https://${GH_TOKEN}@github.com/walaj/SeqLib gh-pages; cd gh-pages; rm -rf doxygen/; mv ../docs/html doxygen/; diff --git a/SeqLib/GenomicRegionCollection.cpp b/SeqLib/GenomicRegionCollection.cpp index a41f42c95..2677ee5e7 100644 --- a/SeqLib/GenomicRegionCollection.cpp +++ b/SeqLib/GenomicRegionCollection.cpp @@ -64,8 +64,6 @@ namespace SeqLib { } } - static bool header_has_chr_string = false; - template void GenomicRegionCollection::CoordinateSort() { @@ -169,54 +167,65 @@ template bool GenomicRegionCollection::ReadBED(const std::string & file, const BamHeader& hdr) { m_sorted = false; - std::ifstream iss(file.c_str()); - if (!iss || file.length() == 0) { - std::cerr << "BED file does not exist: " << file << std::endl; + idx = 0; + + gzFile fp = NULL; + fp = strcmp(file.c_str(), "-")? gzopen(file.c_str(), "r") : gzdopen(fileno(stdin), "r"); + + if (file.empty() || !fp) { + std::cerr << "BED file not readable: " << file << std::endl; return false; } - std::string line; - std::string curr_chr = "-1"; - while (std::getline(iss, line, '\n')) { + // http://www.lemoda.net/c/gzfile-read/ + size_t c = 0; + while (1) { + + int err; + char buffer[GZBUFFER]; + gzgets(fp, buffer, GZBUFFER); + int bytes_read = strlen(buffer); + + // get one line + if (bytes_read < GZBUFFER - 1) { + if (gzeof (fp)) break; + else { + const char * error_string; + error_string = gzerror (fp, &err); + if (err) { + fprintf (stderr, "Error: %s.\n", error_string); + exit (EXIT_FAILURE); + } + } + } + // prepare to loop through each field of BED line size_t counter = 0; std::string chr, pos1, pos2; + std::string line(buffer); std::istringstream iss_line(line); std::string val; + if (line.find("#") != std::string::npos) + continue; - if (line.find("#") == std::string::npos) { - while(std::getline(iss_line, val, '\t')) { - switch (counter) { - case 0 : chr = val; break; - case 1 : pos1 = val; break; - case 2 : pos2 = val; break; - } - if (counter >= 2) - break; - ++counter; - - if (chr != curr_chr) { - //std::cerr << "...reading from BED - chr" << chr << std::endl; - curr_chr = chr; - } - - } - - if (header_has_chr_string) { - //if (chr == "X" || chr == "Y" || std::stoi(chr) < 23) - chr = "chr" + chr; + // loop BED lines + while(std::getline(iss_line, val, '\t')) { + switch (counter) { + case 0 : chr = val; break; + case 1 : pos1 = val; break; + case 2 : pos2 = val; break; } - - // construct the GenomicRegion - T gr(chr, pos1, pos2, hdr); - - if (gr.chr >= 0) { - m_grv->push_back(gr); - } - - //} - } // end "keep" conditional - } // end main while + if (counter >= 2) + break; + ++counter; + } + + // construct the GenomicRegion + T gr(chr, pos1, pos2, hdr); + + if (gr.chr >= 0) + m_grv->push_back(gr); + } return true; } @@ -225,62 +234,14 @@ template bool GenomicRegionCollection::ReadVCF(const std::string & file, const BamHeader& hdr) { m_sorted = false; - std::ifstream iss(file.c_str()); - if (!iss || file.length() == 0) { - std::cerr << "VCF file does not exist: " << file << std::endl; - return false; - } - - std::string line; - - while (std::getline(iss, line, '\n')) { - if (line.length() > 0) { - if (line.at(0) != '#') { // its a valid line - std::istringstream iss_this(line); - int count = 0; - std::string val, chr, pos; - - while (std::getline(iss_this, val, '\t')) { - switch (count) { - case 0 : chr = val; - case 1 : pos = val; - } - count++; - } - if (count < 3) { - std::cerr << "Didn't parse VCF line properly: " << line << std::endl; - return false; - } else { - - // construct the GenomicRegion - T gr(chr, pos, pos, hdr); - if (gr.chr >= 0) { - m_grv->push_back(gr); - } - } - - } - } - } - - return true; - -} - -template -GenomicRegionCollection::GenomicRegionCollection(const std::string &file, const BamHeader& hdr) { - - __allocate_grc(); - idx = 0; - // gzFile fp = NULL; fp = strcmp(file.c_str(), "-")? gzopen(file.c_str(), "r") : gzdopen(fileno(stdin), "r"); if (file.empty() || !fp) { - std::cerr << "Region file not readable: " << file << std::endl; - return; + std::cerr << "VCF file not readable: " << file << std::endl; + return false; } // http://www.lemoda.net/c/gzfile-read/ @@ -292,6 +253,7 @@ GenomicRegionCollection::GenomicRegionCollection(const std::string &file, con gzgets(fp, buffer, GZBUFFER); int bytes_read = strlen(buffer); + // get one line if (bytes_read < GZBUFFER - 1) { if (gzeof (fp)) break; else { @@ -304,9 +266,42 @@ GenomicRegionCollection::GenomicRegionCollection(const std::string &file, con } } - std::cerr << std::string(buffer) << " C " << c << std::endl; + // prepare to loop through each field of BED line + size_t counter = 0; + std::string chr, pos; + std::string line(buffer); + std::istringstream iss_line(line); + std::string val; + if (line.empty() || line.at(0) == '#') + continue; + + // loop VCF lines + while(std::getline(iss_line, val, '\t')) { + switch (counter) { + case 0 : chr = val; break; + case 1 : pos = val; break; + } + if (counter >= 2) + break; + ++counter; + } + + // construct the GenomicRegion + T gr(chr, pos, pos, hdr); + if (gr.chr >= 0) + m_grv->push_back(gr); } + return true; +} + +template +GenomicRegionCollection::GenomicRegionCollection(const std::string &file, const BamHeader& hdr) { + + __allocate_grc(); + + idx = 0; + /* // get the header line to check format std::string header; @@ -316,20 +311,17 @@ GenomicRegionCollection::GenomicRegionCollection(const std::string &file, con } iss.close(); - // MUTECT CALL STATS - if ((header.find("MuTect") != std::string::npos) || - (file.find("call_stats") != std::string::npos) || - (file.find("callstats") != std::string::npos)) - std::cerr << "MuTect reading not currently available" << std::endl; //ReadMuTect(file, hdr); + */ + // BED file - else if (file.find(".bed") != std::string::npos) + if (file.find(".bed") != std::string::npos) ReadBED(file, hdr); // VCF file else if (file.find(".vcf") != std::string::npos) ReadVCF(file, hdr); else // default is BED file ReadBED(file, hdr); - */ + } // reduce a set of GenomicRegions into the minium overlapping set (same as GenomicRanges "reduce") diff --git a/bwa b/bwa index 88a5ca582..fbd4dbc03 160000 --- a/bwa +++ b/bwa @@ -1 +1 @@ -Subproject commit 88a5ca582908ffc5dde06e5843068a6db561cdc5 +Subproject commit fbd4dbc03904eccd71cdca8cac7aa48da749c19c diff --git a/seq_test/seq_test.cpp b/seq_test/seq_test.cpp index e9ebbc2c9..4a550ce92 100644 --- a/seq_test/seq_test.cpp +++ b/seq_test/seq_test.cpp @@ -14,6 +14,7 @@ #include "SeqLib/RefGenome.h" #define GZBED "test_data/test.bed.gz" +#define GZVCF "test_data/test.vcf.gz" #define SBAM "test_data/small.bam" #define OBAM "test_data/small_out.bam" #define OCRAM "test_data/small_out.cram" @@ -36,7 +37,14 @@ BOOST_AUTO_TEST_CASE( read_gzbed ) { br.Open("test_data/small.bam"); SeqLib::GRC g(GZBED, br.Header()); + BOOST_CHECK_EQUAL(g.size(), 3); + + BOOST_CHECK_EQUAL(g[2].chr, 22); + + SeqLib::GRC v(GZVCF, br.Header()); + BOOST_CHECK_EQUAL(v.size(), 31); + BOOST_CHECK_EQUAL(v[29].chr, 22); } BOOST_AUTO_TEST_CASE ( bfc ) { @@ -405,15 +413,15 @@ BOOST_AUTO_TEST_CASE( bam_record ) { size_t count = 0; br.GetNextRecord(r); - BOOST_CHECK_EQUAL(r.asGenomicRegion().chr, 22); - BOOST_CHECK_EQUAL(r.asGenomicRegion().pos1,999901); - BOOST_CHECK_EQUAL(r.asGenomicRegion().pos2,1000002); - BOOST_CHECK_EQUAL(r.asGenomicRegion().strand,'+'); + BOOST_CHECK_EQUAL(r.AsGenomicRegion().chr, 22); + BOOST_CHECK_EQUAL(r.AsGenomicRegion().pos1,999901); + BOOST_CHECK_EQUAL(r.AsGenomicRegion().pos2,1000002); + BOOST_CHECK_EQUAL(r.AsGenomicRegion().strand,'+'); - BOOST_CHECK_EQUAL(r.asGenomicRegionMate().chr, 22); - BOOST_CHECK_EQUAL(r.asGenomicRegionMate().pos1,999993); - BOOST_CHECK_EQUAL(r.asGenomicRegionMate().pos2,1000094); - BOOST_CHECK_EQUAL(r.asGenomicRegionMate().strand,'-'); + BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().chr, 22); + BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().pos1,999993); + BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().pos2,1000094); + BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().strand,'-'); BOOST_CHECK_EQUAL(std::floor(r.MeanPhred()), 34); @@ -987,14 +995,6 @@ BOOST_AUTO_TEST_CASE( bam_write ) { // empty constructor SeqLib::BamWriter w; - // specify bam explicitly - //w = SeqLib::BamWriter(SeqLib::BAM); - - //BOOST_CHECK_THROW(w.WriteHeader(), std::runtime_error); - //BOOST_CHECK_THROW(w.Close(), std::runtime_error); - //BOOST_CHECK_THROW(w.BuildIndex(), std::runtime_error); - //BOOST_CHECK_THROW(w.WriteRecord(rec), std::runtime_error); - BOOST_CHECK(!w.WriteHeader()); BOOST_CHECK(!w.Close()); BOOST_CHECK(!w.BuildIndex()); @@ -1004,8 +1004,6 @@ BOOST_AUTO_TEST_CASE( bam_write ) { // check that set CRAM fails BOOST_CHECK(!w.SetCramReference("dummy")); - - //BOOST_CHECK_THROW(w.WriteHeader(), std::runtime_error); BOOST_CHECK(!w.WriteHeader()); w.SetHeader(h); @@ -1014,13 +1012,10 @@ BOOST_AUTO_TEST_CASE( bam_write ) { size_t count = 0; - while(br.GetNextRecord(rec) && count++ < 10000) { + while(br.GetNextRecord(rec) && count++ < 10000) w.WriteRecord(rec); - } - BOOST_CHECK(!w.BuildIndex()); - //BOOST_CHECK_THROW(w.BuildIndex(), std::runtime_error); w.Close(); w.BuildIndex(); @@ -1438,31 +1433,6 @@ BOOST_AUTO_TEST_CASE( plot_test ) { // BOOST_CHECK_EQUAL(rec1.Qname(), rec2.Qname()); // } - -BOOST_AUTO_TEST_CASE( bed_read ) { - - SeqLib::BamReader r; - r.Open("test_data/small.bam"); - - SeqLib::GRC g(BEDFILE, r.Header()); - - BOOST_CHECK_EQUAL(g.size(), 3); - BOOST_CHECK_EQUAL(g[2].chr, 22); - -} - -BOOST_AUTO_TEST_CASE( vcf_read ) { - - SeqLib::BamReader r; - r.Open("test_data/small.bam"); - - SeqLib::GRC g(VCFFILE, r.Header()); - - BOOST_CHECK_EQUAL(g.size(), 31); - BOOST_CHECK_EQUAL(g[29].chr, 22); - -} - BOOST_AUTO_TEST_CASE (json_parse) { SeqLib::BamReader r; diff --git a/seq_test/test_data/test.bed b/seq_test/test_data/test.bed deleted file mode 100644 index 36a36244f..000000000 --- a/seq_test/test_data/test.bed +++ /dev/null @@ -1,3 +0,0 @@ -1 1 10 -2 1 200 -X 4 500 diff --git a/seq_test/test_data/test.vcf b/seq_test/test_data/test.vcf deleted file mode 100644 index eb751de95..000000000 --- a/seq_test/test_data/test.vcf +++ /dev/null @@ -1,151 +0,0 @@ -##fileformat=VCFv4.2 -##fileDate=20160216 -##source=snowman run -t /xchip/gistic/Jeremiah/Data/HCC1143_mem.bam -n /xchip/gistic/Jeremiah/Data/HCC1143_BL_mem.bam -a full -p 15 -D /xchip/gistic/Jeremiah/Projects/SnowmanFilters/dbsnp_138.b37_indel.vcf -##reference=/seq/references/Homo_sapiens_assembly19/v1/Homo_sapiens_assembly19.fasta -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##contig= -##INFO= -##INFO= -##INFO= -##INFO= -##INFO= -##INFO= -##INFO= -##INFO= -##INFO= -##INFO= -##FILTER= -##FILTER= -##FILTER= -##FILTER= -##FILTER= -##FILTER= -##FILTER= -##FILTER= -##FORMAT= -##FORMAT==0.5"> -##FORMAT= -##FORMAT= -##FORMAT= -##FORMAT= -##FORMAT= -##FORMAT= -##FORMAT= -##FORMAT= -##SAMPLE= -##SAMPLE= -#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT /xchip/gistic/Jeremiah/Data/HCC1143_BL_mem.bam /xchip/gistic/Jeremiah/Data/HCC1143_mem.bam -1 4099803 5093 T TGTG 99 PASS LOD=14.38;MAPQ=60;NM=6;READNAMES=x;REPSEQ=AAAAAAAAAA;SCTG=c_1_4091501_4116501_1728C;SOMATIC;SPAN=3 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:66:0:0:0:0:19.87:0:0 0/1:15:70:0:0:9:15:-8.928:14.38:11.5 -1 11439833 12916 AC A 99 PASS LOD=116.1;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAA;SCTG=c_1_11417001_11442001_1156C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:124:0:0:0:0:37.33:0:0 0/1:37:118:0:0:17:37:-112.5:116.1:116.1 -1 14942339 15565 CATTTCTGA C 99 PASS LOD=85.16;MAPQ=60;NM=0;READNAMES=x;REPSEQ=CC;SCTG=c_1_14920501_14945501_1351C;SOMATIC;SPAN=8 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:120:0:0:0:0:36.12:0:0 0/1:28:108:0:0:26:28:-79.49:85.16:85.16 -1 15682315 16005 CTCCCAGGCAAG C 99 PASS LOD=82.81;MAPQ=60;NM=2;READNAMES=x;REPSEQ=AA;SCTG=c_1_15680001_15705001_1508C;SOMATIC;SPAN=11 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:121:0:0:0:0:36.42:0:0 0/1:27:99:0:0:26:27:-78.2:82.81:77.29 -1 38902271 35676 TA T 99 PASS LOD=15.75;MAPQ=60;NM=1;READNAMES=x;REPSEQ=AAAAAAAA;SCTG=c_1_38881501_38906501_937C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:32:0:0:0:0:9.633:0:0 0/1:8:20:0:0:6:8:-15.57:15.75:15.22 -1 39525575 36160 T TA 99 PASS LOD=14.33;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAAAAAAA;SCTG=c_1_39518501_39543501_1148C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:1:87:0:0:1:0:24.79:-0.3299:-0.3299 0/1:25:77:0:0:23:25:-11.77:14.33:14.33 -1 48468636 44071 CT C 99 PASS LOD=20.56;MAPQ=60;NM=1;READNAMES=x;REPSEQ=TTTTTTTTT;SCTG=c_1_48461001_48486001_214C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:1:40:0:0:1:1:9.518:0.5389:0.5209 0/1:14:66:0:0:11:14:-15.45:20.56:19.87 -1 48815931 44181 AT A 99 PASS DBSNP=TRUE;LOD=24.75;MAPQ=60;NM=0;READNAMES=x;REPSEQ=TTTTTTTTTT;SCTG=c_1_48804001_48829001_365C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:66:0:0:0:0:19.87:0:0 0/1:23:84:0:0:10:23:-20.71:24.75:24.75 -1 51821226 46045 G GT 99 PASS DBSNP=TRUE;LOD=77.19;MAPQ=60;NM=1;READNAMES=x;REPSEQ=TTTTTTTT;SCTG=c_1_51817501_51842501_1497C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:56:0:0:0:0:16.86:0:0 0/1:38:85:0:0:28:38:-76.97:77.19:74.62 -1 52315493 46583 C CA 99 PASS LOD=42.12;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAAAAAA;SCTG=c_1_52307501_52332501_1048C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:108:0:0:0:0:32.51:0:0 0/1:38:130:0:0:22:38:-36.87:42.12:42.12 -1 70815825 59243 C CA 99 PASS LOD=38.84;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAAA;SCTG=c_1_70805001_70830001_385C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:84:0:0:0:0:25.29:0:0 0/1:19:73:0:0:14:19:-35.02:38.84:38.84 -1 72689198 60728 CT C 99 PASS LOD=30.29;MAPQ=60;NM=0;READNAMES=x;REPSEQ=TT;SCTG=c_1_72667001_72692001_694C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:81:0:0:0:0:24.38:0:0 0/1:11:77:0:0:11:11:-20.82:30.29:30.29 -1 78512468 64714 GA G 99 PASS LOD=35.05;MAPQ=60;NM=1;READNAMES=x;REPSEQ=AAAAAAAAA;SCTG=c_1_78498001_78523001_575C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:125:0:0:0:0:37.63:0:0 0/1:24:115:0:0:12:24:-25.93:35.05:33.89 -1 80023744 65624 TG T 99 PASS LOD=71.6;MAPQ=60;NM=1;READNAMES=x;REPSEQ=CC;SCTG=c_1_80017001_80042001_1057C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:89:0:0:0:0:26.79:0:0 0/1:24:104:0:0:20:24:-64.69:71.6:69.22 -1 90549816 73294 GA G 99 PASS LOD=78.6;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AA;SCTG=c_1_90527501_90552501_736C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:112:0:0:0:0:33.72:0:0 0/1:26:104:0:0:20:26:-72.69:78.6:78.6 -1 102149804 81384 GTCTGACCAACGTTTGTAATAACATATCCATGAAAAAAAAATATCAAATT G 99 PASS LOD=46.91;MAPQ=60;NM=1;READNAMES=x;SCTG=c_1_102140501_102165501_731C;SOMATIC;SPAN=49 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:112:0:0:0:0:33.72:0:0 0/1:16:77:0:0:16:0:-40.82:46.91:45.35 -1 105990617 83703 TG T 99 PASS LOD=81.94;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAA;SCTG=c_1_105987001_106012001_210C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:108:0:0:0:0:32.51:0:0 0/1:27:60:0:0:23:27:-81.81:81.94:81.94 -1 109175699 85349 GC G 99 PASS LOD=72.08;MAPQ=60;NM=0;READNAMES=x;SCTG=c_1_109172001_109197001_307C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:91:0:0:0:0:27.39:0:0 0/1:23:74:0:0:16:23:-69.72:72.08:72.08 -1 114428673 89982 TA T 99 PASS LOD=48.95;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAA;SCTG=c_1_114415001_114440001_438C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:54:0:0:0:0:16.26:0:0 0/1:19:47:0:0:6:19:-48.57:48.95:48.95 -1 114540407 90128 AT A 99 PASS DBSNP=TRUE;LOD=70.51;MAPQ=60;NM=1;READNAMES=x;REPSEQ=TTTTTTTT;SCTG=c_1_114537501_114562501_366C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:71:0:0:0:0:21.37:0:0 0/1:32:54:0:0:23:32:-70.51:70.51:68.16 -1 118173006 92693 GCCA G 99 PASS LOD=99.78;MAPQ=60;NM=0;READNAMES=x;REPSEQ=CCC;SCTG=c_1_118163501_118188501_700C;SOMATIC;SPAN=3 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:127:0:0:0:0:38.23:0:0 0/1:31:85:0:0:28:31:-98.41:99.78:99.78 -1 154669743 113124 G GT 99 PASS LOD=43.68;MAPQ=60;NM=0;READNAMES=x;REPSEQ=GGGGG;SCTG=c_1_154668501_154693501_1319C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:118:0:0:0:0:35.52:0:0 0/1:17:93:0:0:13:17:-34.89:43.68:43.68 -1 159511118 116958 CA C 99 PASS LOD=84.44;MAPQ=60;NM=1;READNAMES=x;REPSEQ=AAA;SCTG=c_1_159495001_159520001_1269C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:100:0:0:0:0:30.1:0:0 0/1:27:88:0:0:19:27:-81.51:84.44:81.62 -1 165840409 121508 TGC T 99 PASS LOD=9.597;MAPQ=60;NM=0;READNAMES=x;REPSEQ=GTGTGT;SCTG=c_1_165816001_165841001_1416C;SOMATIC;SPAN=2 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:128:0:0:0:0:38.53:0:0 0/1:6:115:0:0:3:6:14.81:9.597:9.597 -1 170974122 125308 AT A 99 PASS DBSNP=TRUE;LOD=88.79;MAPQ=60;NM=0;READNAMES=x;REPSEQ=TTTTTTTTTT;SCTG=c_1_170961001_170986001_1029C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:1:91:0:0:1:1:25.39:-0.1145:-0.1145 0/1:60:105:0:0:37:60:-88.79:88.79:88.79 -1 174182428 127789 TG T 99 PASS LOD=14.34;MAPQ=60;NM=0;READNAMES=x;REPSEQ=GG;SCTG=c_1_174170501_174195501_1198C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:108:0:0:0:0:32.51:0:0 0/1:6:93:0:0:6:6:3.996:14.34:14.34 -1 176466739 128990 C CA 99 PASS LOD=22.2;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAAAAA;SCTG=c_1_176449001_176474001_1047C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:68:0:0:0:0:20.47:0:0 0/1:15:69:0:0:15:14:-17.07:22.2:22.2 -3 181761320 133129 GA G 99 PASS DBSNP=TRUE;LOD=236.6;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAA;SCTG=c_1_181741001_181766001_1395C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:113:0:0:0:0:34.02:0:0 0/1:83:131:0:0:62:83:-236.6:236.6:236.6 -Y 184963190 135220 TA T 99 PASS LOD=9.293;MAPQ=60;NM=0;READNAMES=x;REPSEQ=TTTTTTTTTT;SCTG=c_1_184950501_184975501_1050C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:102:0:0:0:0:30.71:0:0 0/1:12:84:0:0:8:12:1.287:9.293:9.293 -X 185540320 135702 GA G 99 PASS DBSNP=TRUE;LOD=58.5;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAA;SCTG=c_1_185538501_185563501_562C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:131:0:0:0:0:39.43:0:0 0/1:25:97:0:0:21:25:-53.33:58.5:58.5 -X 1100000 135702 GA G 99 PASS DBSNP=TRUE;LOD=58.5;MAPQ=60;NM=0;READNAMES=x;REPSEQ=AAAAAA;SCTG=c_1_185538501_185563501_562C;SOMATIC;SPAN=1 GT:AD:DP:GQ:PL:SR:CR:LR:LO:SL 0/1:0:131:0:0:0:0:39.43:0:0 0/1:25:97:0:0:21:25:-53.33:58.5:58.5 diff --git a/src/BamWriter.cpp b/src/BamWriter.cpp index bfab0cc27..b2594be8c 100644 --- a/src/BamWriter.cpp +++ b/src/BamWriter.cpp @@ -133,8 +133,14 @@ bool BamWriter::SetCramReference(const std::string& ref) { // need to open reference for CRAM writing char* fn_list = samfaipath(ref.c_str()); // eg ref = my.fa returns my.fa.fai if (fn_list) { - int status = hts_set_fai_filename(fop.get(), fn_list); - if (status != 0) { + + // in theory hts_set_fai_filename should give back < 0 + // if fn_list not there, but it doesnt + if (!read_access_test(std::string(fn_list))) + return false; + + int status = hts_set_fai_filename(fop.get(), fn_list); + if (status < 0) { fprintf(stderr, "Failed to use reference \"%s\".\n", fn_list); return false; }