查找multifasta文件中的所有模式，包括重叠的主题
我有一个multifasta文件，它看起来像这样：
>NP_001002156.1
MKTAVDRRKLDLLYSRYKDPQDENKIGVDGIQQFCDDLMLDPASVSVLIVAWKFRAATQCEFSRQEFLDG
MTDLGCDSPEKLKSLLPRLEQELKDSGKFRDFYRFTFSFAKSPGQKCLDLEMAVAYWNLILSGRFKFLGL
WNTFLLEHHKKSIPKDTWNLLLDFGNMIADDMSNYAEEGAWPVLIDDFVEFARPIVTAENLQTL
>NP_957070.2
MAKDAGLKETNGEIKLFINQSPGKAAGVLQLLTVHPASITTVKQILPKTLTVTGAHVLPHMVVSTPQRPT
IPVLLTSPHTPTAQTQQESSPWSSGHCRRADKSGKGLRHFSMKVCEKVQKKVVTSYNEVADELVQEFSSA
DHSSISPNDAVSSCHVYDQKNIRRRVYDALNVLMAMNIISKDKKEIKWIGFPTNSAQECEDLKAERQRRQ
ERIKQKQSQLQELIVQQIAFKNLVQRNREVEQQSKRSPSANTIIQLPFIIINTSKKTIIDCSISNDKFEY
LFNFDSMFEIHDDVEVLKRLGLALGLESGRCSAEQMKIATSLVSKALQPYVTEMAQGSVNQPMDFSHVAA
ERRASSSTSSRVETPTSLMEEDEEDEEEDYEEEDD
>NP_123456.1
MALLLLLGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
...
虽然有一个很好的python脚本来处理multifasta文件中的主题搜索（https://www.biostars.org/p/14305/），如果使用模式“[KHR] {3}”，它将只返回主题列表和许多空结果：
>NP_001002156.1
:['RRK']
>NP_001002156.1
:[]
>NP_001002156.1
:['HHK']
>NP_957070.2
:[]
>NP_957070.2
:['RRR']
...
并且一些图案（HKK）以相同的顺序泄露。
在这里，我发现了另一个python脚本：
#coding:utf-8
import re
pattern = "[KHR]{3}"
with open('seq.fasta') as fh:
fh.readline()
seq = ""
for line in fh:
seq += line.strip()
rgx = re.compile(pattern)
result = rgx.search(seq)
patternfound = result.group()
span = result.span()
leftpos = span[0]-10
if leftpos < 0:
leftpos = 0
print(seq[leftpos:span[0]].lower() + patternfound + seq[span[1]:span[1]+10].lower())
它返回在上下文中找到的第一个匹配的基序（在匹配的基序后向前10个氨基酸，在匹配的基序之前向后10个）仅对于一个fasta（第一个）序列，对于第一个fasta序列NP_001002156.1使用scirpt ，返回的结果：
mktavdRRKldllysrykd
但它没有文件头“> NP_001002156.1”，上下文中的其他2个主题都被省略了：
glwntfllehHHKksipkdtwnl
lwntfllehhHKKsipkdtwnll
在这里，我希望所需的脚本在multifasta文件中的每个fasta序列的上下文中返回匹配的motif及其postition，并且它将显示如下结果：
>NP_001002156.1_matchnumber_1_(7~9)
mktavdrRRKldllysrykd
>NP_001002156.1_matchnumber_2_(148~150)
glwntfllehHHKksipkdtwnl
>NP_001002156.1_matchnumber_3_(149~151)
lwntfllehhHKKsipkdtwnll
>NP_957070.2_matchnumber_1_(163~165)
chvydqknirRRRvydalnvlma
>NP_123456.1
no match found
注意：匹配模式的位置不是上下文的位置。
有人可以帮帮我吗？提前致谢。
回答 1
Stack Overflow用户
发布于 2019-04-15 17:06:08
这里的“主题”是[HKR]字符的任意三长组合; 图案可能重叠。
通过在正则表达式中使用“前瞻”来解决重叠。详情见下文。引用或显示的资源似乎都没有处理，我也看不出它们如何捕获重叠的图案。
use warnings;
use strict;
use feature 'say';
my $file = shift || die "Usage: $0 fasta-file\n";
open my $fh, '<', $file or die "Can't open $file: $!";
my ($seq, $seq_name);
while (<$fh>) {
chomp;
if (/^>(.*)/) {
# Process the previous assembled sequence
if ($seq) {
proc_seq($seq_name, $seq);
$seq = '';
}
$seq_name = $1;
next;
}
$seq .= $_;
}
# Process the last one
proc_seq($seq_name, $seq);
sub proc_seq {
my ($seq_name, $seq, $multiline) = @_;
# Build output in the loop, as motifs are found. By default, print all
# output for one seq_name in one line. To print each motif on its own
# line instead, invoke this sub with a true third argument (1 will do).
my $output = ">$seq_name";
my $cnt = 0;
while ($seq =~ /([HKR])(?=([HKR]{2}))/g) {
++$cnt;
my $motif = $1 . $2;
my $pos = pos($seq);
my $pre_context = ($pos >= 11)
? substr($seq, $pos-11, 10)
: substr($seq, 0, $pos-1);
my $post_context = substr $seq, $pos+2, 10;
$output .= " n$cnt($pos~" . ($pos+2) . ") ";
$output .= "\n" if $multiline;
$output .= lc($pre_context) . $motif . lc($post_context);
}
say ($cnt > 0 ? $output : $output . ' no match found');
}
关于正则表达式的注意事项：我们需要对第二个和第三个字符进行预测，以便能够捕获重叠的图案。
一个例子。有HHKK第一序列，其中包含在重叠HHK和HKK。如果正则表达式匹配HHK使用/[HKR]{3}/之后，那么正则表达式引擎在字符串中的位置将在第一个之后K，因为它“消耗” HHK。所以它接下来只能看到一个K，所以接下来就不会[HKR]{3}匹配，因此它会错过下一个主题。
这就是为什么我只匹配一个字母并为接下来的两个字母做一个“前瞻”：然后在匹配之后H（并且“看到”确实有HK跟随）只消耗了一个字母并且引擎只是先过去了H，它是H在下一场比赛的第二场之前定位。现在它将能够以HKK相同的方式匹配（因此它可以保持匹配甚至多重叠的图案）。
消耗所有三个字符消耗第一个，预测下两个
HHKK HHKK
| |
[HP] {3} | [HP]（= [HP] {2}？）|
| |
旁边的比赛旁边的比赛
这标识了所需输出中指示的所有内容（具有拼写错误）; 注意注释中要求的变化，在一行上打印一个序列的所有图案。所以它打印
> N1 NP_001002156.1（7〜9）mktavdRRKldllysrykd N2（148〜150）lglwntflleHHKksipkdtwnl N3（149〜151）glwntfllehHKKsipkdtwnll
> NP_957070.2 n1（163~165）schvydqkniRRRvydalnvlma
> NP_bogus_with_no_motifs未找到匹配项
在一行中具有相同序列名称的所有图案，如所需。我添加了一条虚假的输入线，没有图案来测试no match found添加; 这画了上面输出的最后一行。
仍然有一个选项可以在一个单独的行上打印每个主题，就像最初想要的那样：调用proc_seq带有附加的第三个参数的函数，这是真的，就像
proc_seq($seq_name, $seq, 1)
然后它会打印出来
> NP_001002156.1 n1（7~9）
mktavdRRldllight富n2（148~150）
dillfHHKksipkdtwnl n3（149~151）
glwntfllehHKKsipkdtwnll
> NP_957070.2 n1（163~165）
schvydqkniRRRvydalnvlma
> NP_bogus_with_no_motifs未找到匹配项
https://stackoverflow.com/questions/-100006596
复制相似问题