Ошибка при использовании smith-waterman из skbio 0.5.4

Question

Я использую упакованную версию smith-waterman из skbio (0.5.4), но у меня есть неожиданная ошибка:

_, оценка, _ = local_pairwise_align_ssw(белок-список [idx1], белок-список [idx2], substitution_matrix = blosum62)
Файл "/anaconda3/lib/python3.6/site-packages/skbio/alignment/_pairwise.py", строка 732, в local_pairwise_align_ssw
проверить =False)
Файл "/anaconda3/lib/python3.6/site-packages/skbio/alignment /_tabular_msa.py", строка 785, в __init__
reset_index = minter - Нет, а index - Нет.)
Файл "/anaconda3/lib/python3.6/site-packages/skbio/alignment /_tabular_msa.py", строка 1956, в расширении
self._assert_valid_sequences(последовательность)
Файл "/anaconda3/lib/python3.6/site-packages/skbio/alignment /_tabular_msa.py", строка 2035, в _assert_valid_sequence
% (длина, ожидаемая длина)
ValueError: Длина каждой последовательности должна соответствовать количеству позиций в MSA: 232!= 231

Странно то, что иногда ошибка появляется с белковой парой 0-10, а другие с 0-116. Так что я не верю, что это ошибка белка.

Answer 1

У меня аналогичная проблема. Однако мне удалось ограничить ошибку оптимизированной версией SSW. Так что ошибок в форматировании последовательности нет.

import warnings
from skbio.sequence import Protein
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", message="...")
    from Bio.Align import substitution_matrices
from skbio.alignment import local_pairwise_align_ssw
from skbio.alignment import local_pairwise_align

peptide1 = Protein("CGAGDNQAGTALIF")
peptide2 = Protein("CAGEEGGGADGLTF")
gap_open_penalty = 10
gap_extend_penalty = 10
substitution_matrix = substitution_matrices.load("BLOSUM45")

## works correct
rv = local_pairwise_align_ssw(
      sequence1 = peptide1
    , sequence2 = peptide2
    , gap_open_penalty=1
    , gap_extend_penalty=1
    , substitution_matrix=substitution_matrix
)
print(rv)

## but if I swap peptide1 and peptide 2 the ValueError occur
rv = local_pairwise_align_ssw(
      sequence1 = peptide2
    , sequence2 = peptide1
    , gap_open_penalty=1
    , gap_extend_penalty=1
    , substitution_matrix=substitution_matrix
)
print(rv)

## if I do the same with local_pairwise_align it works!
rv = local_pairwise_align(
      seq1=peptide2
    , seq2=peptide1
    , gap_open_penalty=1
    , gap_extend_penalty=1
    , substitution_matrix=substitution_matrix
)
print(rv)