Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/allele_functions.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,8 @@ std::vector<allele> update_allele_depth(char ref,std::string bases, std::string
}
if(ad.size() > 0)
std::sort(ad.begin(), ad.end());
//test code
//print_allele_depths(ad);
return ad;
}

Expand Down
15 changes: 13 additions & 2 deletions src/call_variants.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min
<< std::endl;
int ctr = 0;
int64_t pos = 0;
uint32_t mdepth = 0, pdepth = 0; // mpdepth for mpileup depth and pdeth for ungapped depth at position
uint32_t mdepth = 0, pdepth = 0, gdepth =0; // mpdepth for mpileup depth and pdeth for ungapped depth at position
double pval_left, pval_right, pval_twotailed, *freq_depth, err;
std::stringstream line_stream;
char ref;
Expand Down Expand Up @@ -88,14 +88,19 @@ int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min
}
ctr++;
}

ad = update_allele_depth(ref, bases, qualities, min_qual);
if(ad.size() == 0){
line_stream.clear();
continue;
}
// Get ungapped depth
pdepth = 0;
//get gapped depth
gdepth = 0;
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
gdepth += it->depth;
//commenting these two lines gets the gapped depth
if(it->nuc[0]=='*' || it->nuc[0] == '+' || it->nuc[0] == '-')
continue;
pdepth += it->depth;
Expand All @@ -117,8 +122,14 @@ int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
if((*it == *ref_it) || it->nuc[0]=='*')
continue;
freq_depth = get_frequency_depth(*it, pdepth, mdepth);
//use gapped depth if we're looking at insertion or deletion
if (it->nuc[0] == '+' || it->nuc[0] == '-'){
freq_depth = get_frequency_depth(*it, gdepth, gdepth);
}else{
freq_depth = get_frequency_depth(*it, pdepth, mdepth);
}
if(freq_depth[0] < min_threshold)

continue;
//this only adds the first bit of the information to the tsv
out_str << region << "\t";
Expand Down