Skip to content

Commit

Permalink
Merge pull request #8 from ppillot/gap-penalty-fixes
Browse files Browse the repository at this point in the history 
Gap-penalty-fixes
  • Loading branch information
@ppillot
ppillot committed Jan 7, 2024
2 parents 9126374 + 654e271 commit 157bb91
Show file tree
Hide file tree
Showing 4 changed files with 35 additions and 41 deletions.
25 changes: 12 additions & 13 deletions src/align/align.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,14 +73,14 @@ export function pairwiseAlignment (
isOdd = 0;

// Traceback matrix size: for convenience, an extra column has been added
// but will remain half empty (TODO: fix iteration boundaries).
// but will remain half empty.
// The traceback values are stored in a typed array of (n+1) * m / 2 bytes.
// Using bit values for M, I and D transitions, shifted depending on the matrix
// they encode, it's possible to encode all states using 4 bits.
// As a cell in the array contains 8 bits, each ultimately contains the
// results for 2 successive values in the DP matrix.

const tbM = new Uint8Array(Math.ceil((lSeqALen + 1) * lSeqBLen / 2)); // Trace back matrix
const tbM = new Uint8Array(Math.ceil((lSeqALen + 1) * lSeqBLen / 2) + 1); // Trace back matrix


let gapOpenA = 0,
Expand Down Expand Up @@ -221,9 +221,7 @@ export function pairwiseAlignment (

}

// if (DEBUG) Log.add('End DP computation');

var score = Math.max(lMatch, lLastInsert, lDelArr[lDelArr.length - 1]);
const score = Math.max(lMatch, lLastInsert, lDelArr[lDelArr.length - 1]);

// Traceback
const [lEpathA, lEpathB] = tracebBackToEpaths(tbM, lSeqALen, lSeqBLen, tb);
Expand Down Expand Up @@ -257,7 +255,7 @@ export function MSASeqAlignment(
let lMatch = 0.0, // match score
lMatchArr: number[] = [],
lDelArr : number[] = [];
const tbM = new Uint8Array(Math.ceil((lProfALen + 1) * (lSeqBLen + 1) / 2));
const tbM = new Uint8Array(Math.ceil((lProfALen + 1) * (lSeqBLen + 1) / 2) + 1);
let tbIdx: number;
let isOdd: number;

Expand All @@ -268,6 +266,7 @@ export function MSASeqAlignment(
tb = 0,
lLastInsert = 0.0,
lPrevLastInsert = 0.0,
lPrevDel = 0.,
lPrevMatch = 0.0,
lDeletej_1 = 0.0,
lInserti_1 = 0.0,
Expand All @@ -276,7 +275,7 @@ export function MSASeqAlignment(
/** Position specific gap close penalty of profile A at pos i-1 */
lProfAGapCP = 0.0,
/** Position specific substitution scores of profile A at pos i-1 */
lProfASScores = new Float32Array(params.abSize);
lProfASScores: Float32Array;

const GAP_OPEN = params.gapOP * wB;
const GAP_OPEN_B = GAP_OPEN / 2;
Expand Down Expand Up @@ -328,7 +327,7 @@ export function MSASeqAlignment(

//Delete i,j score computation
lGapOpenA = lMatchArr[j] + lProfAGapOP; //
lGapExtendA = lDelArr[j];
lGapExtendA = lPrevDel = lDelArr[j];
if (j === lSeqBLen) { //terminal penalties are halved
lGapOpenA -= lProfAGapOP / 2;
}
Expand All @@ -355,7 +354,7 @@ export function MSASeqAlignment(
}

//Match i,j score computation
lDeletej_1 = lDelArr[j] + lProfAGapCP; //it should be prev
lDeletej_1 = lPrevDel + lProfAGapCP; //it should be prev
lInserti_1 = lPrevLastInsert + GAP_CLOSE_B;
lMatch = lMatchArr[j - 1] + lProfASScores[sB[j - 1]];

Expand Down Expand Up @@ -421,7 +420,7 @@ export function MSAMSAAlignment(

const lProfALen = nodeA.profile.length,
lProfBLen = nodeB.profile.length,
tbM = new Uint8Array(Math.ceil((lProfALen + 1) * (lProfBLen + 1) / 2));
tbM = new Uint8Array(Math.ceil((lProfALen + 1) * (lProfBLen + 1) / 2) + 1);
let tbIdx: number;
let isOdd: number;

Expand All @@ -436,6 +435,7 @@ export function MSAMSAAlignment(
lLastInsert = 0.0,
lPrevLastInsert = 0,
lPrevMatch = 0.0,
lPrevDel = 0.,
lDeletej_1 = 0.0,
lInserti_1 = 0.0,
lProfAGapOP = 0.0,
Expand Down Expand Up @@ -485,7 +485,7 @@ export function MSAMSAAlignment(

//Delete i,j score computation
lGapOpenA = lMatchArr[j] + lProfAGapOP; //
lGapExtendA = lDelArr[j];
lGapExtendA = lPrevDel = lDelArr[j];
if (j === lProfBLen && !(opt & ALIGNOPT.DISABLE_FAVOR_END_GAP)) {
lGapOpenA -= lProfAGapOP / 2;
}
Expand Down Expand Up @@ -524,8 +524,7 @@ export function MSAMSAAlignment(
k++;
}

//match = Match[j - 1] + _utils.sumOfPairsScorePP3(profAAAScores, profB[j - 1]);
lDeletej_1 = lDelArr[j] + lProfAGapCP;
lDeletej_1 = lPrevDel + lProfAGapCP;
lInserti_1 = lPrevLastInsert + profB.m_ScoreGapClose[j - 1];

if (j === 1 && !(opt & ALIGNOPT.DISABLE_FAVOR_END_GAP)) { //terminal penalties are halved
Expand Down
8 changes: 2 additions & 6 deletions src/sequence/profile.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,11 +322,9 @@ export function seqToProf (pSeq: TSequence, pWeight: number, params: TAlignmentP
// Corrections for end/begin gap penalties
if (!(opt & ALIGNOPT.DISABLE_FAVOR_START_GAP)) {
lProf.m_ScoreGapOpen[0] /= 2;
lProf.m_ScoreGapOpen[l-1] /= 2;
}
if (!(opt & ALIGNOPT.DISABLE_FAVOR_END_GAP)) {
lProf.m_ScoreGapClose[0] /= 2;
lProf.m_ScoreGapClose[l-1] /= 2;
lProf.m_ScoreGapOpen[l-1] /= 2;
}


Expand Down Expand Up @@ -466,11 +464,9 @@ export function mergeProfiles (pProfA: ProfPos, pProfB: ProfPos, pESA: number[],
// Corrections to end/begin gap penalties
if (!(opt & ALIGNOPT.DISABLE_FAVOR_START_GAP)) {
lProf.m_ScoreGapOpen[0] /= 2;
lProf.m_ScoreGapOpen[l-1] /= 2;
}
if (!(opt & ALIGNOPT.DISABLE_FAVOR_END_GAP)) {
lProf.m_ScoreGapClose[0] /= 2;
lProf.m_ScoreGapClose[l-1] /= 2;
lProf.m_ScoreGapOpen[l-1] /= 2;
}


Expand Down
31 changes: 15 additions & 16 deletions src/sequence/sequence.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,28 +387,27 @@ function fractionalIdentity(seqA: string, seqB: string) {
seqA = seqA || '';
seqB = seqB || '';

var match = 0,
let match = 0,
tokenA = '',
tokenB = '',
lA = 0,
lB = 0,
maxI = seqA.length;
lB = 0;
const maxI = seqA.length;

for (var i = 0; i < maxI; i++) {
for (let i = 0; i < maxI; i++) {
tokenA = seqA.charAt(i);
tokenB = seqB.charAt(i);
lA += (tokenA !== '-') ? 1 : 0;
lB += (tokenB !== '-') ? 1 : 0;

if ((tokenA === '-') || (tokenB === '-')) {
lA += (tokenA !== '-') ? 1 : 0;
lB += (tokenB !== '-') ? 1 : 0;
continue;
} else {
lA++;
lB++;
if (tokenA === tokenB) {
match++;
}
if (tokenA === '-' || tokenB === '-') {
continue
}

if (tokenA === tokenB) {
match++;
}

}
return match / Math.min(lA, lB);
}
Expand All @@ -428,15 +427,15 @@ export function distanceKimura (msa: string[]) {
let dk = 0; // Computed distance

//parcours de la matrice
for (var i = 0; i < l; i++) {
for (let i = 0; i < l; i++) {

if ( distMatrix[i] === undefined ) {
distMatrix[i] = new Array(l);
}

distMatrix[i][i] = 0;

for (var j = i + 1; j < l; j++) {
for (let j = i + 1; j < l; j++) {

if (distMatrix[j] === undefined) {
distMatrix[j] = new Array(l);
Expand Down
12 changes: 6 additions & 6 deletions src/sequence/tree.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ export function makeTree(mD: number[][], tSeq: TSequence[]) {
let nbSeq = mD.length;
let tabIdx: number[] = [];

for (var i = 0; i < nbSeq; i++) {
for (let i = 0; i < nbSeq; i++) {
clusters[i] = {
type: NODE_TYPE.LEAF,
seq: tSeq[i],
Expand All @@ -76,7 +76,7 @@ export function makeTree(mD: number[][], tSeq: TSequence[]) {
let lNode: InternalNode;
let lMinX: number; // index in matrix of column with the min value
let lMinY: number; // index in matrix of row with the min value.
for (var i = 0, nbIterations = nbSeq - 1; i < nbIterations; i++) {
for (let i = 0, nbIterations = nbSeq - 1; i < nbIterations; i++) {
// Make a node from the min value in distance matrix

[lNode, lMinX, lMinY] = findMinInDistanceMatrix(mD, tabIdx);
Expand Down Expand Up @@ -166,7 +166,7 @@ function recomputeDistMatrix(matrix: number[][], x: number, y: number, tI: numbe
const l = matrix.length;

// loop matrix rows
for (var i = 0; i < l; i++) {
for (let i = 0; i < l; i++) {
if (i == x || i == y) continue;

// avg computation from MAFFT
Expand Down Expand Up @@ -259,7 +259,7 @@ function recomputeDistMatrix(matrix: number[][], x: number, y: number, tI: numbe
const root = cluster[cluster.length - 1];

function setWeight (node: LeafNode|InternalNode) {
var edgeLength = 0,
let edgeLength = 0,
nbLeaves = 0;

switch (node.type) {
Expand Down Expand Up @@ -301,8 +301,8 @@ function recomputeDistMatrix(matrix: number[][], x: number, y: number, tI: numbe
* treeB.
*/
export function compareTrees (treeA: Tree, treeB: Tree) {
var i = 0,
nbLeaves = (treeA.length + 1) / 2 + 1, //chaque noeud possède deux descendants (peuvent être un noeud ou une feuille). Le nombre de noeud est égal au nombre de feuilles -1
let i = 0;
const nbLeaves = (treeA.length + 1) / 2 + 1, //chaque noeud possède deux descendants (peuvent être un noeud ou une feuille). Le nombre de noeud est égal au nombre de feuilles -1
//nbNodes = nbLeaves - 1,

compareParentNodes = function (nodeA: InternalNode|LeafNode, nodeB: InternalNode|LeafNode) {
Expand Down

0 comments on commit 157bb91

Please sign in to comment.