call_variants.cpp
#include<iostream>
#include<fstream>
#include<sstream>
#include<vector>
#include<algorithm>
#include<string>
#include<regex>
#include<cmath>
#include<htslib/kfunc.h>
#include "htslib/kfunc.h"
#include "allele_functions.h"
const char gap='N';
const float sig_level = 0.01;
std::vector<allele>::iterator get_ref_allele(std::vector<allele> &ad, char ref){
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
if(it->nuc[0] == ref)
return it;
}
// print_allele_depths(ad);
return ad.end();
}
int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min_qual, double min_threshold){
std::string line, cell, bases, qualities, region;
std::ofstream fout(out_file+".tsv");
fout << "REGION\tPOS\tREF\tALT\tAD\tRAD\tDP\tFREQ\tQUAL\tPVAL\tPASS"<<std::endl;
int ctr = 0, pos = 0, mdepth = 0;
double pval_left, pval_right, pval_twotailed, freq, err;
std::stringstream lineStream;
char ref;
std::vector<allele> ad;
std::vector<allele>::iterator ref_it;
while (std::getline(cin, line)){
lineStream << line;
ctr = 0;
while(std::getline(lineStream,cell,'\t')){
switch(ctr){
case 0:
region = cell;
break;
case 1:
pos = stoi(cell);
break;
case 2:
ref = cell[0];
break;
case 3:
mdepth = stoi(cell);
break;
case 4:
bases = cell;
break;
case 5:
qualities = cell;
break;
case 6:
break;
}
ctr++;
}
ad = update_allele_depth(ref, bases, qualities, min_qual);
if(ad.size() == 0){
lineStream.clear();
continue;
}
ref_it = get_ref_allele(ad, ref);
// Get ungapped coverage
mdepth = 0;
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
if(it->nuc[0]=='*' || it->nuc[0] == '+' || it->nuc[0] == '-')
continue;
mdepth += it->depth;
}
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
if((*it == *ref_it) || it->nuc[0]=='*')
continue;
freq = it->depth/(double)mdepth;
if(freq < min_threshold)
continue;
fout << region << "\t";
fout << pos << "\t";
fout << ref << "\t";
fout << it->nuc << "\t";
fout << it->depth << "\t";
fout << it->reverse << "\t";
fout << mdepth << "\t";
fout << freq << "\t";
fout << (uint16_t)it->mean_qual << "\t";
/*
| Var | Ref |
Exp | Error | Err free |
Obs | AD | RD |
*/
err = pow(10, ( -1 * (it->mean_qual)/10));
kt_fisher_exact((err * mdepth), (1-err) * mdepth, it->depth, ref_it->depth, &pval_left, &pval_right, &pval_twotailed);
fout << pval_left << "\t";
if(pval_left <= sig_level){
fout << "TRUE";
} else {
fout << "FALSE";
}
fout << std::endl;
}
lineStream.clear();
}
fout.close();
return 0;
}
// int main_old(int argc, char* argv[]) {
// std::string line, cell;
// std::string out_file = argv[1];
// uint8_t min_qual = 20;
// if(argc > 2)
// min_qual = atoi(argv[2]);
// std::cout << "Min Qual: " << (uint16_t)min_qual << std::endl;
// std::ofstream fout(out_file+".tsv");
// fout << "POS\tREF\tALT\tAD\tRAD\tDP\tQUAL"<<std::endl;
// int ctr = 0, pos = 0, mdepth = 0;
// std::stringstream lineStream;
// char ref;
// std::string bases;
// std::string qualities;
// std::vector<allele> ad;
// while (std::getline(std::cin, line)){
// lineStream << line;
// ctr = 0;
// while(std::getline(lineStream,cell,'\t')){
// switch(ctr){
// case 0:
// break;
// case 1:
// pos = stoi(cell);
// break;
// case 2:
// ref = cell[0];
// break;
// case 3:
// mdepth = stoi(cell);
// break;
// case 4:
// bases = cell;
// break;
// case 5:
// qualities = cell;
// break;
// case 6:
// break;
// }
// ctr++;
// }
// ad = update_allele_depth(ref, bases, qualities, min_qual);
// for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
// if((it->nuc[0] == ref && it->nuc.length() == 1) || it->nuc[0]=='*')
// continue;
// fout << pos << "\t";
// fout << ref << "\t";
// fout << it->nuc << "\t";
// fout << it->depth << "\t";
// fout << it->reverse << "\t";
// fout << mdepth << "\t";
// fout << (uint16_t)it->mean_qual << std::endl;
// }
// lineStream.clear();
// }
// fout.close();
// return 0;
// }