alignment.cpp
#include<iostream>
#include <algorithm>
#include<vector>
#include<fstream>
#include "alignment.h"
/* Substitution Matrix
A T G C N
A 1 -1 -1 -1 0
T -1 1 -1 -1 0
G -1 -1 1 -1 0
C -1 -1 -1 1 0
N 0 0 0 0 0
*/
int get_sub_score(char a, char b){
int i,j;
switch(a){
case 'A':
i = 0;
break;
case 'T':
i = 1;
break;
case 'G':
i = 2;
break;
case 'C':
i = 3;
break;
case 'N':
i = 4;
break;
}
switch(b){
case 'A':
j = 0;
break;
case 'T':
j = 1;
break;
case 'G':
j = 2;
break;
case 'C':
j = 3;
break;
case 'N':
i = 4;
break;
}
return substitution[i][j];
}
int get_gap_penalty(int k, char a){
if(a == 'N') // Allow for i7 and i5 indexesto be specified as NNNNNNNN
return 0;
return k * gap_open; // Linear Gap Penalty
}
void print_matrix(int h[][max_adapter_size], int r, int c, std::string read, std::string adap){
for(int i =0; i< r;i++){
if(i == 0){
std::cout << " ";
for(int k =0; k < c-1; k++) std::cout << adap[k] << " ";
std::cout << "\n";
}
for(int j = 0; j < c;j++){
if(j == 0 && i > 0){
if(i -1 < read.length())
std::cout << read[i-1] << " ";
}
if(i == 0 && j == 0)
std::cout << " " << " ";
std::cout << h[i][j] << " ";
}
std::cout << "\n";
}
}
// Return Value 1 - Diag, 2->Right, 3->Down, 0-> 0
int* get_score_cell(int h[][max_adapter_size], int i, int j, std::string read, std::string adap){
int s[] = {0,0,0};
int *tmp = new int[2];
tmp[0] = 0;
tmp[1] = 0;
s[0] = h[i-1][j-1] + get_sub_score(read[i-1], adap[j-1]);
s[1] = h[i-1][j] - get_gap_penalty(1, adap[j-1]);
s[2] = h[i][j-1] - get_gap_penalty(1, read[i-1]);
int max = (s[0]>=s[1]) ? s[0] : s[1];
max = (max>=s[2]) ? max : s[2];
max = (max > 0) ? max : 0;
tmp[0] = max;
if(max == s[0]){
tmp[1] = 1;
} else if(max == s[1]){
tmp[1] = 3;
} else if(max == s[2]){
tmp[1] = 2;
}
return tmp;
}
void print_alignment(char a[2][max_read_size], int n){
std::cout << "Alignment: " << std::endl << std::endl;
for(int j =0;j<2;j++){
if(j == 0)
std::cout << "Read: ";
else if(j == 1)
std::cout << "Adap: ";
for(int i = n -1; i >= 0;i--){
std::cout << a[j][i] << " ";
}
std::cout << std::endl;
}
}
int* align_seqs(std::string read, std::string adap){
int h[max_read_size][max_adapter_size],
t[max_read_size][max_adapter_size],
m = read.length() + 1,
n=adap.length()+ 1,
max_i = 0,
max_j = 0,
max_v = -1,
max_score,
*rt = new int[2],
*tmp = new int[2];
rt[0] = -1;
rt[1] = read.length();
// Initialize first column to 0
for(int i = 0; i < n; i++) h[0][i] = 0;
// Initialize first row to 0
for(int i = 0; i < m; i++) h[i][0] = 0;
// print_matrix(h, m, n, read, adap);
for(int i = 1; i < m; i++){
for(int j = 1; j < n; j++){
tmp = get_score_cell(h, i, j, read, adap);
h[i][j] = *tmp;
t[i][j] = *(tmp+1);
if(*tmp >= max_v){
max_i = i;
max_j = j;
max_v = *tmp;
}
}
}
max_score = max_v;
// if(max_i != read.length()) // Ensure alignment to end of read.
// return rt;
// std::cout << "Score: " << max_v << " ";
// print_matrix(h, m, n, read, adap);
// print_matrix(t, m, n, read, adap);
// Traceback
int _t, _l, _d, _align_n = 0;
char _align[2][max_read_size];
// 1 - Diag, 2->Right, 3->Down
while(max_v != 0){
_d = h[max_i-1][max_j-1];
_l = h[max_i][max_j-1];
_t = h[max_i-1][max_j];
switch(t[max_i][max_j]){
case 1:
max_v = _d;
max_i = max_i -1;
max_j = max_j - 1;
_align[0][_align_n] = read[max_i];
_align[1][_align_n] = adap[max_j];
_align_n++;
break;
case 2:
max_v = _l;
max_j = max_j - 1;
_align[0][_align_n] = '-';
_align[1][_align_n] = adap[max_j];
_align_n++;
break;
case 3:
max_v = _t;
max_i = max_i - 1;
_align[0][_align_n] = read[max_i];
_align[1][_align_n] = '-';
_align_n++;
break;
case 0:
max_v = 0;
break;
}
}
rt[1] = max_i;
print_alignment(_align, _align_n);
if(max_score <= ((_align_n - 2) * unit_score) - get_gap_penalty(2, 'A') || max_j != 0) // If alignment does not start at beginning of adapter
return rt;
rt[0] = max_score;
// print_alignment(_align, _align_n);
return rt;
}
// int find_adapters_contaminants(std::istream &cin, std::string adp_cntms_file){
// std::vector<std::string> adp = read_adapters_from_fasta(adp_cntms_file, "Read1");
// std::string line, bases, qual, b;
// int ctr = 0, score, *t = new int[2];
// while (std::getline(cin, line)){
// if(ctr % 400000 == 0)
// std::cout << "Processed " << (ctr/4) << " reads ..." << std::endl;
// if(ctr % 2 == 0){
// ctr++;
// continue;
// }
// if((ctr+1) % 2 == 0 && (ctr + 1) % 4 == 0 ){
// qual = line;
// // std::cout << bases << std::endl;
// // std::cout << qual << std::endl;
// score = -1;
// for(std::vector<std::string>::iterator it = adp.begin(); it != adp.end() && score == -1; ++it) {
// b = ((*it).length() > bases.length()) ? bases : bases.substr(bases.length() - (*it).length(), (*it).length());
// t = align_seqs(b, *it);
// score = *t;
// if(score != -1)
// std::cout << "Trimmed: " << bases.substr(0, bases.length() - *(t+1)) << std::endl;
// }
// } else {
// bases = line;
// }
// ctr++;
// }
// return 0;
// }
// int main(int argc, char* argv[]){
// std::string adp = "../data/adapters/NexteraPR.fa";
// // clock_t begin = clock();
// find_adapters_contaminants(std::cin, adp);
// // clock_t end = clock();
// // double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
// // std::cout << "Time Taken: " << elapsed_secs << std::endl;
// return 0;
// }