如何过滤“生产性”氨基酸序列

Question

我有一个具有不同氨基酸序列的fasta文件，例如：

>abc
HSTSDSAQTMFPVALLLLAAGSCVKGEQLTQPTSVTVQPGQRLTITCQVSYSLGTYFTAW
IRQPAGKGLEWIGMRSTGASYYKDSLKNKFSIDLDTSSKTVTLNGQNVQPEDTAVYYCAR
APSRGFDYWGKGTMVTITSATPKGPTVFPL

>def
TARQIQHKPCFL*LCCCWQLDHV*RVNS*HSRPL*LCSQVNV*PSPVRSLILLVPTSQLG
SDSLQEKDWSGLE*DLLELHTTKIH*RTSSVST*TLPAKL*L*MDRMCSLKTLLCITVPE
RPVGVLTTGGKAPWSPSPRPPQRDQLCFL*

>ghi
GSQHVRFSTNHVSCSSAAVGSWIMCEG*TVDTADLCDCAARSTSDHHLSGLLFSW*LLHS
LDQTACRKRTGVDWEQIYWSCILQRFIKEQVQYRLRHFQQNCDSKWTECAA*RHCCVLLC
QTTGSGSWLLGERHHGHHHLGHPKGTNCVSS

我想从“非生产性”序列中过滤掉“生产性”序列。

附加信息：我已将所有 6 帧中的每个 DNA 序列翻译为氨基酸序列。

“非生产性”是指那些不翻译成蛋白质的蛋白质（没有氨基酸 M 和/或有太多终止密码子）。我想在 fasta 文件中过滤掉这些非生产性序列。

至于“高效”的，我还想将每个“高效”序列仅与完整帧一起保存在另一个 fasta 文件中。

我正在尝试用 python 来做，但我被困住了，所以欢迎任何想法。

提前谢谢您

Answer 1

使用

biopython

和终止密码子数量阈值的示例。

# pip install biopython
from Bio import SeqIO

seqs = SeqIO.parse(open('example.fasta'), 'fasta')
productive = {}
for s in seqs:
    productive.setdefault(s.count('*')<3, []).append(s)

print(productive)

输出：

{True:  [SeqRecord(seq=Seq('HSTSDSAQTMFPVALLLLAAGSCVKGEQLTQPTSVTVQPGQRLTITCQVSYSLG...FPL'), id='abc', name='abc', description='abc', dbxrefs=[])],
 False: [SeqRecord(seq=Seq('TARQIQHKPCFL*LCCCWQLDHV*RVNS*HSRPL*LCSQVNV*PSPVRSLILLV...FL*'), id='def', name='def', description='def', dbxrefs=[]),
         SeqRecord(seq=Seq('GSQHVRFSTNHVSCSSAAVGSWIMCEG*TVDTADLCDCAARSTSDHHLSGLLFS...VSS'), id='ghi', name='ghi', description='ghi', dbxrefs=[])]}

您可以通过将

s.count('*')<3

替换为自定义函数来添加模式复杂逻辑：

from Bio import SeqIO

seqs = SeqIO.parse(open('test/stackoverflow/example.fasta'), 'fasta')

def is_productive(s):
    # does the sequence start with M and contain less than 3 stops?
    return s.seq.startswith('M') and (s.count('*')<3)

productive = {}
for s in seqs:
    productive.setdefault(is_productive(s), []).append(s)

写为fasta：

with open('productive_seqs.fasta', 'w') as fw:
    SeqIO.write(productive.get(True, []), fw, 'fasta')
with open('nonproductive_seqs.fasta', 'w') as fw:
    SeqIO.write(productive.get(False, []), fw, 'fasta')

输出：

# productive_seqs.fasta
>abc
HSTSDSAQTMFPVALLLLAAGSCVKGEQLTQPTSVTVQPGQRLTITCQVSYSLGTYFTAW
IRQPAGKGLEWIGMRSTGASYYKDSLKNKFSIDLDTSSKTVTLNGQNVQPEDTAVYYCAR
APSRGFDYWGKGTMVTITSATPKGPTVFPL

# nonproductive_seqs.fasta
>def
TARQIQHKPCFL*LCCCWQLDHV*RVNS*HSRPL*LCSQVNV*PSPVRSLILLVPTSQLG
SDSLQEKDWSGLE*DLLELHTTKIH*RTSSVST*TLPAKL*L*MDRMCSLKTLLCITVPE
RPVGVLTTGGKAPWSPSPRPPQRDQLCFL*
>ghi
GSQHVRFSTNHVSCSSAAVGSWIMCEG*TVDTADLCDCAARSTSDHHLSGLLFSW*LLHS
LDQTACRKRTGVDWEQIYWSCILQRFIKEQVQYRLRHFQQNCDSKWTECAA*RHCCVLLC
QTTGSGSWLLGERHHGHHHLGHPKGTNCVSS

注意，如果只需要文件，可以直接在循环中写入序列：

from Bio import SeqIO

seqs = SeqIO.parse(open('test/stackoverflow/example.fasta'), 'fasta')

def is_productive(s):
    return s.seq.startswith('M') and (s.count('*')<3)

with open('productive_seqs.fasta', 'w') as fw1, open('nonproductive_seqs.fasta', 'w') as fw2:
    for s in seqs:
        SeqIO.write(s, fw1 if is_productive(s) else fw2, 'fasta')