Выравнивание последовательностей ДНК внутри python

Question

У меня есть тысячи последовательностей ДНК в диапазоне от 100 до 5000 bp, и мне нужно выровнять и рассчитать показатель идентификации для указанных пар. Biopython pairwise2 отлично справляется, но только для коротких последовательностей, и когда размер последовательности становится больше 2 КБ, он показывает сильную утечку памяти, которая приводит к «MemoryError», даже если используются опции «score_only» и «one_alignment_only»!

whole_coding_scores={} 
from Bio import pairwise2 
for genes in whole_coding: # whole coding is a <25Mb dict providing DNA sequences 
    alignment=pairwise2.align.globalxx(whole_coding[genes][3],whole_coding[genes][4],score_only=True,one_alignment_only=True) 
    whole_coding_scores[genes]=alignment/min(len(whole_coding[genes][3]),len(whole_coding[genes][4])) 

Результат вернулся из суперкомпьютеров:

Max vmem   = 256.114G #Memory usage of the script 
failed assumedly after job because: 
job 4945543.1 died through signal XCPU (24) 

Я знаю, что есть и другие инструменты для выравнивания, но в основном они могут просто написать счет в выходной файл, который нужно читать и проанализирован снова для получения и используя оценки выравнивания. Есть ли какой-либо инструмент, который может выровнять последовательности и возвращать оценку выравнивания внутри среды python, как это делает parwise2, но без утечки памяти?

Answer 1

Для глобального выравнивания можно попробовать NWalign https://pypi.python.org/pypi/nwalign/. Я не использовал его, но, похоже, вы можете восстановить оценку выравнивания в своем скрипте.

В противном случае, возможно, рельефность и инструменты могут помочь: http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needleall.html

Answer 2

Первый, я использовал иглу BioPython для этого. Хороший HOWTO (игнорировать унаследованный дизайн :-)) можно найти here

Второго: может быть, вы можете избежать всего набора в память с помощью генератора? Я не знаю, откуда идет ваш объект whole_coding. Но если это файл, убедитесь, что вы не читаете весь файл, а затем перебираете объект памяти. Например:

whole_coding = open('big_file', 'rt').readlines() # Will consume memory 

но

for gene in open('big_file', 'rt'):  # will not read the whole thing into memory first 
    process(gene) 

Если обработка потребности, вы могли бы написать генератор Funtion:

def gene_yielder(filename): 
    for line in open('filename', 'rt'): 
     line.strip() # Here you preprocess your data 
     yield line  # This will return 

затем

for gene in gene_yielder('big_file'): 
    process_gene(gene) 

В принципе, вы хотите, чтобы ваши программа действовать как труба: t потоки проходят через него и обрабатываются. Не используйте его в качестве кухонного горшка при приготовлении бульона: добавляйте все и применяйте тепло. Надеюсь, это сравнение не надолго :-)

Answer 3

Biopython может (сейчас). Модуль pairwise2 в Biopython ver. 1,68 (намного) быстрее и может занимать более длинные последовательности. Вот сравнение нового и старого парного2 (на 32-битном Python 2.7.11 с 2 МБ памяти, 64 бит Win7, Intel Core i5, 2.8 ГГц):

from Bio import pairwise2 

length_of_sequence = ... 
seq1 = '...'[:length_of_sequence] # Two very long DNA sequences, e.g. human and green monkey 
seq2 = '...'[:length_of_sequence] # dystrophin mRNAs (NM_000109 and XM_007991383, >13 kBp) 
aln = pairwise2.align.globalms(seq1, seq2, 5, -4, -10, -1) 

1. Старый pairwise2

   • максимальная длина/время: ~ 1900 символов/10 сек

2. Новый pairwise2

   • максимальная длина/время: ~ 7450 символов/12 сек.
   • время для 1900 символов: 1 сек

С score_only набора к True, новый pairwise2 может делать две последовательности ~ 8400 символов в 6 сек.