Demo entry 6789519

python3

   

Submitted by anonymous on Apr 30, 2019 at 04:17
Language: Python 3. Code size: 1.9 kB.

# -*- coding: utf-8 -*-
# @Author: orres
# @Email:  jasonjin22@gmail.com
# @Date:   2019-04-29 11:08:22
# @Last Modified by:   orres
# @Last Modified time: 2019-04-29 13:22:37
import copy
import math

# used to calcualte the score
def score(a, b):
	if (a == b):
		return (1)
	else:
		return (-1)


P = 'GGT'
T = 'ACCGGTGCTGTAA'
penalty = -1

global_max = 0
best_k = 1
best_seq = ''

for k in range(1, math.floor(len(T) / len(P) + 1)):
	# iterate k from 1 to maximum
	string = P
	for i in range(k-1):
		string += P
	P = string
	# init the dp matrix
	dp = []
	trace = []
	for i in range(k*len(P) + 1):
		line = []
		for j in range(len(T) + 1):
			line.append(0)
		dp.append(line)
		trace.append(copy.deepcopy(line))
	for i in range(k*len(P) + 1):
		dp[i][0] = -i

	# calculate the remaining elements in matrix
	for i in range(1, len(P) + 1):
		for j in range(1, len(T) + 1):
			dp[i][j] = max(dp[i-1][j-1] + score(P[i-1], T[j-1]), dp[i-1][j] + penalty, dp[i][j-1]+penalty)
			# the trace matrix used to record when the element comes from
			if (dp[i][j] == dp[i-1][j] + penalty):
				trace[i][j] = 2
			elif (dp[i][j] == dp[i-1][j-1] + score(P[i-1], T[j-1])):
				trace[i][j] = 1
			else:
				trace[i][j] = 3
	# find the highest score in the matrix
	big = 0
	pos = (0, 0)
	for i in range(k*len(P) + 1):
		for j in range(len(T) + 1):
			if dp[i][j] >= big:
				big = dp[i][j]
				pos = (i, j)

	# get the aligned sequence
	x = pos[0]
	y = pos[1]
	seq = ""
	for i in range(len(P)):
		if (x < 0 or y < 0):
			break
		if (trace[x][y] == 1):
			seq = T[y-1] + seq
			x -= 1
			y -= 1
		elif (trace[x][y] == 2):
			x -= 1
			seq = '-' + seq
		else:
			y -= 1

	# choose the best k
	if (big > global_max):
		global_max = big
		best_k = k
		best_seq = seq

# print the best k and the best aligned sequence
print(best_k)
print(best_seq)

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).