如何使用Perl在字符串中搜索生物信息学基序的不同变体?

p5cysglq  于 6个月前  发布在  Perl
关注(0)|答案(3)|浏览(73)

我有一个在不同变体中有一个tandem repeat的程序输出。是否可以搜索(在字符串中)motif并告诉程序找到最多“3”个错配/插入/删除的所有变体?

wnavrhmk

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,或者从头到尾搜索字符串,而不是在固定的偏移量上搜索几个字符串。

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;
    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";
}

字符串
输出量:

1000 random strings generated 10858 characters long
Target:       CGTCGCACAG
Candidate:    CGTCGCACAG
Differences:                   0 differences
Array element: 999 

Target:       CGTCGCACAG
Candidate:    CGTCGCCGCG
Differences:        CGC        3 differences
Array element: 696 

Target:       CGTCGCACAG
Candidate:    CGTCGCCGAT
Differences:        CG T       3 differences
Array element: 851 

Target:       CGTCGCACAG
Candidate:    CGTCGCATGG
Differences:         TG        2 differences
Array element: 986 

Target:       CGTCGCACAG
Candidate:    CATGGCACGG
Differences:   A G    G        3 differences
Array element: 998 

..several cut out.. 

Target:       CGTCGCACAG
Candidate:    CGTCGCTCCA
Differences:        T CA       3 differences
Array element: 568 926

zte4gxcn

zte4gxcn2#

我相信在BioPerl中有这类事情的例程。
无论如何,如果你在BioStar, the bioinformatics stack exchange上问这个问题,你可能会得到更好的答案。

egmofgnx

egmofgnx3#

当我在学习perl的头几年,我写了一个我现在认为是非常低效的(但功能)串联重复序列查找器(以前在我以前工作的公司网站上可以找到)叫做tandyman。几年后我写了一个模糊的版本叫做cottonTandy。如果我今天重写它,我将使用哈希进行全局搜索(给定允许的错误),并利用模式匹配进行局部搜索。
下面是一个如何使用它的示例:

#!/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 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.质量指标为1
1.由空格分隔的重复单元
请注意,很容易提供参数(正如您从上面的输出中看到的那样),这些参数将输出垃圾/无关紧要的“重复”,但是如果您知道如何提供好的参数,它可以找到您设置的内容。
不幸的是,这个软件包还没有公开,我从来没有费心去提供它,因为它太慢了,甚至不适合原核大小的基因组搜索(尽管它对单个基因是可行的)。在我初学编码的日子里,我已经开始添加一个功能,以采取“状态”作为输入,这样我就可以在一个序列的部分上并行运行它,我从来没有完成过,一旦我学会了哈希会使它快得多。到那时,我已经转移到其他项目。但如果它适合你的需要,给我发信息,我可以给你发一份副本。
它只有不到1000行的代码,但它有很多花里胡哨的东西,比如允许IUPAC模糊码(BDHVRYKMSWN)。它对氨基酸和核酸都有效。它过滤掉内部重复(例如,不报告TTTT或ATAT为4 nt共识)。

相关问题