#!/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")}
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共识)。