Python код для оптимальной оценки выравнивания и последовательности, дающей неправильный результат - PullRequest
0 голосов
/ 05 февраля 2020

Это мой первый код, поэтому, пожалуйста, поймите, мой код очень грязный. Я сделал два разных способа получить оптимальный результат и оптимальную последовательность, к сожалению, оба моих ответа неверны. В моем коде я включил способ открыть файл fasta, но так как это, похоже, не сработало, я просто сам включил последовательности в код.

Мой оптимальный результат вычислен, но по какой-то причине не распечатан - это также неправильно, у меня 208, когда я должен получить 275. Я также не получаю правильный счет выравнивания назад. Две последовательности: необходимо выполнить выравнивание оценок, 11 для внутренних зазоров, 8 для концевых зазоров на 5'-конце, 7 для зазоров на 3'-конце, 4 для несовпадений, 0 для совпадений

Мой файл находится на [удаленная ссылка]

     my_file = open("one.fasta","w")

 my_file.write (""">Testseq1
 TCTGGTGTCCTAGGCGTAGAGGAACCACACCAATCCATCCCGAACTCTGGTGGTTAAACTCTACTGCGGTGACGATACT""")

 sequenceone= open("one.fasta","r")

 line = sequenceone.readline()
 header = ""
 seqA = ""
 while line:
    line = line.rstrip("\n")
    if ">" in line:
    header = line 
else :
    seqA = seqA + line
line = sequenceone.readline()
my_file.close()

my_files = open("two.fasta","w")
my_files.write (""">Testseq2
TGGTGCGGTCATACCAGCGCTAATGCACCGGATCCCATCAGAACTCCGCAGTTAAGCGCGCTTGGGCCAGAACAGTACTGGGATGGGTGTCC""")
sequencetwo= open("two.fasta","r")

line = sequencetwo.readline()
header = ""
seqB = ""
while line:
    line = line.rstrip("\n")
     if ">" in line:
        header = line 
    else :
        seqB = seqB + line
    line = sequencetwo.readline()
    my_files.close()

alphabet = ["A","C","G","T"]
score = [[8,8,8,8,8],\
         [0,4,4,4,11],\
         [4,0,4,4,11],\
         [4,4,0,4,11],\
         [4,4,4,0,11],\
         [7,7,7,7,7]]


def Global(a,b):
    D = []
    for i in range(len(a)+1):
        D.append([0]* (len(b)+1))

    for i in range(len(a)+1):
        D[i][0] = D[i-1][0] + score[alphabet.index(a[i-1])][-1]
    for i in range(len(b)+1):
        D[0][i] = D[0][i-1] + score[-1][alphabet.index(b[i-1])]

    for i in range (1, len(a)+1):
        for j in range (1, len(b)+1):
            distHor = D[i][j-1] + score[-1][alphabet.index(b[j-1])]
            distVer = D[i-1][j] + score[alphabet.index(a[i-1])][-1]
            if a[i-1] == b[j-1]:
                distDiag = D[i-1][j-1]
            else:
                distDiag = D[i-1][j-1] + score[alphabet.index(a[i-1])][alphabet.index(b[j-1])]

            D[i][j] = min(distHor, distVer, distDiag)

    return D[-1][-1]

seqA = "TCTGGTGTCCTAGGCGTAGAGGAACCACACCAATCCATCCCGAACTCTGGTGGTTAAACTCTACTGCGGTGACGATACT"
seqB = "TGGTGCGGTCATACCAGCGCTAATGCACCGGATCCCATCAGAACTCCGCAGTTAAGCGCGCTTGGGCCAGAACAGTACTGGGATGGGTGTCC"
row = len(seqA)+1
column = len(seqB)+1
match = 0
mismatch = 4
gap = 11
align1=""
align2=""
matrix=[[[[None] for i in range (2)] for i in range(column)] for i in range(row)]
for i in range(column):
    matrix[0][i][0]=gap*i
    if(i>0):
        matrix[0][i][1]="hor"
for i in range(row):
    matrix[i][0][0]=gap*i
    if(i>0):
        matrix[i][0][1]="ver"

for i in range(1,row):
    for j in range(1,column):
        hor=matrix[i][j-1][0]+gap
        ver=matrix[i-1][j][0]+gap
        if (seqA[i-1]==seqB[j-1]):
            diag=matrix[i-1][j-1][0]+match
        else:
            diag=matrix[i-1][j-1][0]+mismatch
        var = {hor:"hor",ver:"ver",diag:"diag"}
        hvd=[hor,ver,diag]
        matrix[i][j][0]=max(hvd)
        matrix[i][j][1]=var.get(max(var))

k=row
l=column
while(True):
    if(l==1 and k==1):
        break
    else:
        if(matrix[k-1][l-1][1]=="ver"):
            align1+=seqA[k-2]
            align2+="-"
            k-=1
        elif(matrix[k-1][l-1][1]=="hor"):
            align1+="-"
            align2+=seqB[l-2]
            l-=1
        elif(matrix[k-1][l-1][1]=="diag"):
            align1+=seqA[k-2]
            align2+=seqB[l-2]
            k-=1
            l-=1

align1=align1[::-1]
align2=align2[::-1]
print (align1)
print (align2)
Global(seqA,seqB)

Может кто-нибудь подсказать мне, что я делаю неправильно?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...