/usr/include/blasr/algorithms/alignment/FullQVAlign.hpp is in libblasr-dev 0~20151014+gitbe5d1bf-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 | #ifndef _BLASR_FULL_QV_ALIGN_HPP_
#define _BLASR_FULL_QV_ALIGN_HPP_
#include "matrix/Matrix.hpp"
#include "FASTQSequence.hpp"
#include "FASTASequence.hpp"
template<typename T_Query, typename T_Reference>
double FullQVAlign(T_Query &query,
T_Reference &target,
Matrix<double> &alignProb) {
alignProb.Resize(query.length + 1, target.length + 1);
alignProb.Initialize(0);
DNALength q, t;
if (query.length == 0 or target.length == 0) { return 0; }
// Initialize boundaries of ins/del/match probability matrices.
q = 0;
VectorIndex numCols = target.length + 1;
VectorIndex numRows = query.length + 1;
alignProb[0][0] = 1;
for (t = 1; t < numCols; t++ ) {
// cannot match to a gap
alignProb[0][t] = log(target.GetInsertionQV(t-1)) + alignProb[0][t-1];
}
for (q = 1; q < numRows; q++) {
alignProb[q][0] = log(query.GetInsertionQV(q-1)) + alignProb[q-1][0];
}
// Now compute probability of alignment with the Forward algorithm.
for (q = 1; q < numRows; q++) {
for (t = 1; t < numCols; t++) {
// First compute p_ins[q,t] as transitions from match matrix
double logMatchedPulseProb = 0;
double logInsertedPulseProb = 0;
double logDeletedPulseProb = 0;
//
// Use inefficient coding for now.
//
// Compute match, the bases are either the same, in which case
// this is simply the product of the probabilities of the two
// matches. Otherwise, either one of the pulses may be correct,
// and the probability is the union of the two cases.
//
double matchedPulseProb;
if (query.seq[q-1] == target.seq[t-1]) {
matchedPulseProb = (1-query.GetSubstitutionQV(q-1)) * (1-target.GetSubstitutionQV(t-1));
}
else {
matchedPulseProb = (query.GetSubstitutionQV(q-1)/3.0)*(1-target.GetSubstitutionQV(t-1)) +
((1-query.GetSubstitutionQV(q-1)))*(target.GetSubstitutionQV(t-1)/3.0);
}
matchedPulseProb = exp(alignProb[q-1][t-1])*matchedPulseProb;
//
// An insertion in the query can be either a normal extra base
// in the query, or a deletion in the reference.
//
// logInsertedPulseProb = uery.GetInsertionQV(q-1)) + alignProb[q-1][t];
double insertedPulseProb = 0;
if (target.GetDeletionTag(t-1) != 'N') {
//
// The target has a pulse that was not strong enough to call a
// real incorporation. For now assume that the weak pulse is
// the previous nucleotide in the query. So the likelihood of
// the weak pulse is influenced by the likelihood of the
// previous nucleotide in the query.
//
// Also, we only consider the previous base to be a missed
// weak pulse if the current base is a match.
//
if (q > 1) {
insertedPulseProb =
(target.GetPreBaseDeletionQV(t-1, query.seq[q-2]) *target.GetDeletionQV(t-1)
+ query.GetInsertionQV(q-1))
* exp(alignProb[q-1][t]);
}
else {
//
// There can be no pre-base deletion tag here (could probably be an assert statement).
//
insertedPulseProb = query.GetInsertionQV(q-1) * exp(alignProb[q-1][t]);
}
}
else {
insertedPulseProb = (query.GetInsertionQV(q-1) + target.GetDeletionQV(t-1))*exp(alignProb[q-1][t]);
}
//
// An insertion in the target may be either a normal extra base
// in the target, or a deletion in the query.
// logDeletedPulseProb = target.GetInsertionQV(t-1)) + alignProb[q][t-1];
double deletedPulseProb = 0;
if (query.GetDeletionTag(q-1) != 'N') {
if (t > 1) {
deletedPulseProb = (query.GetPreBaseDeletionQV(q-1, target.seq[t-2]) * query.GetDeletionQV(q-1)
+ target.GetInsertionQV(t-1))*exp(alignProb[q][t-1]);
}
else {
// There was a dropped pulse before this position, but nothing to align it to.
deletedPulseProb = target.GetInsertionQV(t-1) * exp(alignProb[q][t-1]);
}
}
else {
deletedPulseProb = (target.GetInsertionQV(t-1) + query.GetDeletionQV(q-1)) *exp(alignProb[q][t-1]);
}
// Determine the total probability of reaching this position.
/* cout << "align prob " << q << " " << t << " " << logMatchedPulseProb << " "
<< logInsertedPulseProb << " " << logDeletedPulseProb << endl;*/
alignProb[q][t] = log(matchedPulseProb + insertedPulseProb + deletedPulseProb);
}
}
double fullAlignProb = alignProb[numRows-1][numCols-1];
alignProb.Free();
return fullAlignProb;
}
#endif // _BLASR_FULL_QV_ALIGN_HPP_
|