从FASTA文件中找到DNA序列中DNA子序列的序列ID

问题描述

我想要一个函数来读取具有DNA序列（可能是模棱两可）的FASTA文件，并输入一个子序列，该子序列返回包含给定子序列的序列的所有序列ID。

为了使脚本更有效，我尝试使用nt_search来使FASTA产生歧义序列的所有可能性成为可能。这似乎比产生所有明确的可能性更有效，尤其是对于较大序列的FASTA文件。

现在，我正在努力查看如何检查子序列是否为nt_search给出的输出的一部分。

我想看看'CGC'（输入子序列）是否是nt_search：['TA[GATC][AT][GT]GCGGT']给出的可能性的一部分，并返回其为真的序列的所有序列ID。 / p>

到目前为止我所拥有的：

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