Skip to content

Commit

Permalink
sga-graphdiff now directly outputs the final set of calls.
Browse files Browse the repository at this point in the history
  • Loading branch information
@jts
jts committed Jun 29, 2012
1 parent 3b8915c commit 21daaeb
Show file tree
Hide file tree
Showing 9 changed files with 79 additions and 15 deletions.
25 changes: 19 additions & 6 deletions src/GraphDiff/DindelRealignWindow.cpp
View file
Expand Up @@ -1435,24 +1435,37 @@ void DindelRealignWindowResult::Inference::outputAsVCF(const DindelVariant & var
record.addComment("SB", this->strandBias);
record.addComment("NumLib", this->numLibraries);
record.addComment("NumFragments", this->numReadNames);
out.push_back(record);

/*
// HistDist
std::stringstream histdist_ss;
std::copy(histDistance.begin(), histDistance.end(), std::ostream_iterator<int>(histdist_ss, ","));
for (size_t x=0;x<histDistance.size();x++) {
if (x>0)
histdist_ss << ",";
histdist_ss << histDistance[x];
}
record.addComment("HistDist", histdist_ss.str());

// HistLik
std::stringstream histlik_ss;
std::copy(histAlignLik.begin(), histAlignLik.end(), std::ostream_iterator<int>(histdist_ss, ","));
for (size_t x=0;x<histAlignLik.size();x++) {
if (x>0)
histlik_ss << ",";
histlik_ss << histAlignLik[x];
}
record.addComment("HistLik", histlik_ss.str());

// HistMapQ
std::stringstream histmap_ss;
std::copy(histMapQ.begin(), histMapQ.end(), std::ostream_iterator<int>(histMapQ_ss, ","));
record.addComment("HistMapQ", histlik_ss.str());
for (size_t x=0;x<histMapQ.size();x++) {
if (x>0)
histmap_ss << ",";
histmap_ss << histMapQ[x];
}
record.addComment("HistMapQ", histmap_ss.str());

out.push_back(record);

/*
out.precision(5);
out.setf(std::ios::fixed,std::ios::floatfield);
out << var.getChrom() << "\t" << var.getPos() << "\t";
Expand Down
20 changes: 18 additions & 2 deletions src/GraphDiff/DindelUtil.cpp
View file
Expand Up @@ -21,7 +21,8 @@ DindelReturnCode DindelUtil::runDindelPairMatePair(const std::string& id,
const StringVector& variant_haplotypes,
const GraphCompareParameters& parameters,
std::ostream& baseOut,
std::ostream& variantOut)
std::ostream& variantOut,
std::ostream& callsOut)
{
PROFILE_FUNC("runDindelPairMatePair")

Expand Down Expand Up @@ -283,14 +284,29 @@ DindelReturnCode DindelUtil::runDindelPairMatePair(const std::string& id,
pPreviousResult = pThisResult;
}

// Copy VCFRecords to output
// Copy raw VCFRecords to output
for(size_t i = 0; i <= 1; ++i)
{
std::ostream& curr_out = i == 0 ? baseOut : variantOut;
for(size_t j = 0; j < vcfRecords[i].size(); ++j)
curr_out << vcfRecords[i][j] << "\n";
}

// Make comparative calls
size_t VARIANT_IDX = 1;
size_t BASE_IDX = 0;
bool has_base_calls = !vcfRecords[BASE_IDX].empty();
for(size_t i = 0; i < vcfRecords[1].size(); ++i)
{
bool not_called_in_base = true;
if(has_base_calls)
not_called_in_base = vcfRecords[BASE_IDX][i].passStr == "NoCall";

bool called_in_variant = vcfRecords[VARIANT_IDX][i].passStr == "PASS";
if(called_in_variant && not_called_in_base)
callsOut << vcfRecords[VARIANT_IDX][i] << "\n";
}

baseOut.flush();
variantOut.flush();

Expand Down
3 changes: 2 additions & 1 deletion src/GraphDiff/DindelUtil.h
View file
Expand Up @@ -36,7 +36,8 @@ namespace DindelUtil
const StringVector& variant_haplotypes,
const GraphCompareParameters& parameters,
std::ostream& baseOut,
std::ostream& variantOut);
std::ostream& variantOut,
std::ostream& callsOut);

// Run a naive caller
DindelReturnCode runNaiveCaller(const std::string& normalString,
Expand Down
15 changes: 13 additions & 2 deletions src/GraphDiff/GraphCompare.cpp
View file
Expand Up @@ -175,12 +175,14 @@ GraphCompareResult GraphCompare::process(const SequenceWorkItem& item)

std::stringstream baseVCFSS;
std::stringstream variantVCFSS;
std::stringstream callsVCFSS;
DindelReturnCode drc = DindelUtil::runDindelPairMatePair(kmer,
build_result.base_haplotypes,
build_result.variant_haplotypes,
m_parameters,
baseVCFSS,
variantVCFSS);
variantVCFSS,
callsVCFSS);

//
if(m_parameters.verbose > 0 || 1)
Expand All @@ -195,6 +197,7 @@ GraphCompareResult GraphCompare::process(const SequenceWorkItem& item)
{
result.baseVCFStrings.push_back(baseVCFSS.str());
result.variantVCFStrings.push_back(variantVCFSS.str());
result.calledVCFStrings.push_back(callsVCFSS.str());
}
}
}
Expand Down Expand Up @@ -498,13 +501,15 @@ void GraphCompare::testKmer(const std::string& kmer)

std::stringstream baseVCFSS;
std::stringstream variantVCFSS;
std::stringstream callsVCFSS;

DindelReturnCode drc = DindelUtil::runDindelPairMatePair(kmer,
build_result.base_haplotypes,
build_result.variant_haplotypes,
m_parameters,
baseVCFSS,
variantVCFSS);
variantVCFSS,
callsVCFSS);

std::cout << "base: " << baseVCFSS.str() << "\n";
std::cout << "variant: " << variantVCFSS.str() << "\n";
Expand Down Expand Up @@ -533,6 +538,7 @@ void GraphCompare::showMappingLocations(const std::string& str)
//
GraphCompareAggregateResults::GraphCompareAggregateResults(const std::string& fileprefix) : m_baseVCFFile(fileprefix + ".base.vcf","w"),
m_variantVCFFile(fileprefix + ".variant.vcf","w"),
m_callsVCFFile(fileprefix + ".calls.vcf","w"),
m_numVariants(0)
{
//
Expand All @@ -541,6 +547,7 @@ GraphCompareAggregateResults::GraphCompareAggregateResults(const std::string& fi
//
m_baseVCFFile.outputHeader("stub", "stub");
m_variantVCFFile.outputHeader("stub", "stub");
m_callsVCFFile.outputHeader("stub", "stub");

// Initialize mutex
int ret = pthread_mutex_init(&m_mutex, NULL);
Expand Down Expand Up @@ -601,6 +608,10 @@ void GraphCompareAggregateResults::process(const SequenceWorkItem& /*item*/, con
m_baseVCFFile.getOutputStream() << result.baseVCFStrings[i];
m_variantVCFFile.getOutputStream() << result.variantVCFStrings[i];
}

// Write out the final calls
for(size_t i = 0; i < result.calledVCFStrings.size(); ++i)
m_callsVCFFile.getOutputStream() << result.calledVCFStrings[i];
}

//
Expand Down
2 changes: 2 additions & 0 deletions src/GraphDiff/GraphCompare.h
View file
Expand Up @@ -82,6 +82,7 @@ struct GraphCompareResult

StringVector baseVCFStrings;
StringVector variantVCFStrings;
StringVector calledVCFStrings;
};

//
Expand Down Expand Up @@ -204,6 +205,7 @@ class GraphCompareAggregateResults
// VCF output files
VCFFile m_baseVCFFile;
VCFFile m_variantVCFFile;
VCFFile m_callsVCFFile;

size_t m_numVariants;
};
Expand Down
8 changes: 6 additions & 2 deletions src/GraphDiff/VCFTester.cpp
View file
Expand Up @@ -19,10 +19,12 @@
VCFTester::VCFTester(const GraphCompareParameters& params) : m_parameters(params),
m_graphComparer(params),
m_baseVCFFile("vtest.base.vcf","w"),
m_variantVCFFile("vtest.variant.vcf","w")
m_variantVCFFile("vtest.variant.vcf","w"),
m_callsVCFFile("vtest.calls.vcf","w")
{
m_baseVCFFile.outputHeader("stub", "stub");
m_variantVCFFile.outputHeader("stub", "stub");
m_callsVCFFile.outputHeader("stub", "stub");
DindelUtil::initializeCodeCounts(m_returnCodes);
}

Expand Down Expand Up @@ -63,6 +65,7 @@ void VCFTester::process(const VCFFile::VCFEntry& record)
// Run dindel
std::stringstream baseSS;
std::stringstream variantSS;
std::stringstream callsSS;
std::string id = ".";

StringVector ref_haplotypes;
Expand All @@ -75,7 +78,8 @@ void VCFTester::process(const VCFFile::VCFEntry& record)
var_haplotypes,
m_parameters,
baseSS,
variantSS);
variantSS,
callsSS);

m_returnCodes[code] += 1;

Expand Down
1 change: 1 addition & 0 deletions src/GraphDiff/VCFTester.h
View file
Expand Up @@ -49,6 +49,7 @@ class VCFTester
// Output file
VCFFile m_baseVCFFile;
VCFFile m_variantVCFFile;
VCFFile m_callsVCFFile;
};

#endif
19 changes: 17 additions & 2 deletions src/Util/VCFUtil.cpp
View file
Expand Up @@ -40,6 +40,16 @@ void VCFRecord::addComment(const std::string& key, const int& value)
addComment(key, v_str.str());
}

void VCFRecord::addComment(const std::string& key, const double& value)
{
std::stringstream v_str;
v_str.precision(5);
v_str.setf(std::ios::fixed,std::ios::floatfield);
v_str << value;
addComment(key, v_str.str());
}


void VCFRecord::setPassStr(const std::string str)
{
passStr = str;
Expand Down Expand Up @@ -71,8 +81,13 @@ std::ostream& operator<<(std::ostream& o, VCFRecord& record)
o << record.varStr << "\t";
o << (int)record.quality << "\t";
o << record.passStr << "\t";
std::copy( record.comments.begin(), record.comments.end(),
std::ostream_iterator<std::string>(o, ";")); // comments

for(size_t i = 0; i < record.comments.size(); ++i)
{
if(i > 0)
o << ";";
o << record.comments[i];
}
return o;
}

Expand Down
1 change: 1 addition & 0 deletions src/Util/VCFUtil.h
View file
Expand Up @@ -42,6 +42,7 @@ struct VCFRecord
VCFRecord();
void addComment(const std::string& key, const std::string& value);
void addComment(const std::string& key, const int& value);
void addComment(const std::string& key, const double& value);

void setPassStr(const std::string str);
void setQuality(const double q);
Expand Down

0 comments on commit 21daaeb

Please sign in to comment.