我有一个在不同变体中有一个tandem repeat的程序输出。是否可以搜索（在字符串中）motif并告诉程序找到最多“3”个错配/插入/删除的所有变体？
wnavrhmk1#
下面的程序在一个1,000个元素的数组中生成了一个由2个字符对组成的字符串，长度为5,428对。我意识到你更可能是从文件中阅读这些字符串，但这只是一个例子。显然，你会用任何来源的实际数据替换随机字符串。我不知道我使用的'AT','CG','TC','CA','TG','GC','GG'是否是合法的碱基对组合。（我在生物学上睡过觉......）只要将map块对编辑为合法的对，并将7更改为对的数量，如果你想生成合法的随机字符串进行测试。如果在offset点的子字符串是3个或更少的差异，数组元素（标量值）存储在哈希值部分的匿名数组中。哈希值的关键部分是接近匹配的子字符串。而不是数组元素，值可以是文件名，Perl数据引用或其他相关引用，您希望与motif关联。虽然我刚刚看了字符串之间的字符差异，但您可以通过将foreach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); }行替换为适用于您的问题的比较逻辑来放置任何需要查看的特定逻辑。我不知道mismatches/insertions/deletions如何在字符串中表示，所以我将其作为练习留给读者。很容易修改这个程序，让键盘输入$target和$offset，或者从头到尾搜索字符串，而不是在固定的偏移量上搜索几个字符串。
'AT','CG','TC','CA','TG','GC','GG'
map
7
offset
foreach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); }
mismatches/insertions/deletions
$target
$offset
use strict; use warnings;my @bps;push(@bps,join('',map { ('AT','CG','TC','CA','TG','GC','GG')[rand 7] } 0..5428)) for(1..1_000);my $len=length($bps[0]);my $s_count= scalar @bps;print "$s_count random strings generated $len characters long\n" ;my $target="CGTCGCACAG";my $offset=832;my $nlen=length $target;my %HoA;my $diffs=0;my @a2=split(//, $target);substr($bps[-1], $offset, $nlen)=$target; #guarantee 1 matchsubstr($bps[-2], $offset, $nlen)="CATGGCACGG"; #anja exampleforeach my $i (0..$#bps) { my $cand=substr($bps[$i], $offset, $nlen); my @a1=split(//, $cand); $diffs=0; foreach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); } next if $diffs > 3; push (@{$HoA{$cand}}, $i); }foreach my $hit (keys %HoA) { my @a1=split(//, $hit); $diffs=0; my $ds=""; foreach my $j (0..$#a1) { if($a1[$j] eq $a2[$j]) { $ds.=" "; } else { $diffs++; $ds.=$a1[$j]; } } print "Target: $target\n", "Candidate: $hit\n", "Differences: $ds $diffs differences\n", "Array element: "; foreach (@{$HoA{$hit}}) { print "$_ " ; } print "\n\n";}
use strict; use warnings;
my @bps;
push(@bps,join('',map { ('AT','CG','TC','CA','TG','GC','GG')[rand 7] }
0..5428)) for(1..1_000);
my $len=length($bps[0]);
my $s_count= scalar @bps;
print "$s_count random strings generated $len characters long\n" ;
my $target="CGTCGCACAG";
my $offset=832;
my $nlen=length $target;
my %HoA;
my $diffs=0;
my @a2=split(//, $target);
substr($bps[-1], $offset, $nlen)=$target; #guarantee 1 match
substr($bps[-2], $offset, $nlen)="CATGGCACGG"; #anja example
foreach my $i (0..$#bps) {
my $cand=substr($bps[$i], $offset, $nlen);
my @a1=split(//, $cand);
$diffs=0;
next if $diffs > 3;
push (@{$HoA{$cand}}, $i);
}
foreach my $hit (keys %HoA) {
my @a1=split(//, $hit);
my $ds="";
foreach my $j (0..$#a1) {
if($a1[$j] eq $a2[$j]) {
$ds.=" ";
} else {
$diffs++;
$ds.=$a1[$j];
print "Target: $target\n",
"Candidate: $hit\n",
"Differences: $ds $diffs differences\n",
"Array element: ";
foreach (@{$HoA{$hit}}) {
print "$_ " ;
print "\n\n";
字符串输出量：
1000 random strings generated 10858 characters longTarget: CGTCGCACAGCandidate: CGTCGCACAGDifferences: 0 differencesArray element: 999 Target: CGTCGCACAGCandidate: CGTCGCCGCGDifferences: CGC 3 differencesArray element: 696 Target: CGTCGCACAGCandidate: CGTCGCCGATDifferences: CG T 3 differencesArray element: 851 Target: CGTCGCACAGCandidate: CGTCGCATGGDifferences: TG 2 differencesArray element: 986 Target: CGTCGCACAGCandidate: CATGGCACGGDifferences: A G G 3 differencesArray element: 998 ..several cut out.. Target: CGTCGCACAGCandidate: CGTCGCTCCADifferences: T CA 3 differencesArray element: 568 926
1000 random strings generated 10858 characters long
Target: CGTCGCACAG
Candidate: CGTCGCACAG
Differences: 0 differences
Array element: 999
Candidate: CGTCGCCGCG
Differences: CGC 3 differences
Array element: 696
Candidate: CGTCGCCGAT
Differences: CG T 3 differences
Array element: 851
Candidate: CGTCGCATGG
Differences: TG 2 differences
Array element: 986
Candidate: CATGGCACGG
Differences: A G G 3 differences
Array element: 998
..several cut out..
Candidate: CGTCGCTCCA
Differences: T CA 3 differences
Array element: 568 926
型
zte4gxcn2#
我相信在BioPerl中有这类事情的例程。无论如何，如果你在BioStar, the bioinformatics stack exchange上问这个问题，你可能会得到更好的答案。
egmofgnx3#
当我在学习perl的头几年，我写了一个我现在认为是非常低效的（但功能）串联重复序列查找器（以前在我以前工作的公司网站上可以找到）叫做tandyman。几年后我写了一个模糊的版本叫做cottonTandy。如果我今天重写它，我将使用哈希进行全局搜索（给定允许的错误），并利用模式匹配进行局部搜索。下面是一个如何使用它的示例：
#!/usr/bin/perluse Tandyman;$sequence = "ATGCATCGTAGCGTTCAGTCGGCATCTATCTGACGTACTCTTACTGCATGAGTCTAGCTGTACTACGTACGAGCTGAGCAGCGTACgTG";my $tandy = Tandyman->new(\$sequence,'n'); #Can't believe I coded it to take a scalar reference! Prob. fresh out of a cpp class when I wrote it.$tandy->SetParams(4,2,3,3,4);#The parameters are, in order:# repeat unit size# min number of repeat units to require a hit# allowed mistakes per unit (an upper bound for "mistake concentration")# allowed mistakes per window (a lower bound for "mistake concentration")# number of units in a "window"while(@repeat_info = $tandy->FindRepeat()) {print(join("\t",@repeat_info),"\n")}
#!/usr/bin/perl
use Tandyman;
$sequence = "ATGCATCGTAGCGTTCAGTCGGCATCTATCTGACGTACTCTTACTGCATGAGTCTAGCTGTACTACGTACGAGCTGAGCAGCGTACgTG";
my $tandy = Tandyman->new(\$sequence,'n'); #Can't believe I coded it to take a scalar reference! Prob. fresh out of a cpp class when I wrote it.
$tandy->SetParams(4,2,3,3,4);
#The parameters are, in order:
# repeat unit size
# min number of repeat units to require a hit
# allowed mistakes per unit (an upper bound for "mistake concentration")
# allowed mistakes per window (a lower bound for "mistake concentration")
# number of units in a "window"
while(@repeat_info = $tandy->FindRepeat())
{print(join("\t",@repeat_info),"\n")}
字符串这个测试的输出看起来像这样（运行时间长达11秒）：
25 32 TCTA 2 0.87 TCTA TCTG58 72 CGTA 4 0.81 CTGTA CTA CGTA CGA82 89 CGTA 2 0.87 CGTA CGTG45 51 TGCA 2 0.87 TGCA TGA65 72 ACGA 2 0.87 ACGT ACGA23 29 CTAT 2 0.87 CAT CTAT36 45 TACT 3 0.83 TACT CT TACT24 31 ATCT 2 1 ATCT ATCT51 59 AGCT 2 0.87 AGTCT AGCT33 39 ACGT 2 0.87 ACGT ACT62 72 ACGT 3 0.83 ACT ACGT ACGA80 88 ACGT 2 0.87 AGCGT ACGT81 88 GCGT 2 0.87 GCGT ACGT63 70 CTAC 2 0.87 CTAC GTAC32 38 GTAC 2 0.87 GAC GTAC60 74 GTAC 4 0.81 GTAC TAC GTAC GAGC23 30 CATC 2 0.87 CATC TATC71 82 GAGC 3 0.83 GAGC TGAGC AGC1 7 ATGC 2 0.87 ATGC ATC54 60 CTAG 2 0.87 CTAG CTG15 22 TCAG 2 0.87 TCAG TCGG70 81 CGAG 3 0.83 CGAG CTGAG CAG44 50 CATG 2 0.87 CTG CATG25 32 TCTG 2 0.87 TCTA TCTG82 89 CGTG 2 0.87 CGTA CGTG55 73 TACG 5 0.75 TAGCTG TAC TACG TACG AG69 83 AGCG 4 0.81 ACG AGCTG AGC AGCG15 22 TCGG 2 0.87 TCAG TCGG
25 32 TCTA 2 0.87 TCTA TCTG
58 72 CGTA 4 0.81 CTGTA CTA CGTA CGA
82 89 CGTA 2 0.87 CGTA CGTG
45 51 TGCA 2 0.87 TGCA TGA
65 72 ACGA 2 0.87 ACGT ACGA
23 29 CTAT 2 0.87 CAT CTAT
36 45 TACT 3 0.83 TACT CT TACT
24 31 ATCT 2 1 ATCT ATCT
51 59 AGCT 2 0.87 AGTCT AGCT
33 39 ACGT 2 0.87 ACGT ACT
62 72 ACGT 3 0.83 ACT ACGT ACGA
80 88 ACGT 2 0.87 AGCGT ACGT
81 88 GCGT 2 0.87 GCGT ACGT
63 70 CTAC 2 0.87 CTAC GTAC
32 38 GTAC 2 0.87 GAC GTAC
60 74 GTAC 4 0.81 GTAC TAC GTAC GAGC
23 30 CATC 2 0.87 CATC TATC
71 82 GAGC 3 0.83 GAGC TGAGC AGC
1 7 ATGC 2 0.87 ATGC ATC
54 60 CTAG 2 0.87 CTAG CTG
15 22 TCAG 2 0.87 TCAG TCGG
70 81 CGAG 3 0.83 CGAG CTGAG CAG
44 50 CATG 2 0.87 CTG CATG
25 32 TCTG 2 0.87 TCTA TCTG
82 89 CGTG 2 0.87 CGTA CGTG
55 73 TACG 5 0.75 TAGCTG TAC TACG TACG AG
69 83 AGCG 4 0.81 ACG AGCTG AGC AGCG
15 22 TCGG 2 0.87 TCAG TCGG
型正如你所看到的，它允许indel和SNP。列的顺序是：1.开始位置1.停止位置1.共有序列1.找到的单位数1.质量指标为11.由空格分隔的重复单元请注意，很容易提供参数（正如您从上面的输出中看到的那样），这些参数将输出垃圾/无关紧要的“重复”，但是如果您知道如何提供好的参数，它可以找到您设置的内容。不幸的是，这个软件包还没有公开，我从来没有费心去提供它，因为它太慢了，甚至不适合原核大小的基因组搜索（尽管它对单个基因是可行的）。在我初学编码的日子里，我已经开始添加一个功能，以采取“状态”作为输入，这样我就可以在一个序列的部分上并行运行它，我从来没有完成过，一旦我学会了哈希会使它快得多。到那时，我已经转移到其他项目。但如果它适合你的需要，给我发信息，我可以给你发一份副本。它只有不到1000行的代码，但它有很多花里胡哨的东西，比如允许IUPAC模糊码（BDHVRYKMSWN）。它对氨基酸和核酸都有效。它过滤掉内部重复（例如，不报告TTTT或ATAT为4 nt共识）。
3条答案
按热度按时间wnavrhmk1#
下面的程序在一个1,000个元素的数组中生成了一个由2个字符对组成的字符串，长度为5,428对。我意识到你更可能是从文件中阅读这些字符串，但这只是一个例子。显然，你会用任何来源的实际数据替换随机字符串。
我不知道我使用的
'AT','CG','TC','CA','TG','GC','GG'是否是合法的碱基对组合。（我在生物学上睡过觉......）只要将map块对编辑为合法的对，并将7更改为对的数量，如果你想生成合法的随机字符串进行测试。如果在
offset点的子字符串是3个或更少的差异，数组元素（标量值）存储在哈希值部分的匿名数组中。哈希值的关键部分是接近匹配的子字符串。而不是数组元素，值可以是文件名，Perl数据引用或其他相关引用，您希望与motif关联。虽然我刚刚看了字符串之间的字符差异，但您可以通过将
foreach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); }行替换为适用于您的问题的比较逻辑来放置任何需要查看的特定逻辑。我不知道mismatches/insertions/deletions如何在字符串中表示，所以我将其作为练习留给读者。很容易修改这个程序，让键盘输入
$target和$offset，或者从头到尾搜索字符串，而不是在固定的偏移量上搜索几个字符串。字符串
输出量：
型
zte4gxcn2#
我相信在BioPerl中有这类事情的例程。
无论如何，如果你在BioStar, the bioinformatics stack exchange上问这个问题，你可能会得到更好的答案。
egmofgnx3#
当我在学习perl的头几年，我写了一个我现在认为是非常低效的（但功能）串联重复序列查找器（以前在我以前工作的公司网站上可以找到）叫做tandyman。几年后我写了一个模糊的版本叫做cottonTandy。如果我今天重写它，我将使用哈希进行全局搜索（给定允许的错误），并利用模式匹配进行局部搜索。
下面是一个如何使用它的示例：
字符串
这个测试的输出看起来像这样（运行时间长达11秒）：
型
正如你所看到的，它允许indel和SNP。列的顺序是：
1.开始位置
1.停止位置
1.共有序列
1.找到的单位数
1.质量指标为1
1.由空格分隔的重复单元
请注意，很容易提供参数（正如您从上面的输出中看到的那样），这些参数将输出垃圾/无关紧要的“重复”，但是如果您知道如何提供好的参数，它可以找到您设置的内容。
不幸的是，这个软件包还没有公开，我从来没有费心去提供它，因为它太慢了，甚至不适合原核大小的基因组搜索（尽管它对单个基因是可行的）。在我初学编码的日子里，我已经开始添加一个功能，以采取“状态”作为输入，这样我就可以在一个序列的部分上并行运行它，我从来没有完成过，一旦我学会了哈希会使它快得多。到那时，我已经转移到其他项目。但如果它适合你的需要，给我发信息，我可以给你发一份副本。
它只有不到1000行的代码，但它有很多花里胡哨的东西，比如允许IUPAC模糊码（BDHVRYKMSWN）。它对氨基酸和核酸都有效。它过滤掉内部重复（例如，不报告TTTT或ATAT为4 nt共识）。