allele_functions.cpp
#include<string>
#include<vector>
#include<algorithm>
#include<iostream>
#include "allele_functions.h"
void print_allele_depths(std::vector<allele> ad){
std::cout << "AD Size: " << ad.size() << " " << std::endl;
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
std::cout << "Nuc: " << it->nuc << std::endl;;
std::cout << "Depth: " << it->depth << std::endl;
std::cout << "Reverse: " << it->reverse <<std::endl;
std::cout << "Qual: " << (uint16_t) it->mean_qual << std::endl;
std::cout << "Beg: " << it->beg << std::endl;
std::cout << "End: " << it->end << std::endl;
}
}
int check_allele_exists(std::string n, std::vector<allele> ad){
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
if(it->nuc.compare(n) == 0){
return it - ad.begin();
}
}
return -1;
}
int find_ref_in_allele(std::vector<allele> ad, char ref){
std::string ref_s(1, ref);
std::vector<allele>::iterator it = ad.begin();
while(it < ad.end()){
if(it->nuc.compare(ref_s) == 0)
return (it - ad.begin());
it++;
}
return -1;
}
std::vector<allele> update_allele_depth(char ref,std::string bases, std::string qualities, uint8_t min_qual){
std::vector<allele> ad;
std::string indel;
int i = 0, n =0, j = 0, q_ind = 0;
bool beg, end;
uint8_t q;
float num = 0, den = 0;
while (i < bases.length()){
beg = false;
end = false;
if(bases[i] == '^'){
i += 2; // Skip mapping quality as well (i+1) - 33
continue;
}
if(bases[i] == '$'){
i++;
continue;
}
q = qualities[q_ind] - 33;
std::string b;
allele tmp;
bool forward= true;
switch(bases[i]){
case '.':
b = ref;
end = (bases[i+1] == '$');
beg = (bases[i+1] == '^');
break;
case ',':
b = ref;
forward = false;
end = (bases[i+1] == '$');
beg = (bases[i+1] == '^');
break;
case '*':
b = bases[i];
break;
case '+': case '-':
j = i+1;
while(isdigit(bases[j])){
j++;
}
j = j - (i+1);
n = stoi(bases.substr(i+1, j));
indel = bases.substr(i+1+j, n);
transform(indel.begin(), indel.end(), indel.begin(),::toupper);
b = bases[i] + indel; // + for Insertion and - for Deletion
i += n + j;
if(indel[0]>=97 && indel[0] <= 122)
forward=false;
q = min_qual; // For insertions and deletion ust use minimum quality.
break;
default:
int asc_val = bases[i];
if(asc_val >= 65 && asc_val <= 90){
b = bases[i];
} else if(bases[i]>=97 && bases[i]<=122){
b = bases[i] - ('a' - 'A');
forward = false;
}
end = (bases[i+1] == '$');
beg = (bases[i+1] == '^');
}
int ind = check_allele_exists(b, ad);
if(q >= min_qual){
if (ind==-1){
tmp.nuc = b;
tmp.depth = 1;
tmp.tmp_mean_qual = q;
if(!forward)
tmp.reverse = 1;
else
tmp.reverse = 0;
if(beg)
tmp.beg = 1;
else
tmp.beg = 0;
if(end)
tmp.end = 1;
else
tmp.end = 0;
ad.push_back(tmp);
} else {
ad.at(ind).tmp_mean_qual = (ad.at(ind).tmp_mean_qual * ad.at(ind).depth + q)/(ad.at(ind).depth + 1);
ad.at(ind).depth += 1;
if(beg)
ad.at(ind).beg += 1;
if(end)
ad.at(ind).end += 1;
if(!forward)
ad.at(ind).reverse += 1;
}
}
i++;
if(b[0] !='+' && b[0]!='-')
q_ind++;
}
for(std::vector<allele>::iterator it = ad.begin(); it!=ad.end(); ++it){
it->mean_qual = (uint8_t) it->tmp_mean_qual;
}
if(ad.size() > 0)
std::sort(ad.begin(), ad.end());
return ad;
}
int get_index(char a){
switch(a){
case'Y':
a = 0;
break;
case'R':
a = 1;
break;
case'W':
a = 2;
break;
case'S':
a = 3;
break;
case'K':
a = 4;
break;
case'M':
a = 5;
break;
case'A':
a = 6;
break;
case'T':
a = 7;
break;
case'G':
a = 8;
break;
case'C':
a = 9;
break;
case'D':
a = 10;
break;
case'V':
a = 11;
break;
case'H':
a = 12;
break;
case'B':
a = 13;
break;
default:
a = -1;
}
return (int) a;
}
// Modified gt2iupac from bcftools.h. Expanded iupac matrix to include all ambigious nucleotides(added Y, R, W, S, K, M, A, T, G, C, D, V, H, B) - https://github.com/samtools/bcftools/blob/b0376dff1ed70603c9490802f37883b9009215d2/bcftools.h#L48
/*
Expanded IUPAC matrix:
Y N H B B H H Y B Y N N H B
N R D V D V R D R V D V N N
H D W N D H W W D H D N H N
B V N S B V V B S S N V N B
B D D B K N D K K B D N N B
H V H V N M M H V M N V H N
H R W V D M A W R M D V H N
Y D W B K H W T K Y D N H B
B R D S K V R K G S D V N B
Y V H S B M M Y S C N V H B
N D D N D N D D D N D N N N
N V N V N V V N V V N V N N
H N H N N H H H N H N N H N
B N N B B N N B B B N N N B
*/
char gt2iupac(char a, char b){
if(a == '*' || b == '*')
return 'N';
static const char iupac[14][14] = {
{'Y', 'N', 'H', 'B', 'B', 'H', 'H', 'Y', 'B', 'Y', 'N', 'N', 'H', 'B'},
{'N', 'R', 'D', 'V', 'D', 'V', 'R', 'D', 'R', 'V', 'D', 'V', 'N', 'N'},
{'H', 'D', 'W', 'N', 'D', 'H', 'W', 'W', 'D', 'H', 'D', 'N', 'H', 'N'},
{'B', 'V', 'N', 'S', 'B', 'V', 'V', 'B', 'S', 'S', 'N', 'V', 'N', 'B'},
{'B', 'D', 'D', 'B', 'K', 'N', 'D', 'K', 'K', 'B', 'D', 'N', 'N', 'B'},
{'H', 'V', 'H', 'V', 'N', 'M', 'M', 'H', 'V', 'M', 'N', 'V', 'H', 'N'},
{'H', 'R', 'W', 'V', 'D', 'M', 'A', 'W', 'R', 'M', 'D', 'V', 'H', 'N'},
{'Y', 'D', 'W', 'B', 'K', 'H', 'W', 'T', 'K', 'Y', 'D', 'N', 'H', 'B'},
{'B', 'R', 'D', 'S', 'K', 'V', 'R', 'K', 'G', 'S', 'D', 'V', 'N', 'B'},
{'Y', 'V', 'H', 'S', 'B', 'M', 'M', 'Y', 'S', 'C', 'N', 'V', 'H', 'B'},
{'N', 'D', 'D', 'N', 'D', 'N', 'D', 'D', 'D', 'N', 'D', 'N', 'N', 'N'},
{'N', 'V', 'N', 'V', 'N', 'V', 'V', 'N', 'V', 'V', 'N', 'V', 'N', 'N'},
{'H', 'N', 'H', 'N', 'N', 'H', 'H', 'H', 'N', 'H', 'N', 'N', 'H', 'N'},
{'B', 'N', 'N', 'B', 'B', 'N', 'N', 'B', 'B', 'B', 'N', 'N', 'N', 'B'}
};
if ( a>='a' ) a -= 'a' - 'A';
if ( b>='a' ) b -= 'a' - 'A';
int _a, _b;
_a = get_index(a);
_b = get_index(b);
if(_a == -1 || _b == -1)
return 'N';
return iupac[_a][_b];
}
allele
Definition: allele_functions.h:7