如何过滤“生产性”氨基酸序列 [已关闭]

How to filter "productive" amino acid sequences [closed]

提问人:Silvia 提问时间:9/29/2023 最后编辑:pippo1980Silvia 更新时间:10/24/2023 访问量:80

问:


想改进这个问题吗?通过编辑这篇文章来更新问题,使其仅关注一个问题。

上个月关闭。

我有一个具有不同氨基酸序列的 fasta 文件,用于

例::example.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 R 序列 FASTA

评论

0赞 PGSA 9/29/2023
这很可能会被关闭,因为要求推荐工具是违反网站规则的——它往往会导致基于意见的答案和论点。我想说最重要的问题之一是数据将由谁使用?如果您可以使用多个软件包,那么选择与该过程的下一步最兼容的软件包可能是有意义的。
0赞 Silvia 9/29/2023
你好!我是新来的,我是生物信息学的应届毕业生,所以编程目前不是我的强项......我问是否有任何工具可以做到这一点,因为它对我来说是目前最快和最简单的......至于编程,我知道可以使用一些循环等来完成,但我仍在努力编程,所以欢迎任何想法......
0赞 pippo1980 10/24/2023
biology.stackexchange.com/questions/56939/......所有的蛋白质都是从蛋氨酸开始的吗?

答:

3赞 mozway 9/29/2023 #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) -> bool:
    # 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')

评论

0赞 Silvia 9/29/2023
我想先过滤掉“非生产性”的,它们不包含任何“M”。然后在“生产性”密码子中,过滤掉中间有一个 M 或最多 2 个且不超过三个终止密码子的密码子(开头和结尾限制为 10 个氨基酸)。非常感谢您的帮助!
0赞 mozway 9/29/2023
@Silvia您可以轻松地使用嵌套条件,然后:if s.seq.startswith('M'): if (s.count('*')<3): # then file1 ; else: # file2 ; else: # file3
0赞 Silvia 9/29/2023
我一直在检查,许多“生产性”序列被过滤掉为“非生产性”。问题在于,M不需要在序列的开头,它可以是大约10个氨基酸......我怎么能这么说呢?我尝试了 s.seq.find('M'),稍微好一点,但还不完全是......谢谢!
0赞 mozway 9/29/2023
对于任何位置,对于前 10 个氨基酸。对于第一个之前:'M' in s.seq'M' in s.seq[:10]M*'M' in s.split('*')[0]
1赞 Silvia 9/29/2023
嗨,我试过了,效果很好。非常感谢您的帮助!'M' in s.seq[:] and (s.count('*')<2)