# Pairwise Alignment Score function: From Python to Julia. Suggestions

**URL:** <https://discourse.julialang.org/t/pairwise-alignment-score-function-from-python-to-julia-suggestions/48089>\
**Category:** Biology, Health, and Medicine\
**Tags:** question, suggestions\
**Created:** [October 10, 2020, 6:07am UTC](https://discourse.julialang.org/t/pairwise-alignment-score-function-from-python-to-julia-suggestions/48089 "2020-10-10T06:07:54Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![aadam](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aadam/32/7113_2.png) [@aadam](https://discourse.julialang.org/u/aadam)\
**Post date:** [October 10, 2020, 6:07am UTC](https://discourse.julialang.org/t/pairwise-alignment-score-function-from-python-to-julia-suggestions/48089/1 "2020-10-10T06:07:54Z")

</div>

I wanted a function to compute the Pairwise Alignment Score in Julia, but didn’t find any in the BioJulia repository. I submitted an [issue](https://github.com/BioJulia/BioAlignments.jl/issues/43), asking for some guidance, but no response till now. So, after waiting a long time, I resorted to porting my Python implementation into Julia. Here’s what I got:

## Python Implementation

```python
def score_pairwise(self, seq1, seq2, matrix, gap=True):
        """Pairwise score generator.

        Keyword arguments:
        seq1 -- first alignment
        seq2 -- second alignment
        matrix -- substitution matrix"""
        for A, B in zip(seq1, seq2):
            diag = ('-' == A) or ('-' == B)
            yield (self.gap_extend_score if gap else self.gap_open_score) if diag else matrix[(A, B)]
            gap = diag

```

## Julia port

```julia
function score(S1, S2; matrix=BLOSUM90, gop=-1, gep=-1)
	score = 0
	prev_gap = false
	for (A, B) in zip(S1, S2)
		diag = ('-' == A) || ('-' == B)
		score += (diag) ? (prev_gap ? gep : gop) : matrix[A, B]
		prev_gap = diag
	end
	score
end

```

**Disclaimer** : I’m pretty new to Julia and I feel this isn’t the best implementation possible and we could improve this a lot. **Any suggestions would be welcome** 🙂

## Problem

The result I’m getting from this method is a little bit off from the one I get from the BioJulia package.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/0/c02f191140d4fc942fb914ac3b07a7dbfd31fc73.png)  
I’ve tried it on a number of sequences, and every time it’s a little bit off, sometimes more than BioJulia score and sometimes less.

Can anyone help me in figuring out the problem? I would be very thankful.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [October 10, 2020, 7:40am UTC](https://discourse.julialang.org/t/pairwise-alignment-score-function-from-python-to-julia-suggestions/48089/2 "2020-10-10T07:40:43Z")

</div>

Hey and welcome to the community 👋  
In general, it’s a good idea to include some data-generating code so that people can copy-paste your code into a repl and have it run. That way it’s much easier to help you.

---

<div class="post-metadata">

**Author:** ![aadam](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aadam/32/7113_2.png) [@aadam](https://discourse.julialang.org/u/aadam)\
**Post date:** [October 10, 2020, 9:32am UTC](https://discourse.julialang.org/t/pairwise-alignment-score-function-from-python-to-julia-suggestions/48089/3 "2020-10-10T09:32:47Z")

</div>

Sure thing.

### Imports required

```julia
using BioSequences
using BioAlignments

```

### Score using BioJulia

```julia-auto
function score_bio(S1, S2; matrix=BLOSUM90, gop=-3, gep=-1, asterisk_penalty=missing)
	scoremodel = AffineGapScoreModel(matrix, gap_open=gop, gap_extend=gep)
	res = pairalign(GlobalAlignment(), S1, S2, scoremodel)
	res
end

score_bio("ACCTAACTAAC", "AACTTTACAATAC")

```

This returns:

```julia-auto
BioAlignments.PairwiseAlignmentResult{Int64,String,String}:
  score: 44
  seq: 1 ACCTA-ACTA-AC 11
          | || || | ||
  ref: 1 AACTTTACAATAC 13

```

The problem with using this method is that if my sequences contain `'-'` then it won’t run and results in an exception.  
I am getting the alignments from another manual method and here, I just want to calculate the score of the generated alignments. Because the generated alignments have `'-'` in them, I can’t use this method to get the score, even using `scoreonly=true` parameter doesn’t help.

So I coded my own scoring function, which is:

```julia-auto
function pairwise_score(S1, S2; matrix=BLOSUM90, gop=-3, gep=-1)
	score = 0
	prev_gap = false
	for (A, B) in zip(S1, S2)
		diag = ('-' == A) || ('-' == B)
		score += (diag) ? (prev_gap ? gep : gop) : matrix[A, B]
		prev_gap = diag
	end
	score
end

```

To test its accuracy, I passed it the alignments that were generated from the above function, `score_bio`:

```julia-auto
pairwise_score("ACCTA-ACTA-AC", "AACTTTACAATAC")

```

This gave the score:

```julia-auto
46

```

while the `score_bio` function was returning a score of `44` for the same alignment.

I tried different input sequences, but the result is always the same, i.e. my output is off by a couple of numbers.

Hopefully, you’ll be able to reproduce the problem now. 🙂  
I’d also appreciate suggestions on how you think the Julia implementation of the above python code could be improved further (made succint, or more Julian-way 🙄). Thanks 🙂

---

<div class="post-metadata">

**Author:** ![jonathanBieler](https://avatars.discourse-cdn.com/v4/letter/j/82dd89/32.png) [@jonathanBieler](https://discourse.julialang.org/u/jonathanBieler)\
**Post date:** [October 10, 2020, 1:23pm UTC](https://discourse.julialang.org/t/pairwise-alignment-score-function-from-python-to-julia-suggestions/48089/4 "2020-10-10T13:23:24Z")

</div>

I think your code is correct, and that the difference is that BioAlignments adds up the opening and extension gap penalties when you have a new gap (here you have two gaps so you end up with -2 of difference) :

[https://github.com/BioJulia/BioAlignments.jl/blob/8f1bdca96013b855665458c1d2ddbd9614d46ba6/src/pairwise/algorithms/needleman\_wunsch.jl#L124](https://github.com/BioJulia/BioAlignments.jl/blob/8f1bdca96013b855665458c1d2ddbd9614d46ba6/src/pairwise/algorithms/needleman_wunsch.jl#L124)

I don’t know if that’s a bug in BioAlignments or if that’s an alternative definition. But to me the opening and extension penalties sounds like they should be exclusive. If I’m correct running your score with `gop=-4` should give you the same score as BioAlignments.

Edit : not a bug, the docs explain the scoring model here :

[https://biojulia.net/BioAlignments.jl/latest/pairalign.html#Alignment-types-and-scoring-models-1](https://biojulia.net/BioAlignments.jl/latest/pairalign.html#Alignment-types-and-scoring-models-1)

You can either adopt their definition or rescale you gap opening penalty.