Блог

Я хочу создать функцию, которая считывает FASTA-файл с последовательностями ДНК (возможно, неоднозначными) и вводит подпоследовательность, которая возвращает все идентификаторы последовательностей последовательностей, содержащих данную подпоследовательность.

Чтобы сделать скрипт более эффективным, я попытался использовать nt_search , чтобы предоставить все возможности неоднозначной последовательности из FASTA. Это казалось более эффективным, чем создание всех однозначных возможностей, особенно для больших последовательностей в FASTA-файлах.

Прямо сейчас я пытаюсь понять, как я могу проверить, является ли подпоследовательность частью выходных данных, предоставляемых nt_search .

Я хочу посмотреть, является ли eg 'CGC' (входная подпоследовательность) частью возможностей, предоставляемых nt_search : ['TA[GATC][AT][GT]GCGGT'] , и вернуть все идентификаторы последовательностей последовательностей, для которых это верно.

Что у меня есть на данный момент:

 def bonus_subsequence(file, unambiguous_sequence):
    seq_records = SeqIO.parse(file,'fasta', alphabet =ambiguous_dna) 
    resultListOfSeqIds = [] 
    print(f'Unambiguous sequence {unambiguous_sequence} could be a subsequence of:')
    for record in seq_records:
        d = Seq.IUPAC.IUPACData.ambiguous_dna_values 
        couldBeSubSequence = False;    
        if unambiguous_sequence in nt_search(unambiguous_sequence,record): 
            couldBeSubSequence = True; 
        if couldBeSubSequence == True:
            print(f'{record.id}')
            resultListOfSeqIds.append({record.id})
  

На втором этапе я хочу иметь возможность также использовать это для неоднозначных подпоследовательностей, но я был бы более чем рад помощи по этому первому вопросу, заранее спасибо!

Я не знаю, хорошо ли я вас понял, но вы можете попробовать это:

Пример файла fasta:

 >seq1
ATGTACGTACGTACNNNNACTG
>seq2
NNNATCGTAGTCANNA
>seq3
NNNNATGNNN
  

Код:

 from Bio import SeqIO
from Bio import SeqUtils
from Bio.Alphabet.IUPAC import ambiguous_dna

if __name__ == '__main__':

  sub_seq = input('Enter a subsequence: ')
  results = []

  with open('test.fasta', 'r') as fh:
      for seq in SeqIO.parse(fh, 'fasta', alphabet=ambiguous_dna):
        if sub_seq in seq:
          results.append((seq.name))

  print(results, sep='n')
  

Результаты (консоль):

 Enter a subsequence: ATG
Results:
seq1    
seq3   
  

 Enter a subsequence: NNNA
Results:
seq1    
seq2    
seq3